Dynamic assessment of balance ability based on time-series data and fall prevention training method

Abstract

The application provides a balance ability dynamic evaluation and fall prevention training method based on time sequence data, relates to the technical field of control theory, and comprises the following steps: collecting user foot pressure time sequence data and joint motion data synchronously, constructing a coupling correlation matrix of pressure center trajectory space shape features and joint coordination features, and extracting a defect feature vector; constructing a multi-objective optimization function, solving a Pareto optimal solution set, extracting balance defect sensitivity feature parameters, performing hierarchical quantitative evaluation, and generating a training demand evaluation report. The application can realize accurate evaluation of user balance ability and personalized training scheme design, and improve the pertinence and effectiveness of fall prevention.

Description

Technical Field

[0001] This invention relates to control theory technology, and more particularly to a method for dynamic assessment of balance ability and fall prevention training based on time-series data. Background Technology

[0002] Falls are a major health risk for the elderly and those with balance disorders, and prevention requires accurate assessment of individual dynamic balance abilities and targeted training programs. Existing balance assessment methods mainly rely on static posture tests or simple functional scales, which fail to capture the complexity of the body's balance regulation mechanisms during dynamic movement and individual variability. Current assessment systems typically employ fixed assessment criteria and recommended training strategies, failing to fully consider the significant differences among individuals in the type of balance deficit, its sensitivity, and adaptability, resulting in insufficient targeting and effectiveness of training. Summary of the Invention

[0003] This invention provides a method for dynamic assessment of balance ability and fall prevention training based on time-series data, which can solve the problems in the prior art.

[0004] A first aspect of the present invention provides a method for dynamic assessment of balance ability and fall prevention training based on time-series data, comprising:

[0005] Simultaneously collect foot pressure time-series data and joint motion data of users on the mobile training platform;

[0006] The pressure center trajectory in the foot pressure time series data is divided into continuous observation periods according to the time dimension. For each observation period, a coupling correlation matrix is ​​constructed between the spatial morphological features of the pressure center trajectory and the joint coordination features of the joint motion data. The coupling correlation matrix is ​​then decomposed by eigenvalue to extract defect feature vectors.

[0007] Construct a multi-objective optimization function, which includes an optimization objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, and solve the Pareto optimal solution set of the multi-objective optimization function;

[0008] Feature parameters characterizing the user's sensitivity to balance defects are extracted from the Pareto optimal solution set. Based on the feature parameters, the user's dynamic balance ability is evaluated in a graded and quantitative manner, and a balance training requirement evaluation report is generated, which includes defect type identification results, sensitivity level classification, and corresponding training suggestion parameters.

[0009] The steps of dividing the pressure center trajectory in the foot pressure time series data into continuous observation periods according to the time dimension, constructing a coupling correlation matrix between the spatial morphological features of the pressure center trajectory and the joint coordination features of the joint motion data for each observation period, and extracting defect feature vectors by eigenvalue decomposition of the coupling correlation matrix include:

[0010] The pressure center trajectory is divided into continuous observation periods according to a preset time window length. For the pressure center trajectory within each observation period, the rate of change of movement direction and the rate of change of curvature are calculated. The rate of change of movement direction and the rate of change of curvature are decomposed in time and frequency, and the frequency domain energy distribution vector is extracted as a spatial morphological feature.

[0011] For the joint motion data corresponding to the observation period, the angular velocity phase difference between adjacent joints is extracted, the angular velocity phase difference is continuously sampled, and the phase synchronization coefficient sequence is calculated as a joint coordination feature;

[0012] A coupling correlation matrix is ​​constructed using the spatial morphological features as row vectors and the joint coordination features as column vectors; the coupling correlation matrix is ​​decomposed into eigenvalues, and the eigenvector corresponding to the largest eigenvalue is extracted as the defect feature vector. The weight distribution of each element in the defect feature vector reflects the degree of influence of the spatial morphological features on the joint coordination features.

[0013] The steps of performing time-frequency decomposition on the rate of change of the direction of movement and the rate of change of curvature to extract the frequency domain energy distribution vector include:

[0014] Wavelet packet decomposition is performed on the rate of change of the direction of movement and the rate of change of curvature respectively, and the frequency band is divided into three sub-bands: low-frequency steady-state component, mid-frequency transition component and high-frequency abrupt component. The energy proportion of each sub-band is calculated to form an initial frequency domain energy distribution vector.

[0015] The trajectory smoothness is evaluated by calculating the autocorrelation coefficient of the low-frequency steady-state component. The autocorrelation coefficient is then weighted and fused with the initial frequency domain energy distribution vector. When the autocorrelation coefficient is lower than a preset smoothing threshold, the energy weight of the high-frequency abrupt component is increased. When the autocorrelation coefficient is higher than the preset smoothing threshold, the energy weight of the low-frequency steady-state component is increased, thereby generating a corrected frequency domain energy distribution vector.

[0016] For the joint motion data corresponding to the observation period, the steps of extracting the angular velocity phase difference between adjacent joints, continuously sampling the angular velocity phase difference, and calculating the phase synchronization coefficient sequence as a joint coordination feature include:

[0017] The angular velocity time-series data of the joint nodes in the lower limb kinetic chain are extracted from the joint motion data, and the angular velocity phase difference between adjacent joint nodes is calculated; the angular velocity phase difference is subjected to analytical signal transformation to extract the instantaneous phase sequence, and the phase synchronization coefficient is calculated based on the time-domain fluctuation characteristics of the instantaneous phase sequence;

[0018] During the observation period, the phase synchronization coefficients are continuously sampled to form a phase synchronization coefficient sequence; the switching frequency between different coordination states in the phase synchronization coefficient sequence is statistically analyzed, a state transition frequency matrix is ​​constructed and normalized, the weight distribution of the dominant transition path representing the stable coordination mode in the normalized transition frequency matrix is ​​extracted, and the weight distribution is fused with the phase synchronization coefficient sequence to generate an enhanced joint coordination feature vector.

[0019] Constructing a multi-objective optimization function, wherein the multi-objective optimization function includes an optimization objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, and the steps of solving the Pareto optimal solution set of the multi-objective optimization function include:

[0020] A multi-objective optimization function is constructed, which includes a first sub-objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, a second sub-objective of minimizing the probability that the pressure center trajectory offset exceeds the safe range, and a third sub-objective of minimizing the fluctuation amplitude of the joint coordination feature.

[0021] During the continuous acquisition of the foot pressure time series data, the time-domain sequence of the offset of the pressure center trajectory relative to the preset reference area is calculated, and the peak fluctuation amplitude and steady-state residual of the offset are extracted as adaptation evaluation indicators through sliding window statistical analysis.

[0022] The weight coefficients among the first, second, and third sub-objectives are dynamically adjusted based on the deviation between the adaptation evaluation index and the preset adaptation threshold. When the adaptation evaluation index is lower than the preset adaptation threshold, the weight coefficient of the first sub-objective is increased; when the adaptation evaluation index is higher than the preset adaptation threshold, the weight coefficient of the first sub-objective is decreased. The Pareto optimal solution set of the multi-objective optimization function is then solved under the updated weight coefficients.

[0023] The steps of dynamically adjusting the weight coefficients among the first, second, and third sub-objectives based on the deviation between the adaptation evaluation index and the preset adaptation threshold, and solving for the Pareto optimal solution set of the multi-objective optimization function under the updated weight coefficients include:

[0024] Calculate the relative deviation rate between the adaptation evaluation index and the preset adaptation threshold and perform normalization processing; set the initial weight coefficients of the first sub-objective, the second sub-objective and the third sub-objective, wherein the initial weight coefficient of the first sub-objective is higher than that of the other sub-objectives;

[0025] The weight adjustment factor is calculated using a piecewise nonlinear function based on the normalized relative deviation rate. The slope of the adjustment factor's response to the deviation rate within each interval of the piecewise nonlinear function satisfies an increasing relationship. The weight coefficients of each sub-objective are updated based on the weight adjustment factor while maintaining the weight normalization constraint.

[0026] A multi-objective evolutionary optimization algorithm is used to solve the multi-objective optimization function after updating the weights. A search direction guidance vector is constructed based on the weight distribution characteristics of the defect feature vector and used as the directional perturbation term for the population mutation operation. The proportion of individuals in the current population whose pressure center trajectory offset exceeds the safe range is statistically analyzed, and the relationship with the preset proportion threshold is used to adjust the population size and mutation probability in reverse. Individuals in the population are screened and retained according to the target space crowding distance.

[0027] The steps of extracting feature parameters characterizing the user's sensitivity to balance defects from the Pareto optimal solution set, performing a graded and quantitative assessment of the user's dynamic balance ability based on the feature parameters, and generating a balance training requirement assessment report containing defect type identification results, sensitivity level classification, and corresponding training suggestion parameters include:

[0028] The ratio of the response magnitude of the defect feature vector of each solution in the Pareto optimal solution set to the degree of constraint violation is calculated as the defect sensitivity index, and the sensitivity level classification is determined based on the distribution of the defect sensitivity index.

[0029] From the Pareto optimal solution set, candidate safe solutions with constraint violation degree below the preset safety threshold are selected. The solution with the largest response amplitude of the defect feature vector in the candidate safe solution set is selected as the optimal evaluation solution. The defect feature vector weight distribution pattern of the optimal evaluation solution is matched with the preset defect pattern library for similarity, and the defect type identifier with the highest similarity is output as the defect type identification result.

[0030] Based on the defect type identification result, the corresponding benchmark training parameters are queried. The benchmark training parameters include the training intensity range and training mode category. The benchmark training parameters are adaptively adjusted according to the sensitivity level classification, and the adjusted parameters are used as training suggestion parameters.

[0031] The defect type identification results, the sensitivity level classification, and the training suggestion parameters are organized to generate a balanced training requirement assessment report.

[0032] A second aspect of the present invention provides an electronic device, comprising:

[0033] processor;

[0034] Memory used to store processor-executable instructions;

[0035] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0036] A third aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0037] This invention achieves dynamic, real-time assessment of a user's balance ability by simultaneously acquiring foot pressure time-series data and joint motion data. Compared to traditional static assessment methods, it can more comprehensively reflect a user's balance control ability under actual movement conditions. Based on a multi-objective optimization function and Pareto optimal solution set, it quantitatively characterizes the sensitivity of the user's balance defects, realizing a graded and quantitative assessment of balance ability, making the assessment results more objective and comparable. Attached Figure Description

[0038] Figure 1 This is a flowchart illustrating the dynamic assessment of balance ability and fall prevention training method based on time-series data, as described in an embodiment of the present invention.

[0039] Figure 2 Flowchart for constructing the coupling correlation matrix and extracting defect features. Detailed Implementation

[0040] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0041] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0042] Figure 1 This is a flowchart illustrating the dynamic assessment of balance ability and fall prevention training method based on time-series data according to an embodiment of the present invention. Figure 1 As shown, the method includes:

[0043] Simultaneously collect foot pressure time-series data and joint motion data of users on the mobile training platform;

[0044] The pressure center trajectory in the foot pressure time series data is divided into continuous observation periods according to the time dimension. For each observation period, a coupling correlation matrix is ​​constructed between the spatial morphological features of the pressure center trajectory and the joint coordination features of the joint motion data. The coupling correlation matrix is ​​then decomposed by eigenvalue to extract defect feature vectors.

[0045] Construct a multi-objective optimization function, which includes an optimization objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, and solve the Pareto optimal solution set of the multi-objective optimization function;

[0046] Feature parameters characterizing the user's sensitivity to balance defects are extracted from the Pareto optimal solution set. Based on the feature parameters, the user's dynamic balance ability is evaluated in a graded and quantitative manner, and a balance training requirement evaluation report is generated, which includes defect type identification results, sensitivity level classification, and corresponding training suggestion parameters.

[0047] In one optional implementation, the pressure center trajectory in the foot pressure time series data is divided into continuous observation periods according to the time dimension. For each observation period, a coupling correlation matrix is ​​constructed between the spatial morphological features of the pressure center trajectory and the joint coordination features of the joint motion data. The step of extracting defect feature vectors by eigenvalue decomposition of the coupling correlation matrix includes:

[0048] The pressure center trajectory is divided into continuous observation periods according to a preset time window length. For the pressure center trajectory within each observation period, the rate of change of movement direction and the rate of change of curvature are calculated. The rate of change of movement direction and the rate of change of curvature are decomposed in time and frequency, and the frequency domain energy distribution vector is extracted as a spatial morphological feature.

[0049] For the joint motion data corresponding to the observation period, the angular velocity phase difference between adjacent joints is extracted, the angular velocity phase difference is continuously sampled, and the phase synchronization coefficient sequence is calculated as a joint coordination feature;

[0050] A coupling correlation matrix is ​​constructed using the spatial morphological features as row vectors and the joint coordination features as column vectors; the coupling correlation matrix is ​​decomposed into eigenvalues, and the eigenvector corresponding to the largest eigenvalue is extracted as the defect feature vector. The weight distribution of each element in the defect feature vector reflects the degree of influence of the spatial morphological features on the joint coordination features.

[0051] Combination Figure 2The flowchart for constructing the coupling correlation matrix and extracting defect features is illustrated below. For example, when a user performs balance training on a mobile training platform, foot pressure distribution data is collected in real time through a foot pressure sensor array embedded in the platform. Simultaneously, a motion capture device based on an inertial measurement unit or optical markers acquires the three-dimensional spatial position and posture information of the hip, knee, and ankle joints of the lower limbs, and calculates kinematic parameters such as angles and angular velocities of each joint.

[0052] The pressure center trajectory is extracted from foot pressure distribution data. Specifically, at each sampling moment, the pressure values ​​of all sensing units on the sensor array are used as weights to calculate the weighted centroid coordinates of the spatial positions of these units. This centroid is the pressure center position at the current moment. The two-dimensional planar trajectory line formed by the continuous change of the pressure center position over time is the pressure center trajectory, which reflects the dynamic projection of the body's center of gravity onto the support surface. Under normal gait, the pressure center trajectory exhibits a regular forward and backward and left and right swing pattern, while the trajectory of individuals with impaired balance shows irregular swinging, excessive deviation, or frequent abrupt changes in direction.

[0053] The acquired pressure center trajectory was divided into multiple continuous observation periods along the time dimension. The initial time window length was set to 1 second, and a 50% overlap sliding strategy was used between windows to ensure temporal continuity. At the same time, the user's step frequency was monitored. When the step frequency was lower than 90 steps per minute, the window length was automatically adjusted to 1.5 seconds to ensure that each window contained at least one complete gait cycle, providing a sufficient data foundation for subsequent feature extraction.

[0054] For the trajectory of the pressure center within each observation period, two geometric dynamic parameters are calculated: the rate of change of its direction of movement and the rate of change of its curvature. The rate of change of direction of movement describes the speed of change in the trajectory's direction. During calculation, sampling points are extracted on the trajectory at 10-millisecond intervals. The angle between two vector segments formed by three consecutive points is calculated, and the change in angle is divided by the corresponding time interval to obtain the time series of the rate of change of direction. The rate of change of curvature characterizes the evolution of the trajectory's curvature. A sliding three-point circle fitting method is used: for any point on the trajectory and its adjacent points, a unique fitting circle is determined using three points. The reciprocal of the radius of this circle is the curvature of the trajectory at that point. The rate of change of curvature is obtained by performing time difference calculation on the curvature sequence.

[0055] The time series of the rate of change of movement direction and the rate of change of curvature were subjected to time-frequency domain transformation. Short-time Fourier transform was used to decompose the time-domain signal into a superposition of different frequency components. A Hanning window was chosen to reduce spectral leakage, with a window length of 128 sampling points and an overlap of 64 sampling points between adjacent windows. The transformed spectrum was divided into three frequency bands according to physiological movement characteristics: the low-frequency band (0-2 Hz) corresponds to the basic gait rhythm, the mid-frequency band (2-5 Hz) corresponds to posture adjustment and gait transition processes, and the high-frequency band (5-10 Hz) corresponds to fine neuromuscular control responses. The spectral energy was integrated within each frequency band to obtain the energy value for each band. The energy values ​​of the three bands were normalized to form an energy proportion vector. Since the rate of change of movement direction and the rate of change of curvature each generate a three-dimensional energy proportion vector, these two vectors were concatenated to obtain a six-dimensional spatial morphological feature vector.

[0056] Joint coordination features are extracted from joint motion data corresponding to the observation period, obtaining angular velocity time series for the hip, knee, and ankle joints. These angular velocity signals reflect the rotational speed of each joint. A Hilbert transform is applied to the angular velocity signal of each joint, which expands the real signal into a complex analytic signal, with the complex-valued argument representing the instantaneous phase. By extracting the instantaneous phase series of each joint's angular velocity, the phase difference between adjacent joint pairs is calculated. Specifically, the phase difference time series for the hip and knee joints, and the knee and ankle joints, are calculated. The phase difference values ​​reflect the temporal lead or lag relationship between the two joint movements.

[0057] Phase synchronization coefficients are calculated based on phase difference sequences to quantify the degree of coordination between joints. A phase-locking value method is used, mapping the phase difference value at each moment in the phase difference sequence to a unit vector on the complex plane. The average of these unit vectors over the observation period is calculated, and the magnitude of this average vector is the phase-locking value. A phase-locking value close to 1 indicates highly synchronized and coordinated movement between the two joints, while a value close to 0 indicates chaotic and uncoordinated movement. The phase-locking values ​​for the hip-knee joint pair and the knee-ankle joint pair are calculated separately, and the two values ​​are combined to form a two-dimensional joint coordination feature vector.

[0058] A coupling correlation matrix is ​​constructed using a six-dimensional spatial morphological feature vector and a two-dimensional joint coordination feature vector. The spatial morphological feature vector is transposed into a six-row, one-column column vector, while the joint coordination feature vector is preserved as a one-row, two-column row vector. Matrix multiplication is then performed between the two to obtain the six-row, two-column coupling correlation matrix. The element in the i-th row and j-th column of this matrix represents the correlation contribution of the i-th spatial morphological feature component to the j-th joint coordination feature component. The matrix as a whole encodes the coupling relationship between the plantar pressure center motion characteristics and the lower limb joint coordination patterns.

[0059] Singular value decomposition (SVD) is performed on the coupling correlation matrix. The decomposition yields three matrices: a left singular vector matrix, a singular value diagonal matrix, and a right singular vector matrix. The singular values ​​are arranged in descending order, with the largest singular value corresponding to the principal component direction with the highest energy concentration in the matrix. The left singular vector corresponding to the largest singular value is extracted; this six-dimensional vector is defined as the defect feature vector. Each element of the defect feature vector corresponds to one of the six spatial morphological feature components, and the absolute value of each element reflects its weight in influencing joint coordination. If the weight of the fifth element is significantly higher than the others, and the fifth element corresponds to the rate of change of movement direction in the high-frequency band, it indicates that the drastic directional change of the pressure center trajectory in the high-frequency range is the dominant factor leading to decreased joint coordination. This suggests that the user may lack fine control during rapid posture adjustments, and subsequent training should focus on improving neural reaction speed and dynamic balance adjustment capabilities.

[0060] This invention constructs a coupled correlation matrix between the spatial morphological features of the pressure center trajectory and the joint coordination features, and extracts defect feature vectors. This enables the accurate identification of the dominant factors affecting a user's balance ability, providing a scientific basis for the development of personalized training programs.

[0061] In one optional implementation, the step of performing time-frequency decomposition on the rate of change of the direction of movement and the rate of change of curvature to extract the frequency domain energy distribution vector includes:

[0062] Wavelet packet decomposition is performed on the rate of change of the direction of movement and the rate of change of curvature respectively, and the frequency band is divided into three sub-bands: low-frequency steady-state component, mid-frequency transition component and high-frequency abrupt component. The energy proportion of each sub-band is calculated to form an initial frequency domain energy distribution vector.

[0063] The trajectory smoothness is evaluated by calculating the autocorrelation coefficient of the low-frequency steady-state component. The autocorrelation coefficient is then weighted and fused with the initial frequency domain energy distribution vector. When the autocorrelation coefficient is lower than a preset smoothing threshold, the energy weight of the high-frequency abrupt component is increased. When the autocorrelation coefficient is higher than the preset smoothing threshold, the energy weight of the low-frequency steady-state component is increased, thereby generating a corrected frequency domain energy distribution vector.

[0064] For example, after obtaining two time series, the rate of change of movement direction and the rate of change of curvature, wavelet packet decomposition is performed on them respectively to extract the frequency domain energy distribution features. The wavelet packet decomposition uses the Daubechies wavelet basis function, with a decomposition level of three, recursively dividing the spectral space into eight sub-bands. Based on the physiological frequency characteristics of human gait movement, these eight sub-bands are recombined into three physiologically significant frequency bands: the low-frequency steady-state component corresponds to the 0-2 Hz range, covering the basic gait cycle and the slow adjustment process of posture maintenance; the mid-frequency transition component corresponds to the 2-5 Hz range, reflecting gait transition, direction change, and dynamic balance recovery processes; and the high-frequency abrupt change component corresponds to the 5-10 Hz range, reflecting neural reflex regulation and rapid posture correction actions. During the wavelet packet decomposition process, the original signal is decomposed using a filter bank. Each level of decomposition uses a low-pass filter and a high-pass filter to divide the signal into approximate coefficients and detail coefficients, recursively until the preset number of levels is reached. After decomposition, the wavelet packet coefficient sequence of each sub-band is extracted, and the sum of squares of the coefficients within each sub-band is calculated as the energy value of that frequency band. The low-frequency energy is obtained by summing the energies of the subbands corresponding to the low-frequency steady-state components, the mid-frequency energy by summing the energies of the subbands corresponding to the mid-frequency transition components, and the high-frequency energy by summing the energies of the subbands corresponding to the high-frequency abrupt components. The energy values ​​of the three frequency bands are normalized so that their sum is 1, resulting in a three-dimensional initial frequency domain energy distribution vector. The three components represent the proportions of low-frequency, mid-frequency, and high-frequency energy to the total energy. Since the rate of change of the direction of movement and the rate of change of curvature each generate a three-dimensional energy distribution vector, a total of two three-dimensional vectors are obtained.

[0065] To further enhance the sensitivity of the energy distribution vector to gait anomalies, a trajectory smoothness assessment mechanism based on autocorrelation coefficients is introduced to correct the energy distribution. Autocorrelation characteristics are extracted from the wavelet packet coefficient sequence of the low-frequency steady-state component, and the first-order autocorrelation coefficient of this sequence is calculated. Specifically, the sequence is cross-correlated with its sequence delayed by one sampling point, and the autocorrelation coefficient is obtained after normalization. An autocorrelation coefficient close to 1 indicates that the low-frequency component is highly smooth and has strong periodicity, reflecting stable gait rhythm; an autocorrelation coefficient close to zero or negative indicates irregular fluctuations in the low-frequency component, reflecting disordered basic gait patterns. A preset smoothness threshold of 0.7 is set. This threshold is determined based on statistical analysis of gait data from normal adults and represents the typical autocorrelation level of the low-frequency component of a healthy gait. When the calculated autocorrelation coefficient is below 0.7, the trajectory smoothness is deemed insufficient, indicating a deficiency in basic gait control. In this case, the energy weight of the high-frequency abrupt component should be increased to highlight the impact of the lack of fine control ability. When the autocorrelation coefficient is above 0.7, the basic gait rhythm is deemed normal, and the energy weight of the low-frequency steady-state component should be increased to strengthen the expression of periodic characteristics.

[0066] Weight adjustment is achieved through a weighted fusion mechanism. The adjustment factor is defined as the difference between the autocorrelation coefficient and the smoothing threshold. When the difference is negative, its absolute value is multiplied by 1.5 as the high-frequency weight gain coefficient. The high-frequency components in the initial frequency domain energy distribution vector are multiplied by this gain coefficient, while the low-frequency and mid-frequency components are multiplied by 0.8 for suppression. When the difference is positive, its value is multiplied by 1.2 as the low-frequency weight gain coefficient. The low-frequency components are multiplied by this gain coefficient, while the mid-frequency and high-frequency components are multiplied by 0.85 for suppression. After adjustment, the three components are renormalized so that their sum is restored to 1, generating the corrected frequency domain energy distribution vector. This correction mechanism enables the energy distribution characteristics to adaptively reflect different types of gait abnormality patterns: for groups with impaired low-frequency rhythms, such as Parkinson's disease patients, the high-frequency components of the corrected vector are significantly enhanced; for individuals with overly rigid gait lacking adaptive regulation, the enhanced weight of the low-frequency components highlights their lack of dynamic regulation.

[0067] This invention, through wavelet packet decomposition and autocorrelation coefficient correction mechanism, can adaptively enhance the ability of frequency domain energy distribution to distinguish different types of gait anomalies, and improve the accuracy of balance defect feature extraction.

[0068] In one optional implementation, the steps of extracting the angular velocity phase difference between adjacent joints from the joint motion data corresponding to the observation period, continuously sampling the angular velocity phase difference, and calculating a phase synchronization coefficient sequence as a joint coordination feature include:

[0069] The angular velocity time-series data of the joint nodes in the lower limb kinetic chain are extracted from the joint motion data, and the angular velocity phase difference between adjacent joint nodes is calculated; the angular velocity phase difference is subjected to analytical signal transformation to extract the instantaneous phase sequence, and the phase synchronization coefficient is calculated based on the time-domain fluctuation characteristics of the instantaneous phase sequence;

[0070] During the observation period, the phase synchronization coefficients are continuously sampled to form a phase synchronization coefficient sequence; the switching frequency between different coordination states in the phase synchronization coefficient sequence is statistically analyzed, a state transition frequency matrix is ​​constructed and normalized, the weight distribution of the dominant transition path representing the stable coordination mode in the normalized transition frequency matrix is ​​extracted, and the weight distribution is fused with the phase synchronization coefficient sequence to generate an enhanced joint coordination feature vector.

[0071] For example, the lower limb kinetic chain includes three main nodes: the hip, knee, and ankle joints. A motion capture device records the flexion-extension angle changes of each joint in the sagittal plane at a sampling frequency of 100 Hz. The angular velocity time series data of each joint are obtained by numerical differentiation of the angular time series. The numerical differentiation uses the central difference method; for the angular velocity at time t, the angle difference between time t+10 milliseconds and time t-10 milliseconds is taken and divided by 20 milliseconds. The phase difference of angular velocities between adjacent joint nodes is calculated. Specifically, the first joint pair consisting of the hip and knee joints and the second joint pair consisting of the knee and ankle joints are extracted. For the first joint pair, the hip joint angular velocity sequence and the knee joint angular velocity sequence are acquired synchronously during the observation period; the two sequences are of equal length and time-aligned.

[0072] The analytic signal transformation of the angular velocity phase difference is achieved through Hilbert transform. Applying Hilbert transform to the hip joint angular velocity sequence yields its orthogonal signal. The original angular velocity sequence is used as the real part, and the orthogonal signal as the imaginary part to construct a complex analytic signal. The amplitude information is extracted from the complex analytic signal to obtain the instantaneous phase sequence of the hip joint angular velocity. This sequence monotonically increases with time and cycles within the range of 0 to 360 degrees. The same transformation is performed on the knee joint angular velocity sequence to obtain the instantaneous phase sequence of the knee joint angular velocity. The difference between the two instantaneous phase sequences at each sampling time is calculated to obtain the phase difference sequence. The phase difference value ranges from -180 degrees to +180 degrees; a positive value indicates that the hip joint movement leads the knee joint movement, and a negative value indicates lag. The same processing procedure is performed on the knee and ankle joint angular velocity sequences of the second joint pair to obtain the second set of phase difference sequences.

[0073] Phase synchronization coefficients are calculated based on the temporal fluctuation characteristics of instantaneous phase sequences. The phase difference sequence of the first joint pair is segmented into time windows within the observation period, with a window length of 200 milliseconds and an overlap of 100 milliseconds between adjacent windows. Within each time window, each phase difference value in the phase difference sequence is converted into a unit vector on the complex plane. Specifically, for the phase difference value θ, a complex number is constructed with its real part being the cosine of the angle and its imaginary part being the sine of the angle. The average of all complex vectors within the window is calculated to obtain an average complex vector, and the magnitude of this average complex vector is the phase lock value for that window. The phase lock value is between zero and 1; a value close to 1 indicates that the phase difference remains stable within the window and the two joints are highly synchronized, while a value close to zero indicates that the phase difference fluctuates wildly and the coordination is poor. The phase lock values ​​of each time window are arranged in chronological order to form the phase synchronization coefficient sequence of the first joint pair. The same calculation process is performed on the second joint pair to obtain the phase synchronization coefficient sequence of the second joint pair.

[0074] A phase synchronization coefficient sequence was formed by continuously sampling the phase synchronization coefficients during the observation period. The sampling interval of this sequence was 100 milliseconds, consistent with the time window overlap strategy. The switching frequency between different coordination states in the phase synchronization coefficient sequence was statistically analyzed, and three coordination state thresholds were defined: a phase synchronization coefficient greater than 0.8 was a high coordination state, between 0.5 and 0.8 was a medium coordination state, and less than 0.5 was a low coordination state. The phase synchronization coefficient sequence was traversed, and each sampling point was assigned to a corresponding state according to its value. The state category of two adjacent sampling points was recorded. If the state changed, it was recorded as a state transition. Nine possible state transition types and their occurrence frequency were statistically analyzed during the observation period: high coordination to high coordination, high coordination to medium coordination, high coordination to low coordination, medium coordination to high coordination, medium coordination to medium coordination, medium coordination to low coordination, low coordination to high coordination, low coordination to medium coordination, and low coordination to low coordination. The frequencies of the nine transition types were filled into a three-row, three-column state transition frequency matrix, where the element in the i-th row and j-th column represents the number of times the state transitioned from the i-th state to the j-th state.

[0075] The state transition frequency matrix is ​​normalized by summing all elements and dividing each element by the sum to obtain the normalized transition frequency matrix. The normalized element values ​​represent the probability of the corresponding transition type occurring. The weight distribution of the dominant transition path representing the stable coordination mode is extracted from the normalized transition frequency matrix. The stable coordination mode corresponds to three transition paths: high coordination to high coordination, medium coordination to medium coordination, and low coordination to low coordination, i.e., state self-transition paths. The values ​​of the three elements on the diagonal of the matrix are extracted, representing the probability of maintaining stability in the three states, and these three values ​​are combined to form a weight distribution vector. This weight distribution vector reflects the stability distribution characteristics of gait at different coordination levels. Healthy individuals typically show a dominant weight in the high coordination state, while patients with balance disorders may show a higher weight in the low coordination state or a dispersed weight distribution across states.

[0076] The weight distribution and phase synchronization coefficient sequence are fused to generate an enhanced joint coordination feature vector. Three statistics—sequence mean, standard deviation, and maximum value—are calculated for the phase synchronization coefficient sequence of the first joint pair, forming a three-dimensional basic feature vector. This basic feature vector is then multiplied element-wise with the weight distribution vector: specifically, the first element of the basic feature vector is multiplied by the weight of the high-coordination state, the second element by the weight of the medium-coordination state, and the third element by the weight of the low-coordination state, resulting in a weighted three-dimensional vector. This weighted vector is concatenated with the weight distribution vector to form a six-dimensional enhanced feature vector. The same fusion process is performed on the second joint pair to obtain a second six-dimensional enhanced feature vector. The two six-dimensional vectors are merged to finally generate a twelve-dimensional enhanced joint coordination feature vector. This vector simultaneously encodes the statistical characteristics of phase synchronization and the state stability distribution, comprehensively characterizing the coordination patterns between lower limb joints.

[0077] This invention, through state transition frequency analysis and weight distribution fusion mechanism, can further enhance the characterization ability of joint coordination features on gait stability based on phase synchronization coefficient.

[0078] In one optional implementation, a multi-objective optimization function is constructed, the multi-objective optimization function including an optimization objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, and the step of solving the Pareto optimal solution set of the multi-objective optimization function includes:

[0079] A multi-objective optimization function is constructed, which includes a first sub-objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, a second sub-objective of minimizing the probability that the pressure center trajectory offset exceeds the safe range, and a third sub-objective of minimizing the fluctuation amplitude of the joint coordination feature.

[0080] During the continuous acquisition of the foot pressure time series data, the time-domain sequence of the offset of the pressure center trajectory relative to the preset reference area is calculated, and the peak fluctuation amplitude and steady-state residual of the offset are extracted as adaptation evaluation indicators through sliding window statistical analysis.

[0081] The weight coefficients among the first, second, and third sub-objectives are dynamically adjusted based on the deviation between the adaptation evaluation index and the preset adaptation threshold. When the adaptation evaluation index is lower than the preset adaptation threshold, the weight coefficient of the first sub-objective is increased; when the adaptation evaluation index is higher than the preset adaptation threshold, the weight coefficient of the first sub-objective is decreased. The Pareto optimal solution set of the multi-objective optimization function is then solved under the updated weight coefficients.

[0082] For example, a multi-objective optimization function is constructed, which includes three mutually constraining sub-objectives. The first sub-objective is to maximize the response amplitude of the defect feature vector under perturbation excitation. This is achieved by convolving each dimensional component of the defect feature vector with the perturbation excitation signal to obtain response time-series data. The peak-to-peak value of the response time-series data during the observation period is calculated as the response amplitude index. A larger index value indicates a higher sensitivity of the system to perturbations, which is beneficial for early detection of the degradation trend of balance capability. The second sub-objective is to minimize the probability that the pressure center trajectory offset exceeds the safe range. The safe range is defined as a rectangular area with the geometric center of the foot as the origin, ±30 mm in the front-back direction, and ±20 mm in the left-right direction. The probability of exceeding the safe range is obtained by dividing the number of sampling points where the pressure center trajectory point falls outside the safe range by the total number of sampling points during the observation period. A smaller probability indicates better attitude control stability. The third sub-objective is to minimize the fluctuation amplitude of joint coordination features. The difference between the maximum and minimum values ​​of each dimensional component of the joint coordination feature vector during the observation period is calculated, and the sum of all dimensional differences is obtained as the total fluctuation amplitude. A smaller amplitude indicates stronger consistency of motion patterns.

[0083] During the continuous acquisition of foot pressure time-series data, pressure values ​​from each sensing unit were collected at a frequency of 50 Hz using a plantar pressure sensor array. At each sampling time, the weighted sum of the products of all sensing unit pressure values ​​and their corresponding position coordinates was divided by the total pressure value to obtain the pressure center coordinates. The pressure center coordinate sequence was compared with a preset reference area to calculate the offset time-domain sequence. The preset reference area was the average pressure center position measured when the subject was statically standing. The offset was defined as the Euclidean distance between the current pressure center coordinates and the reference position coordinates. A sliding window statistical analysis was used to extract the peak fluctuation amplitude and steady-state residual of the offset as adaptation evaluation indicators. The sliding window length was set to 2 seconds, and the step interval was 0.5 seconds. Within each window, the peak fluctuation amplitude was calculated as the maximum value minus the minimum value of the offset sequence, and the steady-state residual was calculated as the standard deviation of the offset sequence. The comprehensive adaptation evaluation index for that window was obtained by multiplying the peak fluctuation amplitude by 0.6 and the steady-state residual by 0.4. This index reflects the subject's balance control quality during the current time period.

[0084] The preset adaptation threshold is determined based on the median of the subject's historical data; for healthy adults, this threshold is typically 15 mm. The deviation between the current adaptation evaluation index and the preset adaptation threshold is calculated: Deviation = Adaptation Evaluation Index - Preset Adaptation Threshold. When the deviation is less than zero, meaning the adaptation evaluation index is below the preset adaptation threshold, it indicates that the subject's balance ability still has a margin of error. In this case, the weight coefficient of the first sub-target is increased to enhance the intensity of the disturbance stimulus and uncover potential defect characteristics. The weight adjustment rule is: increase the weight coefficient of the first sub-target from the initial value of 0.4 to 0.4 + absolute deviation value × 0.01, while proportionally decreasing the weight coefficients of the second and third sub-targets, keeping the sum of the three weight coefficients constant at 1.0. When the deviation is greater than zero, meaning the adaptation evaluation index is above the preset adaptation threshold, it indicates that the subject's balance control is approaching its upper limit. In this case, the weight coefficient of the first sub-target is decreased to reduce the intensity of the disturbance and avoid the risk of falling. The weight adjustment rule is: decrease the weight coefficient of the first sub-target from the current value to the current value - deviation value × 0.01, with a lower limit of 0.2, and increase the weight coefficient of the second sub-target to strengthen safety constraints. The weight coefficient of the third sub-objective remained relatively stable throughout the adjustment process, with a value ranging from 0.2 to 0.3.

[0085] The Pareto optimal solution set of the multi-objective optimization function is obtained under the updated weight coefficients. A non-dominated sorting genetic algorithm is used to perform the solution process, with a population size of 100 individuals. Each individual is encoded as a set of perturbation excitation parameters, including three decision variables: perturbation amplitude, frequency, and duration. The perturbation amplitude ranges from 0 to 50 Newtons, the frequency ranges from 0.5 to 5 Hz, and the duration ranges from 0.5 to 3 seconds. For each individual in the population, its fitness function value under the three sub-objectives is calculated: fitness for the first sub-objective = response amplitude × weight coefficient of the first sub-objective; fitness for the second sub-objective = probability of exceeding the target × weight coefficient of the second sub-objective × (-1); fitness for the third sub-objective = fluctuation amplitude × weight coefficient of the third sub-objective × (-1). Non-dominated sorting is performed to divide the individuals into different levels of Pareto front layers. The first layer contains non-dominated solutions not dominated by any other individual. The offspring population was generated through tournament selection, simulated binary crossover, and polynomial mutation operations, with the crossover probability set to 0.9 and the mutation probability set to 0.1. After 200 generations of iterative evolution, the process converged, and all individuals in the first Pareto front were extracted as the Pareto optimal solution set, which contained 15 to 25 uniformly distributed non-dominated solutions.

[0086] This invention, through a dynamic weight adjustment mechanism and Pareto optimal solution set solving, can adaptively optimize the perturbation excitation strategy while ensuring safety, thereby improving the accuracy and personalization of balance function evaluation.

[0087] In an optional implementation, the step of dynamically adjusting the weight coefficients among the first sub-objective, the second sub-objective, and the third sub-objective based on the deviation between the adaptation evaluation index and a preset adaptation threshold, and solving for the Pareto optimal solution set of the multi-objective optimization function under the updated weight coefficients includes:

[0088] Calculate the relative deviation rate between the adaptation evaluation index and the preset adaptation threshold and perform normalization processing; set the initial weight coefficients of the first sub-objective, the second sub-objective and the third sub-objective, wherein the initial weight coefficient of the first sub-objective is higher than that of the other sub-objectives;

[0089] The weight adjustment factor is calculated using a piecewise nonlinear function based on the normalized relative deviation rate. The slope of the adjustment factor's response to the deviation rate within each interval of the piecewise nonlinear function satisfies an increasing relationship. The weight coefficients of each sub-objective are updated based on the weight adjustment factor while maintaining the weight normalization constraint.

[0090] A multi-objective evolutionary optimization algorithm is used to solve the multi-objective optimization function after updating the weights. A search direction guidance vector is constructed based on the weight distribution characteristics of the defect feature vector and used as the directional perturbation term for the population mutation operation. The proportion of individuals in the current population whose pressure center trajectory offset exceeds the safe range is statistically analyzed, and the relationship with the preset proportion threshold is used to adjust the population size and mutation probability in reverse. Individuals in the population are screened and retained according to the target space crowding distance.

[0091] For example, the relative deviation rate between the adaptation evaluation index and the preset adaptation threshold is calculated and normalized. The adaptation evaluation index is obtained through sliding window statistics, and the preset adaptation threshold is determined based on the median of the subjects' historical data. For healthy adults, this threshold is typically 15 mm. The calculation steps for the relative deviation rate are as follows: First, calculate the difference between the adaptation evaluation index and the preset adaptation threshold; then divide this difference by the preset adaptation threshold to obtain the relative deviation rate. This ratio is a dimensionless parameter; a positive value indicates that the current balance control quality is lower than the expected level, and a negative value indicates that it is higher than the expected level. The normalization process uses a hyperbolic tangent mapping function. The specific steps are: first, multiply the relative deviation rate by a scaling factor of 2.5 to obtain an intermediate variable; then use this intermediate variable as the input independent variable of the hyperbolic tangent function for calculation, and the output value range is compressed to the interval from -1 to +1. The selection of the scaling factor of 2.5 allows the nonlinear segment of the mapping curve to be fully expanded within the range of -0.4 to +0.4 for the relative deviation rate. The normalized relative deviation rate is retained to three significant figures, and the error tolerance is controlled within 0.001.

[0092] Initial weight coefficients are set for the first, second, and third sub-objectives, with the first sub-objective having a higher initial weight coefficient than the others. The first sub-objective corresponds to maximizing the response amplitude of the defect feature vector under perturbation excitation, and its initial weight coefficient is set to 0.4. The second sub-objective corresponds to minimizing the probability of the pressure center trajectory deviation exceeding the safe range, and its initial weight coefficient is set to 0.35. The third sub-objective corresponds to minimizing the fluctuation amplitude of joint coordination characteristics, and its initial weight coefficient is set to 0.25. The sum of the three initial weight coefficients equals 1.0, satisfying the normalization constraint. The initial weight configuration reflects the strategy of prioritizing the detection of defect features in the early stages of the balance function assessment. The first sub-objective has a dominant weight, while the second and third sub-objectives provide constraint balance from the perspectives of safety and stability, respectively.

[0093] The weight adjustment factor is calculated using a piecewise nonlinear function based on the normalized relative deviation rate. The piecewise nonlinear function divides the domain of the normalized relative deviation rate from -1 to +1 into five intervals, with the boundary points of each interval being -1, -0.5, -0.2, +0.2, +0.5, and +1, respectively. Within the interval from -1 to -0.5, the calculation steps for the weight adjustment factor are as follows: first, calculate the sum of the normalized relative deviation rate and 1; then, multiply this sum by 0.15; finally, add the product to 0.80 to obtain the weight adjustment factor. The slope of the adjustment factor in this interval is 0.15. Within the interval from -0.5 to -0.2, the calculation steps for the weight adjustment factor are as follows: first, calculate the sum of the normalized relative deviation rate and 0.5; then, multiply this sum by 0.20; finally, add the product to 0.95 to obtain the weight adjustment factor. The slope of the response increases to 0.20. Within the range of -0.2 to +0.2, the weight adjustment factor is calculated as follows: First, multiply the normalized relative deviation rate by 0.30; then add the product to 1.00 to obtain the weight adjustment factor, further increasing the response slope to 0.30. Within the range of +0.2 to +0.5, the weight adjustment factor is calculated as follows: First, calculate the difference between the normalized relative deviation rate and 0.2; then multiply this difference by 0.45; finally, add the product to 1.06 to obtain the weight adjustment factor, increasing the response slope to 0.45. Within the range of +0.5 to +1, the weight adjustment factor is calculated as follows: First, calculate the difference between the normalized relative deviation rate and 0.5; then multiply this difference by 0.60; finally, add the product to 1.20 to obtain the weight adjustment factor, increasing the response slope to 0.60. The slopes of the adjustment factors in response to the deviation rate within each interval are 0.15, 0.20, 0.30, 0.45, and 0.60, respectively, showing a progressively increasing characteristic. This realizes an adaptive mechanism in which the weight adjustment amplitude increases as the balance capability deviates further from the threshold.

[0094] The weight coefficients of each sub-objective are updated based on a weight adjustment factor while maintaining weight normalization constraints. The updated weight coefficient of the first sub-objective is calculated as the product of the initial weight coefficient 0.4 and the weight adjustment factor. When the weight adjustment factor is greater than 1.0, the weight of the first sub-objective increases; when the weight adjustment factor is less than 1.0, the weight of the first sub-objective decreases. The updated weight coefficient of the second sub-objective is calculated as follows: first, calculate the difference between 2 and the weight adjustment factor; then multiply this difference by the initial weight coefficient 0.35 to obtain the updated weight coefficient of the second sub-objective. This calculation rule ensures that the weights of the second and first sub-objectives show an inverse adjustment trend. The updated weight coefficient of the third sub-objective is obtained by reverse deduction through normalization constraints. The calculation steps are: first, subtract the updated weight coefficient of the first sub-objective from 1.0; then subtract the updated weight coefficient of the second sub-objective from the resulting difference to obtain the updated weight coefficient of the third sub-objective, ensuring that the sum of the three updated weight coefficients is always equal to 1.0. The updated weight coefficients retain four decimal places, and the boundary conditions are constrained. The weight coefficient of the first sub-objective is limited to the range of 0.2 to 0.7, the weight coefficient of the second sub-objective is limited to the range of 0.15 to 0.50, and the weight coefficient of the third sub-objective is limited to the range of 0.15 to 0.35. This avoids the excessive expansion or contraction of the weight of a certain sub-objective, which would cause the multi-objective optimization to degenerate into a single-objective optimization.

[0095] A multi-objective evolutionary optimization algorithm is employed to solve the multi-objective optimization function after weight updates. The algorithm is implemented based on a non-dominated sorting genetic algorithm framework, with an initial population size of 100 individuals. Each individual is encoded as a set of perturbation-incentive decision variables, including three dimensions: perturbation amplitude, frequency, and duration. The perturbation amplitude is encoded as real numbers, ranging from 0 to 50 Newtons; the frequency ranges from 0.5 to 5 Hz; and the duration ranges from 0.5 to 3 seconds. A search direction guidance vector is constructed based on the weight distribution characteristics of the defect feature vector and used as the directional perturbation term for population mutation operations. The defect feature vector contains 12 dimensional components. The absolute values ​​of each component are sorted in descending order, and the dimensional indices corresponding to the top four largest components are extracted. Guiding vectors are constructed in the decision variable space based on these dimensional indices. The calculation methods for each dimension component of the guiding vector are as follows: the perturbation amplitude dimension component equals the product of the absolute value of the maximum component and 8 Newtons per unit; the frequency dimension component equals the product of the absolute value of the maximum component and 0.4 Hertz per unit; and the duration dimension component equals the product of the absolute value of the maximum component and 0.18 seconds per unit. In the population mutation operation, the selected individual decision variables are updated. The update steps are: first, calculate the product of the guiding vector and the random perturbation intensity coefficient; then, add this product to the original decision variable to obtain the new decision variable. The perturbation intensity coefficient follows a truncated normal distribution with a mean of 0.25 and a standard deviation of 0.08, and its value range is limited to the interval between 0.1 and 0.4.

[0096] The proportion of individuals in the current population whose pressure center trajectory offset exceeds the safe range is statistically analyzed, and the relationship between this proportion and a preset threshold is used to inversely adjust the population size and mutation probability. For each individual in the population, the maximum pressure center trajectory offset under the perturbation excitation scheme is calculated. This maximum offset is compared with the safe range boundary value, defined as ±30 mm in the forward / backward direction and ±20 mm in the left / right direction. The proportion of individuals whose maximum offset exceeds the safe range is calculated by dividing the number of individuals by the total number of individuals in the population. The preset threshold for this proportion is set to 0.20. When the proportion of individuals exceeding the safe range is less than the preset threshold, it indicates that the overall safety of the current population is relatively good. In this case, the population size and mutation probability are increased to expand the search space. The steps for adjusting the population size are as follows: First, calculate the difference between 0.20 and the proportion of individuals exceeding the safe range; then multiply this difference by 50 to obtain the increment coefficient; then add the increment coefficient to 1.0 to obtain the scaling factor; finally, multiply the current size by the scaling factor to obtain the adjusted population size. The steps for adjusting the mutation probability are as follows: First, calculate the difference between 0.20 and the excess individual proportion; then multiply this difference by 0.25 to obtain the increment coefficient; next, add the increment coefficient to 1.0 to obtain the scaling factor; finally, multiply the current mutation probability 0.1 by the scaling factor to obtain the adjusted mutation probability. When the excess individual proportion is greater than or equal to a preset proportion threshold, it indicates that the current population has a high risk. At this time, the population size and mutation probability are reduced to converge to a safe region. The steps for adjusting the population size are as follows: First, calculate the difference between the excess individual proportion and 0.20; then multiply this difference by 40 to obtain the reduction coefficient; next, calculate the difference between 1.0 and the reduction coefficient to obtain the scaling factor; finally, multiply the current size by the scaling factor to obtain the adjusted population size. The steps for adjusting the mutation probability are as follows: First, calculate the difference between the excess individual proportion and 0.20; then multiply this difference by 0.20 to obtain the reduction coefficient; next, calculate the difference between 1.0 and the reduction coefficient to obtain the scaling factor; finally, multiply the current mutation probability 0.1 by the scaling factor to obtain the adjusted mutation probability. The range of values ​​for population size after adjustment is limited to 80 to 140, and the range of values ​​for mutation probability after adjustment is limited to 0.06 to 0.18.

[0097] Individuals in the population are selected and retained based on their crowding distance in the target space. The crowding distance of each individual in the population is calculated in the three-dimensional target space. The crowding distance is defined as the normalized sum of the differences between the individual's target value and those of its neighboring individuals in each target dimension. For the first sub-target dimension, the calculation steps for the normalized distance component are as follows: First, obtain the target value of the individual's predecessor and successor in that dimension; then, calculate the difference between these two target values ​​as the numerator; next, obtain the maximum and minimum target values ​​of all individuals in that dimension and calculate their difference as the denominator; finally, divide the numerator by the denominator to obtain the normalized distance component for the first sub-target dimension. The same calculation process is performed for the second and third sub-target dimensions, summing the normalized distance components of the three dimensions to obtain the individual's crowding distance. Boundary individuals, i.e., those with the maximum or minimum value in a certain target dimension, have their crowding distance set to infinity to ensure priority retention. The population individuals are sorted in descending order of their crowding distance, and the top N individuals with the largest crowding distances are retained for the next generation of the population, where N is the adjusted population size.

[0098] This invention achieves adaptive safety optimization of the perturbation excitation strategy through piecewise nonlinear weight adjustment and defect feature-guided evolutionary optimization, thereby improving the accuracy and personalization of balance function evaluation.

[0099] In one optional implementation, the steps of extracting feature parameters characterizing the user's sensitivity to balance defects from the Pareto optimal solution set, performing a graded and quantitative assessment of the user's dynamic balance ability based on the feature parameters, and generating a balance training requirement assessment report containing defect type identification results, sensitivity level classification, and corresponding training suggestion parameters include:

[0100] The ratio of the response magnitude of the defect feature vector of each solution in the Pareto optimal solution set to the degree of constraint violation is calculated as the defect sensitivity index, and the sensitivity level classification is determined based on the distribution of the defect sensitivity index.

[0101] From the Pareto optimal solution set, candidate safe solutions with constraint violation degree below the preset safety threshold are selected. The solution with the largest response amplitude of the defect feature vector in the candidate safe solution set is selected as the optimal evaluation solution. The defect feature vector weight distribution pattern of the optimal evaluation solution is matched with the preset defect pattern library for similarity, and the defect type identifier with the highest similarity is output as the defect type identification result.

[0102] Based on the defect type identification result, the corresponding benchmark training parameters are queried. The benchmark training parameters include the training intensity range and training mode category. The benchmark training parameters are adaptively adjusted according to the sensitivity level classification, and the adjusted parameters are used as training suggestion parameters.

[0103] The defect type identification results, the sensitivity level classification, and the training suggestion parameters are organized to generate a balanced training requirement assessment report.

[0104] For example, the defect sensitivity index of each solution in the Pareto optimal solution set is calculated, and the sensitivity level classification is determined. The Pareto optimal solution set contains 15 to 25 non-dominated solutions, each corresponding to a set of perturbation excitation parameters and its performance under three sub-objectives. The response amplitude of the defect feature vector is obtained by convolving each dimension component of the defect feature vector with the perturbation excitation signal corresponding to the solution to obtain the response time series data. The peak-to-peak value of the response time series data within the observation period is extracted as the response amplitude value. The larger the value, the higher the sensitivity to perturbation. The degree of constraint violation is calculated by comprehensively evaluating the deviations of the second and third sub-objectives. The specific steps are as follows: First, obtain the probability value of the pressure center trajectory offset exceeding the safe range and the fluctuation amplitude value of the joint coordination characteristics; then, divide the excess probability value by a preset safe probability upper limit of 0.25 to obtain the first normalized violation degree, and divide the fluctuation amplitude value by a preset fluctuation upper limit of 8.0 to obtain the second normalized violation degree; next, multiply the first normalized violation degree by a weight of 0.6, and multiply the second normalized violation degree by a weight of 0.4; finally, add the two weighted values ​​to obtain the degree of constraint violation for the solution. The calculation steps for the defect sensitivity index are as follows: divide the defect feature vector response amplitude by the sum of the constraint violation degrees and add 1 (this addition avoids a zero denominator); perform a logarithmic transformation on the quotient to compress the numerical range, thus obtaining the defect sensitivity index. Statistical analysis was performed on the defect sensitivity indices of all solutions in the Pareto optimal solution set. The quartiles Q1, Q2, and Q3 of the index sequence were calculated, dividing the index space into four levels: index less than Q1 indicates low sensitivity; index between Q1 and Q2 indicates low to medium sensitivity; index between Q2 and Q3 indicates medium to high sensitivity; and index greater than Q3 indicates high sensitivity. The sensitivity level classification directly reflects the response characteristics of the user's equilibrium system to specific disturbance patterns. A high sensitivity level indicates a significant risk of functional degradation in the corresponding defect dimension.

[0105] From the Pareto optimal solution set, candidate safe solutions with constraint violation levels below a preset safety threshold are selected. The preset safety threshold is set to 0.35 based on clinical safety standards. Each solution in the Pareto optimal solution set is iterated, and its constraint violation level is compared to 0.35. Solutions with constraint violation levels strictly less than 0.35 are retained to form the candidate safe solution set. The candidate safe solution set typically contains 8 to 15 solutions, which possess sufficient defect detection capabilities while ensuring training safety. The solution with the largest defect feature vector response amplitude is selected as the optimal evaluation solution from the candidate safe solution set. The selection method involves iterating through the candidate safe solution set, calculating the response amplitude of each solution, recording the solution index corresponding to the largest response amplitude, and extracting the complete solution parameters corresponding to that index. The optimal evaluation solution simultaneously satisfies the requirements of safety constraints and sensitivity maximization; its perturbation excitation parameters can most fully stimulate the balanced defect features without threatening user safety. The defect feature vector weight distribution pattern of the optimal evaluation solution is extracted. The defect feature vector is a six-dimensional vector, which is normalized so that the sum of the absolute values ​​of the weights of each dimension equals 1. The normalized vector retains the sign information of each dimension. The pre-defined defect pattern library stores standard weight distribution patterns for 12 typical balance defect types, including vestibular dysfunction, proprioceptive loss, cerebellar coordination disorder, muscle weakness, joint stiffness, and gait asymmetry. Each pattern corresponds to a six-dimensional standard weight vector. Similarity matching employs cosine similarity calculation, performing a dot product operation between the normalized weight vector of the optimal evaluation solution and each standard pattern vector in the pre-defined defect pattern library. The dot product result is the similarity score for that pattern, ranging from -1 to +1, with values ​​closer to 1 indicating a higher matching degree. The similarity scores of the 12 standard patterns are iterated, and the pattern index corresponding to the highest similarity score is extracted. The defect type identifier for this pattern is output as the defect type identification result. The defect type identifier includes the type name, feature description, and clinical significance, providing a pathological basis for subsequent training program customization.

[0106] The corresponding baseline training parameters are retrieved from the training parameter database based on the defect type identification results. The training parameter database maintains a mapping relationship between 12 defect types and training schemes, with each defect type corresponding to a set of baseline training parameters. The baseline training parameters include two core elements: training intensity range and training mode category. The training intensity range is defined as the lower and upper limits of the perturbation amplitude, in Newtons, with a typical range of 10 to 40 Newtons. The training mode category is defined as a combination of the time-frequency characteristics of the perturbation signal, including five basic categories: constant perturbation mode, intermittent pulse mode, sinusoidal modulation mode, random white noise mode, and multi-frequency superposition mode. The baseline training parameters are adaptively adjusted based on sensitivity level classification, as follows: For low-sensitivity users, the lower limit of the training intensity range remains unchanged, while the upper limit is increased by the product of the initial upper limit and 0.2. The training mode is selected as a more aggressive random white noise mode or multi-frequency superposition mode to increase stimulation intensity. For low-to-medium sensitivity users, the training intensity range remains unchanged, and the training mode is selected as an intermittent pulse mode or sinusoidal modulation mode to achieve moderate stimulation. For medium-to-high sensitivity users, the upper limit of the training intensity range is decreased by the product of the initial upper limit and 0.15, while the lower limit remains unchanged. The training mode is selected as a constant perturbation mode or intermittent pulse mode to reduce training risk. For high-sensitivity users, the upper limit of the training intensity range is decreased by the product of the initial upper limit and 0.3, while the lower limit is decreased by the product of the initial lower limit and 0.2. The training mode is forcibly selected as a constant perturbation mode with a progressive intensity increase strategy. The adjusted upper and lower limits of the training intensity range are all rounded down to the nearest integer. If the calculated result is not an integer, it is rounded down. The training recommendations also include training frequency recommendations and single training session duration recommendations. The training frequency is set to 3 to 7 times per week based on the sensitivity level, with a lower frequency corresponding to a higher sensitivity level and a higher frequency corresponding to a lower sensitivity level. The single training session duration is set to 15 to 45 minutes based on the degree of constraint violation, with shorter durations for higher constraint violation levels to control fatigue accumulation.

[0107] The evaluation report uses a structured text format and includes four required fields: User Identifier (recording the subject's unique ID and evaluation timestamp); Deficiency Type (recording a complete text description of the deficiency type identifier and its similarity score); Sensitivity Level (recording the four-level classification results and the corresponding deficiency sensitivity index); and Training Recommendation (recording the adjusted training intensity range, training mode category, training frequency recommendations, single training duration recommendations, and special precautions). The report generation process utilizes a template engine, filling the content of each field into a predefined report template. The template includes a title area, a data area, and a recommendation area. The data area presents key numerical parameters in tabular form, while the recommendation area describes the key points of the training program and expected results in paragraph form. Report output supports two formats: structured JSON for inter-system data exchange and a human-readable PDF for clinicians and patient communication. Reports are stored in an evaluation record database with a version management mechanism. A new version report is generated for each evaluation, while historical reports are retained to support longitudinal comparative analysis.

[0108] This invention achieves precise matching of personalized training parameters through defect sensitivity index quantification and safety constraint screening mechanism, thereby improving the pertinence and safety of balance rehabilitation training.

[0109] A second aspect of the present invention provides an electronic device, comprising:

[0110] processor;

[0111] Memory used to store processor-executable instructions;

[0112] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0113] A third aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0114] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0115] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for dynamic assessment of balance ability and fall prevention training based on time-series data, characterized in that, include: Simultaneously collect foot pressure time-series data and joint motion data of users on the mobile training platform; The pressure center trajectory in the foot pressure time series data is divided into continuous observation periods according to the time dimension. For each observation period, a coupling correlation matrix is ​​constructed between the spatial morphological features of the pressure center trajectory and the joint coordination features of the joint motion data. The coupling correlation matrix is ​​then decomposed by eigenvalue to extract defect feature vectors. Construct a multi-objective optimization function, which includes an optimization objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, and solve the Pareto optimal solution set of the multi-objective optimization function; Feature parameters characterizing the user's sensitivity to balance defects are extracted from the Pareto optimal solution set. Based on the feature parameters, the user's dynamic balance ability is evaluated in a graded and quantitative manner, and a balance training requirement evaluation report is generated, which includes defect type identification results, sensitivity level classification and corresponding training suggestion parameters. The step of constructing a coupling correlation matrix between the spatial morphological features of the pressure center trajectory and the joint coordination features of the joint motion data for each observation period, and extracting defect feature vectors by eigenvalue decomposition of the coupling correlation matrix includes: For the trajectory of the pressure center within each observation period, the rate of change of the direction of movement and the rate of change of curvature are calculated. The rate of change of the direction of movement and the rate of change of curvature are then decomposed into time-frequency components, and the frequency domain energy distribution vector is extracted as a spatial morphological feature. For the joint motion data corresponding to the observation period, the angular velocity phase difference between adjacent joints is extracted, the angular velocity phase difference is continuously sampled, and the phase synchronization coefficient sequence is calculated as a joint coordination feature; A coupling correlation matrix is ​​constructed using the spatial morphological features as row vectors and the joint coordination features as column vectors; The coupling correlation matrix is ​​decomposed into eigenvalues, and the eigenvector corresponding to the largest eigenvalue is extracted as the defect feature vector. The weight distribution of each element in the defect feature vector reflects the degree of influence of the spatial morphological features on the joint coordination features.

2. The method according to claim 1, characterized in that, The steps of performing time-frequency decomposition on the rate of change of the direction of movement and the rate of change of curvature to extract the frequency domain energy distribution vector include: Wavelet packet decomposition is performed on the rate of change of the direction of movement and the rate of change of curvature respectively, and the frequency band is divided into three sub-bands: low-frequency steady-state component, mid-frequency transition component and high-frequency abrupt component. The energy proportion of each sub-band is calculated to form an initial frequency domain energy distribution vector. The trajectory smoothness is evaluated by calculating the autocorrelation coefficient of the low-frequency steady-state component. The autocorrelation coefficient is then weighted and fused with the initial frequency domain energy distribution vector. When the autocorrelation coefficient is lower than a preset smoothing threshold, the energy weight of the high-frequency abrupt component is increased. When the autocorrelation coefficient is higher than the preset smoothing threshold, the energy weight of the low-frequency steady-state component is increased, thereby generating a corrected frequency domain energy distribution vector.

3. The method according to claim 1, characterized in that, For the joint motion data corresponding to the observation period, the steps of extracting the angular velocity phase difference between adjacent joints, continuously sampling the angular velocity phase difference, and calculating the phase synchronization coefficient sequence as a joint coordination feature include: Extract the angular velocity time-series data of the joint nodes in the lower limb kinetic chain from the joint motion data, and calculate the angular velocity phase difference between adjacent joint nodes; The instantaneous phase sequence is extracted by performing analytical signal transformation on the angular velocity phase difference, and the phase synchronization coefficient is calculated based on the time-domain fluctuation characteristics of the instantaneous phase sequence. During the observation period, the phase synchronization coefficients are continuously sampled to form a phase synchronization coefficient sequence; The switching frequency between different coordination states in the phase synchronization coefficient sequence is statistically analyzed, a state transition frequency matrix is ​​constructed and normalized, the weight distribution of the dominant transition path representing the stable coordination mode in the normalized transition frequency matrix is ​​extracted, and the weight distribution is fused with the phase synchronization coefficient sequence to generate an enhanced joint coordination feature vector.

4. The method according to claim 1, characterized in that, Constructing a multi-objective optimization function, wherein the multi-objective optimization function includes an optimization objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, and the steps of solving the Pareto optimal solution set of the multi-objective optimization function include: A multi-objective optimization function is constructed, which includes a first sub-objective of maximizing the response amplitude of the defect feature vector under perturbation excitation, a second sub-objective of minimizing the probability that the pressure center trajectory offset exceeds the safe range, and a third sub-objective of minimizing the fluctuation amplitude of the joint coordination feature. During the continuous acquisition of the foot pressure time series data, the time-domain sequence of the offset of the pressure center trajectory relative to the preset reference area is calculated, and the peak fluctuation amplitude and steady-state residual of the offset are extracted as adaptation evaluation indicators through sliding window statistical analysis. The weight coefficients among the first sub-target, the second sub-target, and the third sub-target are dynamically adjusted based on the deviation between the adaptation evaluation index and the preset adaptation threshold. When the adaptation evaluation index is lower than the preset adaptation threshold, the weight coefficient of the first sub-target is increased; when the adaptation evaluation index is higher than the preset adaptation threshold, the weight coefficient of the first sub-target is decreased. The Pareto optimal solution set of the multi-objective optimization function is obtained under the updated weighting coefficients.

5. The method according to claim 4, characterized in that, The steps of dynamically adjusting the weight coefficients among the first, second, and third sub-objectives based on the deviation between the adaptation evaluation index and the preset adaptation threshold, and solving for the Pareto optimal solution set of the multi-objective optimization function under the updated weight coefficients include: Calculate the relative deviation rate between the adaptation evaluation index and the preset adaptation threshold and perform normalization processing; Set initial weight coefficients for the first sub-target, the second sub-target, and the third sub-target, wherein the initial weight coefficient of the first sub-target is higher than that of the other sub-targets; The weight adjustment factor is calculated based on the normalized relative deviation rate using a piecewise nonlinear function. The slope of the adjustment factor's response to the deviation rate within each interval of the piecewise nonlinear function satisfies an increasing relationship. The weight coefficients of each sub-objective are updated based on the weight adjustment factor while maintaining the weight normalization constraint. A multi-objective evolutionary optimization algorithm is used to solve the multi-objective optimization function after updating the weights. A search direction guidance vector is constructed based on the weight distribution characteristics of the defect feature vector and used as the directional perturbation term for the population mutation operation. The population size and mutation probability are adjusted inversely by statistically analyzing the proportion of individuals in the current population whose stress center trajectory deviation exceeds the safe range and comparing it with a preset proportion threshold. Individuals in the population are selected and retained based on the target spatial crowding distance.

6. The method according to claim 1, characterized in that, The steps of extracting feature parameters characterizing the user's sensitivity to balance defects from the Pareto optimal solution set, performing a graded and quantitative assessment of the user's dynamic balance ability based on the feature parameters, and generating a balance training requirement assessment report that includes defect type identification results, sensitivity level classification, and corresponding training suggestion parameters include: The ratio of the response magnitude of the defect feature vector of each solution in the Pareto optimal solution set to the degree of constraint violation is calculated as the defect sensitivity index, and the sensitivity level classification is determined based on the distribution of the defect sensitivity index. From the Pareto optimal solution set, candidate safe solutions with constraint violation degree below the preset safety threshold are selected. The solution with the largest response amplitude of the defect feature vector in the candidate safe solution set is selected as the optimal evaluation solution. The defect feature vector weight distribution pattern of the optimal evaluation solution is matched with the preset defect pattern library for similarity, and the defect type identifier with the highest similarity is output as the defect type identification result. Based on the defect type identification result, the corresponding benchmark training parameters are queried. The benchmark training parameters include the training intensity range and training mode category. The benchmark training parameters are adaptively adjusted according to the sensitivity level classification, and the adjusted parameters are used as training suggestion parameters. The defect type identification results, the sensitivity level classification, and the training suggestion parameters are organized to generate a balanced training requirement assessment report.

7. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

8. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images