A dynamic gait parameter detection, evaluation and analysis method for postoperative rehabilitation of lower limbs

Abstract

This application relates to the field of postoperative rehabilitation assessment technology, and discloses a method for dynamic gait parameter detection, evaluation, and analysis in lower limb postoperative rehabilitation. This method simultaneously acquires acceleration, angular velocity, and plantar pressure sequences; uses plantar pressure to define the double support phase and generates a support phase impact gradient sequence; extracts full-plantar impact events for zero-velocity updates, inversely compensates for inertial sensor integral closure errors, and calculates the interlimb energy transfer ratio; constructs a dynamic weighted envelope using the support phase impact gradient, and performs dynamic weighted cross-correlation calculations on angular velocity within an extended data extraction buffer window to extract intralimb phase coupling delay; constructs a two-dimensional phase plane with energy transfer ratio and phase coupling delay as coordinate axes, analyzes the spatial geometric relationship between the current state point and the historical baseline, and calculates the compensatory attenuation index. This method eliminates integral drift, extracts the transient characteristics of the weight-bearing phase, and achieves an objective quantitative assessment of the patient's lower limb rehabilitation evolution trajectory.

Description

Technical Field

[0001] This invention relates to the field of postoperative rehabilitation assessment technology, specifically a method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation. Background Technology

[0002] Long-term gait monitoring is necessary for lower limb postoperative rehabilitation to assess the patient's neuromuscular recovery. Currently, wearable inertial sensors are commonly used to continuously collect patients' kinematic data. However, inertial sensors suffer from low-frequency cumulative drift errors during the time integration of acceleration. Because existing methods lack a reliable boundary to continuously correct for this cumulative error using the physical contact state between the foot and the ground, the calculated total mechanical energy of the lower limb and interlimb energy transfer parameters are distorted.

[0003] Based on this, existing technologies, when analyzing the phase coordination relationship between lower limb joints, typically assign equal weights to the temporal data of the entire gait cycle for cross-correlation calculations. This approach does not distinguish the force differences between different phase states. The transient motion changes caused by pain or muscle weakness at the moment of weight-bearing loading are masked by the data of the conventional swing phase, making it difficult to separate the true intralimb phase coupling delay characteristics of the affected lower limb.

[0004] Furthermore, existing gait assessment methods often examine single dynamic or kinematic parameters in isolation, lacking a coordinate reference system that integrates energy transfer with movement timing. Because it is impossible to make a unified spatial geometric comparison between the current test state and historical baselines and healthy states, the system struggles to objectively quantify the direction and extent of gait evolution during a patient's rehabilitation process. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a dynamic gait parameter detection, evaluation, and analysis method for lower limb postoperative rehabilitation. This method solves the problems of integral drift error in long-term monitoring by inertial sensors, difficulty in extracting phase delay features at the moment of weight-bearing in weight-bearing analysis of the overall gait cycle, and difficulty in objectively quantifying the evolution of rehabilitation trajectory due to isolated analysis of multidimensional gait parameters.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation, comprising the following steps: Simultaneously acquire the three-dimensional raw acceleration sequence, angular velocity sequence, and plantar pressure time sequence of the test subject's gait cycle; The original three-dimensional acceleration sequence is transformed into a spatial coordinate system and subjected to gravity stripping to generate a linear acceleration sequence. At the same time, the plantar pressure time series is used to calibrate plantar kinematic events to define the double support phase interval, and the support phase impact gradient sequence is generated based on the difference operation. The core integration window is defined based on the waveform characteristics of the supporting phase impact gradient sequence, and the data is extended and the buffer window is extracted based on the set maximum phase delay constant. Drift correction is performed on the linear acceleration sequence based on plantar kinematic events to solve the interlimb energy transfer ratio within the double support phase interval; a dynamic weight envelope is constructed using the support phase impact gradient sequence within the core integration window; and dynamic weighted cross-correlation is performed on the angular velocity sequence using the dynamic weight envelope within the data extraction buffer window to extract the intralimb phase coupling delay. A two-dimensional phase plane is constructed with the interlimb energy transfer ratio and intralimb phase coupling delay as coordinate axes. The spatial evolution geometric relationship of the current test state point relative to the centroid of the historical baseline is analyzed in the two-dimensional phase plane, and the compensatory decay index is calculated and output.

[0007] Furthermore, before synchronously acquiring sequence data, the inertial measurement unit is fixed at the outer center of the thigh, the outer center of the calf, and the instep or shoe upper of both lower limbs of the tester, respectively, and the plantar pressure sensor is placed inside both shoes of the tester; a synchronization interrupt is triggered using a global common hardware clock source at a preset discrete time step to perform sampling frequency matching and time axis alignment between the inertial measurement unit and the plantar pressure sensor, and the aligned multimodal time series data is pushed onto the stack and written into the circular buffer.

[0008] Furthermore, the process of calibrating plantar kinematic events to define the dual support phase interval includes: dividing the sole into the heel region, arch region, forefoot region, and toe region, and extracting the pressure time series of the corresponding regions and the total plantar pressure time series; when the pressure in the heel region crosses a preset noise threshold from bottom to top, it is defined as the heel contact event moment; when the total plantar pressure time series crosses the full contact pressure threshold from bottom to top and the plantar pressure in each region meets the preset plantar pressure distribution stability condition, it is defined as the start moment of the full plantar contact event; when the pressure in the toe region drops to the noise threshold and the total plantar pressure time series drops synchronously, it is defined as the toe lift-off event moment; using the heel contact event moment of one lower limb as the starting boundary and the toe lift-off event moment of the contralateral lower limb as the ending boundary, the dual support phase interval is formed.

[0009] Furthermore, the process of generating the support phase impact gradient sequence based on differential operation includes: extracting the plantar pressure time series of the affected side, calculating the transient loading rate between adjacent sampling points using the central difference algorithm; performing boundary condition processing at the beginning and end of the sequence using forward and backward differences respectively; and smoothing the differentiated sequence using the moving average algorithm to generate the support phase impact gradient sequence.

[0010] Furthermore, the process of defining the core integration window and extending the data extraction buffer window includes: traversing the support phase impact gradient sequence backward from the moment of the heel strike event on the affected side, defining the discrete time point where the gradient value crosses the gradient noise threshold from bottom to top as the starting point of the physical window; searching backward from the starting point of the physical window, defining the discrete time point where the waveform first changes from positive to negative after the peak as the ending point of the physical window; the core integration window is formed by the starting point and the ending point of the physical window; and extending the time axis of the core integration window to both sides by the maximum phase delay constant to form the data extraction buffer window.

[0011] Furthermore, the process of performing drift correction based on plantar kinematic events to calculate the interlimb energy transfer ratio includes: during the duration of a full-plantar landing event, performing zero-velocity updates on the integral velocity state of the inertial measurement unit at the dorsum of the foot or the upper of the shoe to obtain the foot velocity closure error and vertical displacement closure error; using a linear time decay function to inversely distribute the closure error, and continuously compensating for the integral velocity and vertical displacement of the thigh and lower leg segments of the ipsilateral lower limb; combining the velocity modulus and vertical displacement of each segment after compensation to calculate the mechanical energy of each lower limb segment, and summing the mechanical energy of the thigh segment, lower leg segment, and foot segment to obtain the total mechanical energy of the unilateral lower limb; within the double support period, performing time integration on the absolute value of the rate of change of the total mechanical energy of the affected lower limb and the healthy lower limb respectively, and calculating the ratio to calculate and output the interlimb energy transfer ratio.

[0012] Furthermore, the process of constructing the dynamic weight envelope includes: extracting the support phase impact gradient sequence within the core integration window and performing nonnegation processing; obtaining the local maximum absolute value of the gradient within the core integration window; and using the local maximum absolute value to perform normalization operation on each nonnegated gradient scalar within the core integration window to generate a dynamic weight envelope with values ​​ranging from zero to one.

[0013] Furthermore, the process of extracting the intralimb phase coupling delay includes: extracting the coronal axis rotation component orthogonal to the sagittal plane of the human body from the angular velocity sequence of the affected side to form the angular velocity scalar sequence of the thigh side and the angular velocity scalar sequence of the lower leg side; calculating the product integral of the angular velocity scalar sequence of the thigh side, the time-shifted angular velocity scalar sequence of the lower leg side, and the dynamic weight envelope at the corresponding time with the time translation step within the traversal range of the positive and negative maximum phase delay constants; obtaining the time translation term corresponding to the maximum value of the weighted cross-correlation coefficient as the intralimb phase coupling delay.

[0014] Furthermore, the process of constructing the two-dimensional phase plane includes: using the interlimb energy transfer ratio as the horizontal axis coordinate, multiplying the intralimb phase coupling delay by a preset scaling constant as the vertical axis coordinate, constructing a dimensionless two-dimensional space with scale equilibrium as the two-dimensional phase plane; extracting the decoupling features of the current effective gait cycle, generating the current test state point in the two-dimensional phase plane; extracting the coordinates of the effective historical state points within a fixed-length historical sliding window before the current gait cycle, calculating the arithmetic mean, and generating the historical baseline centroid.

[0015] Furthermore, the process of analyzing the spatial evolution geometric relationship and calculating the output compensatory decay index includes: constructing a trajectory vector from the centroid of the historical baseline to the current test state point; constructing a convergence vector from the centroid of the historical baseline to a preset healthy ideal state point, with the horizontal axis coordinate of the healthy ideal state point set to one and the vertical axis coordinate set to zero; calculating the Euclidean modulus parameter of the trajectory vector and the cosine similarity parameter between the trajectory vector and the convergence vector; multiplying the Euclidean modulus parameter with an index adjustment term constructed based on the cosine similarity parameter to nonlinearly amplify the weights deviating from the ideal state evolution direction, thereby generating the compensatory decay index. This invention provides a method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation. It has the following beneficial effects: 1. This invention calibrates the entire foot landing event using a plantar pressure time series. During this period, zero-velocity updates are performed on the integral velocity of the inertial sensor to obtain the closure error, which is then inversely compensated to the thigh and lower leg segments using a linear time decay function. By using the physical state of stable foot contact with the ground as a reference, the low-frequency integral drift error caused by long-term monitoring by the inertial sensor is eliminated, resulting in higher numerical accuracy for the subsequent calculation of the total mechanical energy of the lower limb and the interlimb energy transfer ratio.

[0016] 2. This invention utilizes the local maximum absolute value of the impact gradient of the support phase to construct a dynamic weighted envelope, and introduces it into the weighted cross-correlation calculation of the angular velocity sequence. Compared to processing the entire gait cycle data with equal weight, this technique focuses the calculation on the force loading stage, effectively separating the intralimb phase coupling delay generated by the affected lower limb at the moment of weight-bearing, and reducing the interference of conventional swing data in the non-impact stage on the cross-correlation results.

[0017] 3. This invention constructs a two-dimensional phase plane with interlimb energy transfer ratio and intralimb phase coupling delay as coordinate axes, and constructs a spatial vector by combining the historical baseline centroid and the ideal healthy state point. A compensatory attenuation index is generated by calculating the Euclidean modulus and cosine similarity of the included angle between the vector corresponding to the current test state point. This technical feature transforms isolated dynamic and kinematic parameters into a unified coordinate reference system, enabling objective quantification of gait evolution during the patient's rehabilitation process through the direction and distance of movement of the state point within the plane. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the system architecture of the present invention; Figure 2 This is a schematic diagram of the method flow of the present invention; Figure 3 This is a schematic diagram of the spatial topology and timing synchronization principle for multimodal signal acquisition according to the present invention. Figure 4 This is a schematic diagram illustrating the principle of attitude space calculation and gravity component stripping in this invention. Figure 5 This is a schematic diagram illustrating the principle of dynamic event calibration and impact feature extraction of the present invention. Figure 6 A timing diagram is captured from the dynamic adaptive analysis window of this invention; Figure 7 This is a flowchart of the cross-sensor data constraint and feature decoupling calculation process of the present invention; Figure 8 This is a schematic diagram of the phase plane evolution analysis and evaluation index output principle of the present invention; Figure 9 The graph shows the comparison between the compensatory attenuation index and the traditional gait asymmetry index in the long-distance walking fatigue test of the present invention. (a) is a schematic diagram of the multi-parameter evolution compensatory attenuation index of the present invention, and (b) is a schematic diagram of the pure kinematic stride asymmetry index of the traditional scheme.

[0019] Among them, 10 is the data acquisition module; 11 is the inertial measurement unit; 12 is the plantar pressure sensor; 20 is the spatial calculation module; 30 is the event calibration module; 40 is the window capture module; 50 is the feature decoupling module; and 60 is the state evaluation module. Detailed Implementation

[0020] The technical solutions in 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.

[0021] See attached document Figure 1 This invention provides a dynamic gait parameter detection, evaluation, and analysis system for lower limb postoperative rehabilitation, the system comprising: The data acquisition module 10 is used to simultaneously acquire kinematic state data and plantar dynamics distribution data of the test subject during walking. The data acquisition module 10 includes inertial measurement units 11 fixed to the outside of both thighs and calves, and to the insteps or shoe uppers of both feet, and plantar pressure sensors 12 placed inside both shoes. The inertial measurement units 11 fixed to the insteps or shoe uppers of both feet are used in conjunction with the plantar pressure sensors 12 to perform zero-velocity constraint correction of the foot.

[0022] The spatial calculation module 20 is used to timestamp-align the signals output by the inertial measurement unit 11 and the plantar pressure sensor 12, and store the aligned data stream in a circular buffer. The spatial calculation module 20 updates the attitude quaternion of the inertial measurement unit 11 relative to the global navigation coordinate system, rotates the original three-dimensional acceleration sequence to the global navigation coordinate system, and removes the gravitational acceleration component to generate a linear acceleration sequence. Specifically, the spatial calculation module 20 performs attitude quaternion updates and gravitational component removal processing on the inertial measurement units 11 at the thigh, calf, and foot positions respectively to form linear acceleration sequences for the corresponding segments.

[0023] The event calibration module 30 is used to extract time-domain transition points from the plantar pressure time series and calibrate heel strike events, full plantar strike events, and toe-off events. The event calibration module 30 defines the time interval of the double support phase based on the heel strike event time of one lower limb and the toe-off event time of the contralateral lower limb. The event calibration module 30 calculates the first-order time derivative of the plantar pressure time series of the affected side to generate a support phase impact gradient sequence.

[0024] The window truncation module 40 is used to obtain the moment when the support phase impact gradient sequence crosses the noise floor threshold as the starting point of the physical window, and to obtain the moment when the first zero-crossing point after the peak falls back as the ending point of the physical window. The window truncation module 40 delineates the core integration window based on the above two points. The window truncation module 40, in conjunction with the set maximum phase delay constant, extracts a data extraction buffer window with a time range greater than the core integration window from the annular buffer.

[0025] The feature decoupling module 50 performs zero-velocity updates on the foot inertial measurement unit 11 during the duration of the full-foot impact event, and corrects the drift of the motion integral results of the ipsilateral thigh, lower leg, and foot segments based on the closure error obtained from the foot zero-velocity update. The feature decoupling module 50 calculates the interlimb energy transfer ratio using the rate of change of mechanical energy during the double-support phase. The feature decoupling module 50 normalizes the scalar of the support phase impact gradient within the core integration window to construct a dynamic weighted envelope, and performs dynamic weighted cross-correlation calculations on the angular velocity sequences of the affected thigh and lower leg segments within the data extraction buffer window, extracting the corresponding time lag as the intralimb phase coupling delay. This intralimb phase coupling delay is used to characterize the degree of motion phase coordination of adjacent knee segments during the ground impact loading phase.

[0026] The state assessment module 60 is used to construct a two-dimensional phase plane with the interlimb energy transfer ratio as the horizontal axis and the intralimb phase coupling delay as the vertical axis, and to map the current state coordinate point and the historical centroid baseline coordinate point in the phase plane. The state assessment module 60 constructs a trajectory vector pointing to the current state coordinate point and a convergence vector pointing to the health reference point, calculates the magnitude and cosine of the included angle of the trajectory vector, and outputs the compensatory attenuation index after processing by the mapping function.

[0027] See attached document Figure 2 This invention provides a method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation, the method comprising the following steps: S10, synchronously acquire the three-dimensional acceleration sequence, angular velocity sequence and plantar pressure time sequence of the tester's gait cycle; S20, perform time alignment and circular buffering on the acquired sequence data, update the attitude quaternion to transform the three-dimensional acceleration sequence to the global navigation coordinate system, and strip the gravitational acceleration component to generate a linear acceleration sequence. S30 uses the jump points of the plantar pressure time series to identify plantar kinematic events to define the double support phase, and performs first-order derivative operation on the plantar pressure time series of the affected side to generate the support phase impact gradient sequence. S40, the core integration window is defined based on the waveform abrupt change points and zero crossing points of the supporting phase impact gradient sequence, and the data extraction buffer window is extracted in the annular buffer based on the set phase delay constant. S50 triggers the zero velocity update of the foot inertial measurement unit based on the whole foot landing event, and uses the zero velocity update result to correct the drift of the ipsilateral lower limb segment motion integral result in order to solve the interlimb energy transfer ratio during the double support phase. It constructs a dynamic weight envelope using the support phase impact gradient sequence in the core integration window, and performs dynamic weighted cross-correlation calculation on the angular velocity sequence of the affected thigh segment and the angular velocity sequence of the lower leg segment in the data extraction buffer window to obtain the intralimb phase coupling delay. S60 maps evolution coordinate points in a two-dimensional phase plane composed of interlimb energy transfer ratio and intralimb phase coupling delay through state evaluation module 60, analyzes the spatial geometric relationship between the current trajectory vector and the convergence target vector, and calculates the output compensatory decay index.

[0028] To further clarify the implementation of each technical aspect of the present invention, the following will provide a detailed description of the implementation of each functional module involved above and its internal processing flow.

[0029] See attached document Figure 3 The data acquisition module 10 is used to synchronously acquire the three-dimensional acceleration sequence, angular velocity sequence, and plantar pressure time sequence during the test subject's gait cycle. Specifically, this step may include the following sub-steps.

[0030] S101, a spatial deployment reference for configuring the inertial measurement unit 11 and the plantar pressure sensor 12.

[0031] The inertial measurement unit 11 is fixed to the outer center of the thigh, the outer center of the calf, and the instep or shoe upper of both lower limbs of the test subject. The inertial measurement unit 11 at the instep or shoe upper is used to acquire the foot motion state and provide physical constraints for zero velocity updates during full-foot landing events.

[0032] In practical applications, a unified spatial reference is the physical prerequisite for subsequent phase coupling analysis of the thigh segment angular velocity sequence and the lower leg segment angular velocity sequence. During the binding process, the axes of the sensor body coordinate system of the inertial measurement unit 11 are set to be parallel to the anatomical axes of the human lower limb. For the left and right lower limbs, since the lateral directions are opposite, the system performs a unified anatomical coordinate transformation on the outputs of the two inertial measurement units 11 during the data preprocessing stage, so that the forward axis, proximal axis, and lateral axis of the two lower limbs have a consistent sign and direction under the unified coordinate convention.

[0033] Specifically, the X-axis of the sensor's body coordinate system is set to point directly forward, the Y-axis to point upward along the long axis of the lower limb, and the Z-axis to point outward. To minimize relative slippage errors that may occur during the test subject's gait impact, the inertial measurement unit 11 is typically fixed with a non-elastic fabric strap. The plantar pressure sensor 12 is placed flat inside both shoes of the test subject.

[0034] The plantar pressure sensor 12 consists of a thin-film pressure sensing array whose physical distribution continuously covers the heel region, arch region, forefoot region, and toe region. This topology is used to comprehensively capture the dynamic load distribution across the entire plantar region during the support phase.

[0035] S102, Set the global sampling parameters of the data acquisition module 10.

[0036] To ensure complete capture of transient impact signals and numerical accuracy in subsequent discrete integration calculations, the system configures a uniform sampling frequency for all inertial measurement units 11 and plantar pressure sensors 12. The system's global discrete sampling frequency is set to [value missing]. Discrete time step between adjacent sampling points satisfy .

[0037] Based on the Nyquist sampling theorem and the frequency characteristics of human gait, in practical implementation, the global discrete sampling frequency can be... The sampling frequency should be configured within the range of 100Hz to 500Hz. If the sampling frequency is too low, the high-frequency impact waveform at the moment the heel touches the ground may be lost; if the sampling frequency is too high, unnecessary processor power consumption and memory overhead will be increased.

[0038] The system at each discrete time step The next step triggers a multimodal data sampling. Regarding the data transmission communication architecture between the data acquisition module 10 and the external microprocessor, those skilled in the art can implement it using standard I2C bus, SPI bus, or low-power wireless radio frequency communication protocols. The communication hardware pin connections and underlying register read / write mechanisms are well-known technologies in the field and will not be elaborated upon here.

[0039] S103 performs sampling frequency matching and time axis alignment operations for multi-source signals.

[0040] Since subsequent algorithms rely on dynamic characteristics to weight the kinematic cross-correlation sequence, if the two sensors are not synchronized, the calculated phase delay will lose its practical physical meaning. The data acquisition module 10 is configured with a hardware common clock source as a global master metronome. At each discrete time step... On the rising edge of the clock signal, the hardware common clock source synchronously sends an interrupt trigger command to the inertial measurement unit 11 and the plantar pressure sensor 12. The system eliminates the clock accumulation drift that occurs when multiple sensors operate independently through a hardware interrupt synchronization mechanism.

[0041] At any discrete time point At this point, the inertial measurement unit 11 synchronously outputs the current three-dimensional raw acceleration vector. With angular velocity vector The inertial measurement unit 11 includes at least two thigh inertial measurement units, two lower leg inertial measurement units, and two foot inertial measurement units; the thigh and lower leg inertial measurement units are used to characterize the kinematic coupling state of adjacent segments of the knee joint, and the foot inertial measurement units are used to characterize the velocity constraint state of the foot in the support state.

[0042] Three-dimensional primitive acceleration vector Includes the acceleration components along the X, Y, and Z axes in the sensor's body coordinate system. Angular velocity vector. It includes the angular velocity components that rotate around the three axes mentioned above.

[0043] At the same discrete time point The plantar pressure sensor 12 synchronously outputs the corresponding plantar pressure time series. .

[0044] Specifically, assuming the plantar pressure sensor 12 is built into The discrete pressure sensor array unit, the first Each unit in The local pressure value output at any time .when When the local normal force value is calibrated, the system directly sums the outputs of each array element; when When the pressure value is localized, the system first multiplies it by the effective force-bearing area of ​​the corresponding array element to convert it into a local normal force value, and then performs a summation operation. This results in the plantar pressure time series. We obtain the result by summing the following formulas: ; In the formula, The total number of array cells. This is the spatial distribution index of the sensor array units. The summation value directly represents the magnitude of the overall vertical reaction force generated at the interface between the foot and the ground.

[0045] After acquiring the aforementioned synchronous raw data, in order to suppress the thermoelectric noise of the sensor itself and the high-frequency mechanical vibration interference of human gait, the system performs preliminary alignment of the three-dimensional raw acceleration vector. Angular velocity vector and plantar pressure time series Digital low-pass filtering is applied. Considering that the fundamental frequency of active leg movements during normal walking is low, while the compensatory tremor on the affected side contains some high-frequency components, the cutoff frequency of the digital low-pass filter can be set between 15Hz and 30Hz.

[0046] For denoising and filtering of digital signals, those skilled in the art can use a second-order Butterworth low-pass filter or a moving average filtering algorithm. The design of the filter transfer function and the logic of the difference equation implementation are well-known technologies in the field and will not be elaborated here. After the above deployment and synchronization operations, the data acquisition module 10 generates a digital multimodal time series with a strict temporal mapping relationship, providing underlying data support for subsequent spatial coordinate system calculation and gait feature extraction.

[0047] See attached document Figure 4The spatial calculation module 20 performs temporal alignment and circular buffering on the acquired sequence data, updates the attitude quaternions, and transforms the original three-dimensional acceleration sequence to the global navigation coordinate system, thereby stripping the gravitational acceleration component to generate a linear acceleration sequence. The specific implementation of this step further includes the following sub-steps.

[0048] S201 performs circular buffer storage management for continuous multimodal time series.

[0049] Subsequent steps involving weighted cross-correlation calculations require time-shifting operations with lead and lag on the time-series data. If a linear memory structure is used directly, the time-shifting operations can easily exceed the current data write boundary, causing array out-of-bounds errors and program crashes.

[0050] The spatial calculation module 20 allocates a fixed-length circular buffer block in the system memory. The capacity of the circular buffer is set to... The minimum limit of this capacity must cover the duration of a complete gait cycle and a bidirectional time shift margin to ensure uninterrupted data flow within the retrieval range. In a specific embodiment, this capacity is determined based on the global discrete sampling frequency. The preset maximum estimated duration of the gait cycle With the preset maximum phase delay constant They are jointly determined to satisfy the following relation: ; The spatial calculation module 20 pushes the synchronized and aligned original 3D acceleration sequence, angular velocity sequence, and plantar pressure time sequence onto the circular buffer block in ascending order of time index. When the write pointer reaches the end of physical memory, it automatically wraps back to the beginning address of the buffer block and overwrites the oldest historical data. This rolling storage mechanism ensures that a continuous data field containing the current moment, historical moments, and extended moments is always maintained in the system memory, providing safe storage support for subsequent adaptive window truncation.

[0051] S202 updates the sensor attitude quaternion in real time based on multi-sensor data fusion. The inertial measurement unit 11 will undergo attitude deflection during motion, and its output three-dimensional raw acceleration vector is always relative to its own sensor body coordinate system.

[0052] The system establishes a global navigation coordinate system, with its Z-axis parallel to the direction of gravity and perpendicularly upward. Arbitrary discrete time points are defined. The attitude quaternion of the time sensor body coordinate system relative to the global navigation coordinate system is: . It consists of one scalar part and three vector parts, represented as The space solution module 20 utilizes angular velocity vectors. The attitude is updated using a first-order Runge-Kutta integral, and the discretized quaternion differential equation is expressed as follows: ; In the formula, It is the transpose operator. For discrete time steps, For the previous moment The posture quaternion, This represents quaternion multiplication. This is to expand the angular velocity vector into a pure quaternion with a real part of zero.

[0053] Since relying solely on the integral of angular velocity will result in cumulative drift error over time, the system introduces gravity direction observations from the original three-dimensional acceleration vector to perform gradient descent compensation for the quaternions. For complementary filtering algorithms or Kalman filtering algorithms based on the fusion of accelerometer and gyroscope data, those skilled in the art can use the standard Mahony attitude calculation algorithm to perform real-time correction of the quaternions. The design of its PI error compensator and gradient optimization process are well-known techniques in the field and will not be elaborated here.

[0054] S203 performs coordinate system mapping transformation and explicitly removes the gravitational acceleration component. During human walking, the inertial measurement unit 11 outputs the three-dimensional raw acceleration vector. In reality, it is the vector superposition of the human body's translational acceleration and the gravitational reaction force.

[0055] If the displacement is directly calculated by performing a quadratic integral on the original data, the gravity component that changes dynamically with the attitude will cause the integral result to diverge as a quadratic function series, resulting in the calculated mechanical energy lacking accurate physical reference meaning.

[0056] From a physical standpoint, gravity is a constant downward stationary vector in the global navigation coordinate system. However, the acceleration signal output by the inertial measurement unit (IMU) includes the projection of the gravity vector introduced by the sensor's own tilt. Therefore, it is necessary to first unify the measurement data to the global navigation coordinate system before performing vector subtraction.

[0057] Spatial solution module 20 uses the attitude quaternions updated in sub-step S202. The three-dimensional original acceleration vector in the sensor body coordinate system The rotation is mapped to the global navigation coordinate system. The 3D vector portion of the rotated pure quaternion is extracted, and then the gravitational acceleration constant vector in the global coordinate system is directly subtracted using vector subtraction. The calculation formulas for the above mapping and stripping operations are as follows: ; In the formula, This is the linear acceleration sequence obtained after solving, which purely represents the translational motion of the lower limbs. For attitude quaternions The conjugate quaternion is represented as . This represents the vector extraction operator that extracts the last three terms of a quaternion. The gravitational acceleration vector in the global navigation coordinate system is defined as follows: , is the gravitational constant.

[0058] After this step, the spatial solution module 20 can remove the gravity component from the kinematic integral link and generate a linear acceleration sequence with pure dynamic meaning, providing a data basis for subsequent calculation of the interlimb energy transfer ratio under the foot zero velocity update constraint and segmental drift correction constraint.

[0059] See attached document Figure 5 The event labeling module 30 uses the jump characteristics of the plantar pressure time series to label plantar kinematic events to define the double support phase, and performs a first-order derivative operation on the plantar pressure time series of the affected side to generate a support phase impact gradient sequence. The specific implementation of this step further includes the following sub-steps.

[0060] S301 uses the amplitude abrupt change of the plantar pressure time series to define the timestamp of the plantar kinematic event.

[0061] When a person walks, the contact state between the sole of the foot and the ground changes periodically. The event calibration module 30 presets a foot pressure noise threshold. and full contact pressure threshold Foot pressure noise threshold The threshold pressure is typically set to the average background noise of the plantar pressure sensor 12 under unloaded conditions plus a preset margin, or the peak background noise under unloaded conditions plus a preset margin, to distinguish between foot suspension and initial ground contact. Full contact pressure threshold. The parameters are set according to the proportion of the test subject's body weight to define the state in which the sole of the foot bears the full load of the body.

[0062] Specifically, regarding the full contact pressure threshold The specific value is determined by a static calibration procedure performed by the system before gait data acquisition. Specifically, the average steady-state pressure output from the plantar pressure sensor is collected when the test subject is standing still with one foot fully on the ground, and this value is recorded as the baseline reference value for body weight. Subsequently, the full contact pressure threshold was set. Set as the benchmark reference value for this weight A preset constant between 70% and 85%. Using this calibration ratio can effectively avoid localized load fluctuations when the heel first touches the ground, so as to accurately define the state of the sole bearing the full body load.

[0063] The event labeling module iterates through the plantar pressure time series 30 times. To improve the spatial specificity of foot kinematic event recognition, the event calibration module 30 divides the foot pressure sensor array into heel, arch, forefoot, and toe regions based on the physical distribution of the array, and calculates the pressure in the heel region separately. Pressure in the arch area Forefoot area pressure and pressure in the toe area .

[0064] The contact threshold for each region can be determined by adding a preset margin to the noise level of the corresponding region under unloaded conditions, or by a preset proportion of the average pressure of the corresponding region during single-foot static calibration. Total plantar pressure. For heel area pressure Pressure in the arch area Forefoot area pressure and pressure in the toe area sum.

[0065] The preset auxiliary contact condition can be pressure in the arch area. or pressure in the toe area Exceeding the contact threshold in the corresponding area can also cause pressure in the arch area. or pressure in the toe area Within a preset time window, the pressure distribution either continuously increases or remains in a non-zero contact state. The stable conditions for plantar pressure distribution include: total plantar pressure. Maintain the full contact pressure threshold within the preset stable time window. The above, and the pressure in the heel area Forefoot area pressure Furthermore, the pressure variation in the areas meeting the auxiliary contact conditions is less than the corresponding pressure fluctuation threshold.

[0066] When the current heel pressure value is greater than the noise floor threshold and the previous heel pressure value is less than or equal to the noise floor threshold, that is... and At that time, the event labeling module 30 will record the discrete time point. Defined as the moment of heel strike. .

[0067] When total plantar pressure Crossing the full contact pressure threshold from bottom to top And pressure in the heel area Pressure on the forefoot area All exceeded their respective area contact thresholds, resulting in pressure in the arch area. or pressure in the toe area When the preset auxiliary contact conditions are met and the above contact state is maintained within a preset time window, the corresponding crossing point time is defined as the start time of the full-foot strike event. By simultaneously constraining the total pressure threshold and multiple plantar zone pressure thresholds, it is possible to avoid misjudging a full plantar strike state based solely on local heel impact or forefoot force, and to reduce the probability of missing full plantar strike events caused by postoperative abnormal gait or insufficient toe force.

[0068] When the current pressure in the toe area is satisfied Less than or equal to the noise floor threshold, and total plantar pressure The pressure in the toe area decreased synchronously to near the noise threshold, and the previous moment's pressure... If the noise level is still greater than the noise floor threshold, the event calibration module 30 defines the discrete time point as the moment when the toes leave the ground. If abnormal gait after surgery results in weak pressure signals in the toe area, the event labeling module 30 can combine this with pressure signals in the forefoot area. pressure at the descent edge and toe area The descending edge and total plantar pressure The falling edge is used to determine the moment of toe-off event, thereby improving the stability of toe-off event recognition.

[0069] 302, the time interval of the double support phase across limbs.

[0070] In the gait cycle, the double support phase represents the period when one lower limb has just touched the ground while the other lower limb has not yet left the ground. This phase is the physical interval during which the body's center of gravity shifts and energy transfer occurs between limbs. The event labeling module 30 combines the kinematic event timestamps of both lower limbs to construct the time interval of the double support phase. .

[0071] In practice, the first time interval is formed by using the moment when the heel of the affected lower limb strikes the ground as the starting boundary and the moment when the toe of the unaffected lower limb leaves the ground as the ending boundary. Similarly, the second time interval is formed by using the moment when the heel of the unaffected lower limb strikes the ground to the moment when the toe of the affected lower limb leaves the ground.

[0072] The above time intervals are collectively referred to as the double support period intervals, denoted as... Where subscripts A and B represent two lower limbs that are opposite to each other. After the double support phase interval is formed, the event labeling module 30 performs a validity check on the interval; if Earlier than or equal to If the interval length is less than the preset minimum dual-support time threshold, then the matching of this group of plantar events is deemed invalid, and the matching process for the next group of heel-toe contact events and contralateral toe-off events is skipped. This interval provides a constraint boundary for the subsequent calculation of the energy transfer ratio.

[0073] S303 performs differential operations on the plantar pressure time series of the affected side to generate the support phase impact gradient sequence.

[0074] In the initial stage of lower limb surgery, patients experience a difference in the loading rate of the ground reaction force on the affected side due to pain compensation or insufficient muscle strength compared to healthy patients. This transient change in loading rate objectively reflects the dynamic cushioning capacity of the lower limb. The event labeling module 30 extracts the time series of plantar pressure on the affected side. Perform a first-order discrete difference operation on it in the time domain.

[0075] To suppress high-frequency numerical oscillations introduced by difference operations, a central difference algorithm is used to obtain the transient loading rate in the middle data segment of the sequence. The calculation formula is as follows: ; In the formula, Discrete time points The corresponding impact gradient value of the supporting phase, For discrete time steps, and These are the plantar pressure scalar values ​​of the adjacent sampling points before and after the current time.

[0076] For the first endpoint of the sequence ( (and tail endpoints, to avoid access) When an array goes out of bounds due to an illegal address, forward difference is used respectively. Boundary condition processing is performed with backward difference.

[0077] Because the original signal from the pressure sensor contains minute disturbances, direct differentiation may cause waveform glitches. After the differential operation, the event calibration module 30 additionally uses a sliding window with a width of... The moving average algorithm for Perform smoothing. Typically, the number of sampling points is 5 to 10.

[0078] Through the above traversal calculations, the event calibration module 30 converts the discrete pressure scalar data into a support phase impact gradient sequence characterizing the rate of force change. This gradient sequence is sensitive to dynamic anomalies at the moment of ground contact. The steepness of its waveform and the position of its peaks contain real compensatory information of the affected joint, providing parameters for subsequent delineation of the physical adaptive window and construction of dynamic weights.

[0079] See attached document Figure 6 The window extraction module 40 delineates the core integration window based on the waveform abrupt change points and zero-crossing points of the supporting phase impact gradient sequence, and extracts the data extraction buffer window from the annular buffer based on the set phase delay constant. The specific implementation of this step further includes the following sub-steps.

[0080] S401, adaptively truncates the core integral window based on the characteristics of the supporting phase impact gradient sequence.

[0081] The actual impact range of the lower limbs in the initial stage of ground contact is not a fixed static time period, but rather fluctuates dynamically with walking speed and the degree of joint compensation. From a biomechanical perspective, the impact of the early stage of the support phase is not a nonlinear, uniform loading process. By using impact gradient sequences, the core time period in which the lower limb muscle groups actively eccentrically contract to absorb the impact energy from the ground can be accurately located, thereby eliminating interference data from the subsequent stable support phase.

[0082] Window capture module 40 preset gradient noise floor threshold In practice, this threshold can be determined by collecting the root mean square error of the impact gradient sequence during the no-load period of the swing phase and multiplying it by an engineering margin factor of 2 to 3.

[0083] Window capture module 40 captures the moment of impact of the affected heel. Nearby backward traversal of support phase impact gradient sequence .

[0084] When it is determined that the gradient value of the current sampling point crosses the gradient noise floor threshold from bottom to top, it satisfies... and At that time, the window capture module 40 records the discrete time point as the starting point of the physical window. .

[0085] To prevent the algorithm from getting stuck in an infinite search loop due to insufficient force on the sensor or slippage, the window capture module 40 sets a maximum search time step limit. If no suitable crossing point is found after crossing the maximum limit, the fault tolerance mechanism is activated, and the process is directly reversed. It is marked as the starting point of the physical window.

[0086] Starting from the beginning of this physical window, the impact gradient sequence of the supporting phase gradually climbs to its peak and then begins to fall back. The window capture module 40 continues to search along the time axis and obtains the first zero-crossing point after the waveform peaks, where the value changes from positive to negative.

[0087] When satisfied and At that time, record the discrete time point as the end point of the physical window. To avoid the algorithm logic getting stuck in a dead zone, such as due to the hysteresis effect of the sensor piezoresistive material or severe pathological gait causing the sequence to fail to show a zero-crossing drop for a long period of time, the window truncation module 40 is internally configured with a maximum truncation duration. As a fault-tolerant constraint.

[0088] If we start counting from the beginning of the physical window, the traversal time span exceeds... If the zero-point condition has not yet been triggered, then force... Set as the physical window endpoint. The value of this parameter is typically set to be between 15% and 25% of the total duration of a single gait cycle. To achieve this dynamic constraint, the total duration of a single gait cycle can be estimated in real time by statistically analyzing the average time span of three consecutive normal gait cycles in the test subject's history, thereby adapting to changes in walking speed.

[0089] Based on the two time boundaries obtained through the optimization process described above, the window extraction module 40 delineates a core integral window representing the actual force impact range, denoted as: ,in This is the end point of the physical window. If Earlier than or equal to If the core integration window length is less than the preset minimum window length, the window truncation module 40 determines that the core integration window is invalid and abandons the feature decoupling calculation of the current gait cycle.

[0090] S402 extends the data extraction buffer window within the circular buffer.

[0091] In subsequent feature decoupling and cross-correlation calculations, the system needs to use the data within the core integration window as a base template to perform a time-shift comparison of the angular velocity sequence sliding with lead and lag. If the data array is extracted directly according to the length of the core integration window, edge data indexes may be lost due to out-of-bounds errors during time-shift stepping.

[0092] Window capture module 40 calls the system's preset maximum phase delay constant. This constant represents the maximum expected neuromuscular lag time in the affected limb post-surgery due to compensatory pain or muscle weakness. Based on clinical statistics from lower limb rehabilitation, The specific value range is usually set between 0.1 seconds and 0.3 seconds.

[0093] Window capture module 40 captures the core points window. The maximum phase delay constant is extended on both sides of the time axis. The start time boundary of the data extraction buffer window was calculated. and termination time boundary Before performing memory data extraction, the window capture module 40 performs a validity check on the extended time boundary.

[0094] If the calculation yields If the timestamp of the data is earlier than the oldest data currently remaining in the circular buffer, the starting boundary will be forcibly clamped to the time point corresponding to the oldest valid data in the buffer. In real-time online computing mode, when the physical window terminates... Once determined, the window capture module 40 delays the maximum phase delay constant. After the corresponding time period, perform the dynamic weighted cross-correlation calculation to ensure the data extraction buffer window. All data has been written to the circular buffer. If the system is in low-latency, fast output mode, or if data transmission is interrupted... If the latest write pointer is exceeded, the termination boundary will be clamped to the current latest valid data time point, and the time shift term will be dynamically reduced according to the actual time range of the data already written in the circular buffer. The search scope.

[0095] Using the aforementioned time boundary constraints that have undergone validity verification, the window truncation module 40 returns the circular buffer block established in step S20. Through memory pointer address mapping, it directly extracts the multimodal time series block covering the entire extended interval. This extended truncation operation constructs a safety boundary at the physical memory level, avoiding the risk of array out-of-bounds crashes caused by time-domain shift operations, and ensuring the execution integrity of the underlying cross-correlation matching algorithm.

[0096] See attached document Figure 7 The feature decoupling module 50 triggers the zero-velocity update of the foot inertial measurement unit 11 using the plantar contact event, calculates the interlimb energy transfer ratio based on the anti-drift result, and constructs a dynamic weight envelope to extract the physically constrained intralimb phase coupling delay. A specific implementation of this step further includes the following sub-steps.

[0097] S501 utilizes a zero-speed update mechanism triggered by a full-foot landing event.

[0098] Due to the inherent white noise and low-frequency flicker noise of microelectromechanical system (MEMS) sensors, directly integrating a linear acceleration sequence quadratically will result in a displacement accumulation error that diverges quadratically over time. The feature decoupling module 50 introduces physical topological constraints based on human gait to perform anti-drift calculations.

[0099] Based on the start time of the full-foot strike event defined in step S30 When the sole of the foot is in full contact and the above-mentioned stable conditions for sole pressure distribution are met, the macroscopic translational velocity of the foot inertial measurement unit 11 fixed at the instep or shoe upper relative to the ground is physically close to zero.

[0100] The feature decoupling module 50 intervenes in the integral link during the full-foot contact period, performs zero-velocity update on the integral velocity state of the foot inertial measurement unit 11, and corrects the vertical displacement drift of the foot as a constraint error of vertical displacement change.

[0101] For the inertial measurement unit 11 fixed at the outer center position of the thigh and the outer center position of the lower leg, the feature decoupling module 50 does not directly force its segment velocity to zero, but uses the closure error obtained by the foot inertial measurement unit 11 during the full foot contact as a drift correction reference for the same side gait interval, and performs continuous error compensation for the integral velocity and vertical displacement of the thigh segment and the lower leg segment.

[0102] In the process of drift correction for the thigh and lower leg segments, the foot velocity closure error is not directly equated to the actual velocity error of the corresponding segment. Instead, the zero velocity update result of the foot is used as the drift constraint reference for the gait cycle on the same side, and continuous error compensation is performed by combining the time synchronization integration results of the inertial measurement units 11 of each segment.

[0103] To avoid a direct zeroing of the time-domain curve that could cause a step abrupt change and lead to infinite dead zone glitches in subsequent mechanical energy derivative calculations, the feature decoupling module 50 extracts the velocity closure error vector and vertical displacement closure error of the foot inertial measurement unit 11 at the trigger moment of the full-foot landing event. This error is then inversely distributed and compensated for using a linear time decay function, and integrated into the integral result of the entire preceding swing phase, thus ensuring the continuous differentiability of the velocity and displacement curves. The thigh and lower leg inertial measurement units 11 are not directly zeroed during the full-foot landing period; their integral results are drift-corrected based on the closure error of the same side foot and the time synchronization relationship of the corresponding segments.

[0104] For the construction of specific zero-velocity triggered Kalman filter observation equations, those skilled in the art can use the standard ZUPT algorithm architecture for processing. The design of its state-space equations and the update process of the error covariance matrix are well-known techniques in this field and will not be elaborated here.

[0105] S502, based on the anti-drift correction results, calculates the integral and ratio of the interlimb energy transfer ratio.

[0106] Gait energy alternation is concentrated during the double support phase. The feature decoupling module 50 uses the corrected linear velocity and vertical displacement, combined with human segment mass parameters, to calculate the mechanical energy of one lower limb. Specifically, the mechanical energy of one lower limb can be obtained by superimposing the mechanical energy of the ipsilateral thigh segment, lower leg segment, and foot segment.

[0107] Specifically, the velocity and vertical displacement of the thigh segment are obtained by the thigh inertial measurement unit after attitude calculation, gravity stripping, integration and drift correction. The velocity and vertical displacement of the lower leg segment are obtained by the lower leg inertial measurement unit through the same process. The velocity and vertical displacement of the foot segment are obtained by the foot inertial measurement unit after zero velocity update correction.

[0108] Suppose that one lower limb includes the thigh segment, the lower leg segment, and the foot segment, and the segment index is denoted as . Quality of each segment The weight distribution was determined based on anthropometric proportions derived from the test subject's weight, with the affected and unaffected lower limbs using the same weight distribution rules. For any segment... Its discrete time points The velocity modulus is denoted as The vertical displacement is denoted as The feature decoupling module 50 first calculates the segment in Segmental mechanical energy at time Then, the segmental mechanical energies of the thigh segment, lower leg segment, and foot segment are summed to obtain the unilateral lower limb at discrete time points. Total mechanical energy The calculation formula is: ; In the formula, is the gravitational acceleration constant.

[0109] The feature decoupling module 50 calculates the rate of change of mechanical energy between adjacent sampling points. When using the segmental mechanical energy superposition method, the feature decoupling module 50 first calculates the segmental mechanical energy of the ipsilateral thigh segment, lower leg segment, and foot segment at each sampling time, and then sums them to obtain the total mechanical energy of the lower limb on that side. Then, the rate of change of mechanical energy of the lower limb at the corresponding sampling time is calculated based on the difference in total mechanical energy of the lower limb between adjacent sampling points. The affected lower limb and the healthy lower limb are treated with the same segment division, mass allocation rules, and mechanical energy change rate calculation method.

[0110] The dual support period interval defined in step S30 Within this process, the absolute values ​​of the rate of change of mechanical energy of the affected lower limb (labeled A) and the unaffected lower limb (labeled B) are integrated over time. Finally, the interlimb energy transfer ratio is calculated. The formula is expressed as: ; This ratio objectively quantifies the difference in compensation between the affected limb and the healthy limb in absorbing and transmitting ground reaction forces during the alternating foot support phase.

[0111] S503 constructs a dynamic weight envelope based on impact gradient normalization.

[0112] Feature decoupling module 50 extracts the core integral window extracted in step S40. Internal support phase impact gradient sequence To ensure the non-negativity of the dynamic weight envelope, in a preferred embodiment, the feature decoupling module 50 performs positive loading screening on the support phase impact gradient sequence, setting gradient values ​​less than zero to zero; in other embodiments, absolute value processing can also be performed on the support phase impact gradient sequence.

[0113] The feature decoupling module iterates 50 times to obtain the local maximum absolute value of the gradient within the window. By using this local maximum absolute value, a division normalization operation is performed on each nonnegative gradient scalar within the core integration window, generating a dimensionless weighted function sequence. .

[0114] To ensure the integrity of the algorithm logic and avoid errors caused by the tester being in a static position or the sensor experiencing minimal force. Approaching zero, thus triggering a program crash dead zone when dividing by zero, the feature decoupling module 50 introduces a positive floating-point constant into the denominator. The constant The value is determined based on the precision setting of the system data type, and is usually set to 10. -6 Order of magnitude. In normalization operations, This represents the nonnegative impact gradient of the support phase after positive loading screening or absolute value processing. The normalized mathematical transformation formula is: ; This transformation generates a scalar envelope with values ​​constant between 0 and 1. Since this scalar envelope is obtained by normalizing the nonnegative impact gradient of the support phase, it avoids negative gradient values ​​from negatively weakening the weight contribution of the high-impact phase when participating in cross-correlation calculations. This envelope shape refines the dynamic distribution weight of the ground contact impact load in the time domain.

[0115] S504 performs dynamic weighted cross-correlation to extract intralimb phase coupling delay.

[0116] Conventional kinematic cross-correlation directly performs equal-weight sliding multiplication on two angular velocity sequences, which may mask the local joint motion lag caused by the high-load impact at the moment of gait landing. The feature decoupling module 50 encapsulates the generated dynamic weights. It is embedded as an integral kernel in the cross-correlation analysis of the angular velocity on the affected side.

[0117] The feature decoupling module 50 extracts the coronal axis rotation component orthogonal to the sagittal plane of the human body from the three-dimensional angular velocity sensor data of the affected side, and forms angular velocity scalar sequences for the thigh side respectively. Scalar sequence of angular velocities on the lower leg side The computational closed loop is constructed as follows: ; In the formula, For time translation terms The weighted cross-correlation coefficient is calculated. To ensure that the data flow for correlation calculation is closed-loop, the translation term is used. The traversal range is strictly limited to the range set in step S40. Between, and the translation step size is equal to the discrete time step size. .

[0118] To avoid the backward translation of the time axis If front-end data is missing or memory exceeds the limit, the data index of the specified shift sequence is extended using the safe extraction buffer window established in step S402.

[0119] Feature decoupling module 50 searches within the above traversal range for the weighted cross-correlation number. The translation term corresponding to the maximum value is the extracted intralimb phase coupling delay parameter. In this embodiment, The sign direction is defined according to the time shift direction in the cross-correlation calculation formula; when A value greater than zero indicates that the angular velocity sequence of the lower leg segment lags behind the angular velocity sequence of the upper leg segment. A value less than zero indicates that the angular velocity sequence of the lower leg segment leads the angular velocity sequence of the upper leg segment. If the maximum weighted cross-correlation coefficient is lower than the preset confidence threshold, the feature decoupling module 50 determines that the intralimb phase coupling delay of the current gait cycle is invalid.

[0120] The underlying computational logic forces the alignment of the kinematic phase within the dynamic deterioration range, and outputs parameters characterizing the degree of motion phase lag or neuromuscular control lag generated by the affected joint during the weight-bearing impact phase.

[0121] See attached document Figure 8 The status assessment module 60 constructs a two-dimensional feature phase plane, uses geometric vectors to analyze the spatiotemporal evolution trend of feature points, and finally outputs a compensatory attenuation index characterizing the rehabilitation process or gait deterioration. The specific implementation of this step further includes the following sub-steps.

[0122] S601, Construction of two-dimensional phase plane baseline and coordinate mapping.

[0123] The state assessment module 60 establishes a two-dimensional Euclidean geometric space. Considering the energy transfer ratio... With phase delay The physical dimensions and numerical magnitudes are inconsistent, and directly constructing orthogonal axes would lead to geometric distance calculations being dominated by a single dimension. The state assessment module 60 uses a preset dimensional scaling factor to perform coordinate normalization mapping. (Setting the horizontal axis coordinates...) Since it is a dimensionless ratio, it can be directly mapped. Set the ordinate. , where scaling constant This is used to linearly map the numerical range of time delay to an order of magnitude similar to the energy ratio on the horizontal axis, based on the characteristic that the phase delay of human gait is usually in the range of tens to hundreds of milliseconds. The value is typically taken to be 5 to 10, thus constructing a scale-balanced dimensionless two-dimensional geometric space. In practical implementation, Use seconds as the input unit. A scaling unit that is the inverse of the second is used to ensure that the mapped vertical axis coordinates are dimensionless values.

[0124] The state evaluation module 60 extracts the decoupling features of the current gait cycle and generates the current test state point in the phase plane. Specifically, when the current gait cycle is a valid gait cycle, at least the following conditions must be met: heel strike event, full plantar strike event, and toe lift event must all be successfully calibrated; the dual support phase must pass the validity check; and the core integral window length must meet the preset minimum window length.

[0125] To capture the test subject's own gradual baseline during continuous walking, the state evaluation module 60 maintains a length of A sliding window showing the history of the state, in which Typically, the value is taken as 10 to 15 gait cycles. The system extracts the coordinates of all historical state points cached within the window prior to the current gait cycle and calculates the arithmetic mean.

[0126] To avoid memory dead zones due to insufficient data in the early stages of testing, the cumulative number of steps during system startup should be less than [a certain threshold]. When the system has not yet formed any historical state points, the state assessment module 60 automatically switches to using the total number of historical gait data collected as the denominator to calculate the dynamic mean. When the system has not yet formed any historical state points, the state assessment module 60 uses the state point of the first valid gait cycle as the initial historical baseline and does not output the compensatory attenuation index corresponding to that first valid gait cycle. From the second valid gait cycle onwards, the historical baseline centroid is preferably calculated from the valid historical state points before the current gait cycle to avoid the current state point participating in its own baseline calculation and weakening the deviation representation. This is used as the origin coordinate to record the historical baseline centroid calculated by sliding. .

[0127] S602, Geometric analytical calculation of spatial evolution vectors.

[0128] The state assessment module 60 utilizes vector subtraction operations in spatial analytic geometry to construct an evolutionary trend representation. It constructs trajectory vectors representing the actual deviations from the trajectory. Its direction is from the centroid of the historical baseline to the current test state point, that is This vector represents the direction and magnitude of the transient deviation of the current gait characteristics from the test subject's recent average state.

[0129] Meanwhile, the status assessment module 60 presets ideal health status points. In biomechanical principles, the ideal reference gait assumes symmetrical energy work between limbs, i.e., the horizontal axis... Furthermore, there is no abnormal phase lag caused by causal compensation, i.e., the vertical axis... .

[0130] Using the above definitions, the state evaluation module 60 constructs a convergence vector representing the ideal recovery path. Its direction points from the historical baseline centroid to the ideal health state point, and the calculation process is as follows: When the historical baseline centroid relative to the ideal state of health When the distance between them is less than the preset minimum threshold, the state evaluation module 60 determines that the current historical baseline is close to the ideal healthy state, and clamps the magnitude of the convergence vector to the minimum threshold to avoid numerical instability in the subsequent calculation of the cosine of the included angle.

[0131] S603, the structured mapping output of the compensatory attenuation index.

[0132] The state assessment module 60 comprehensively analyzes the geometric relationship between the trajectory vector and the convergence vector. The system calculates the Euclidean magnitude of the trajectory vector. The modulus parameter quantifies the degree of gait fluctuation. Subsequently, the cosine similarity parameter between the trajectory vector and the convergence vector is calculated. : ; In the formula, This is a preset positive floating-point constant for the system, typically with a value of 10. -6 This is used to avoid the algorithm dead zone where the denominator is zero when the centroid of the historical baseline approaches the ideal state point.

[0133] In a physical sense, when When this occurs, it indicates that the current compensatory state is evolving towards an ideal healthy state; when... This indicates that gait control ability is deviating from a healthy trajectory and is trending towards deterioration.

[0134] The state evaluation module 60 inputs the obtained modulus parameter and the cosine similarity parameter of the included angle into a preset mapping function, and finally generates the compensatory attenuation index value. The specific structured mapping logic takes the following form: ; In the formula, The base of the natural logarithm, The empirical penalty coefficient for adjusting the sensitivity to the decay trend is usually calibrated between 1.5 and 2.5 based on clinical experience.

[0135] The physical derivation logic of this structured mapping is that the simple gait fluctuation modulus cannot distinguish between good and bad evolution. By introducing an exponential adjustment term containing angle cosine, when the gait evolution direction points to the ideal state ( When the overall exponential function domain contracts, it indicates a benign restorative fluctuation; when the evolutionary direction deviates from the ideal state ( When the exponential function is used, it will nonlinearly amplify the weight of the current gait fluctuation.

[0136] Through this mathematical transformation logic, the output compensatory attenuation index provides an assessment indicator for the degree of lower limb fatigue accumulation or neuromuscular compensatory decline in clinical rehabilitation monitoring. In specific implementation, the state assessment module 60 can also perform moving average or median filtering on the compensatory attenuation index and set valid output conditions; when plantar events are missing, double support phase is illegal, core integration window is invalid, or intralimb phase coupling delay is determined to be invalid, the state assessment module 60 suspends the output of the compensatory attenuation index of the current gait cycle, or uses the assessment result of the previous valid gait cycle.

[0137] To enable those skilled in the art to more clearly understand the purpose, technical solution, and advantages of this invention, the present invention will be further described in detail below with reference to specific application embodiments, real experimental test data, and corresponding drawings. It should be noted that the embodiments described herein are only for explaining the present invention and are not intended to limit the scope of protection of the present invention.

[0138] I. Specific Application Examples Taking a patient who underwent right lower extremity anterior cruciate ligament reconstruction 3 months prior as an example, a continuous walking test on flat ground was performed. The system's global discrete sampling frequency was configured as follows: Hz, or discrete time step .

[0139] Within a certain test gait cycle, the event calibration module 30 determines the timing of the heel strike event on the affected side (right side). Timing of toe lift-off event on the unaffected side (left side) s, based on this, defines the double support period interval. s.

[0140] The feature decoupling module 50 integrates the rate of change of mechanical energy within this interval. The calculated absolute value of the integral of the mechanical energy change in the affected lower limb is 12.5 J, and in the healthy lower limb it is 15.0 J. Based on the algorithm, the interlimb energy transfer ratio for this gait cycle is generated. This value is lower than the ideal value of 1, indicating that there is a compensatory attenuation of about 16.7% in the energy buffering work on the affected side.

[0141] The window extraction module 40 searches for abrupt changes and zero-crossing points in the impact gradient sequence on the affected side, and extracts the core integral window. s. Extract the impact gradient sequence of the support phase within the window, perform local maximum normalization, and generate a dynamic weight envelope ranging from 0 to 1. .

[0142] Set the maximum phase delay constant. s, within the securely captured data extraction buffer window, the above The cross-correlation calculation of the angular velocities of the affected thigh and lower leg is embedded. When the time translation term... At time s, the weighted cross-correlation coefficient The peak value is reached. Therefore, the intralimb phase coupling delay of this gait cycle is extracted. s represents an abnormal neuromuscular retardation of 60 ms in the initial stage of loading on the knee joint.

[0143] Entering the state evaluation module 60, the dimensional scaling constant is called. Construct a two-dimensional dimensionless geometric space. Map the coordinates of the current test state point. ,in , .

[0144] Extract historical state sliding window data from the system cache, assuming the recent historical baseline centroid coordinates are... The preset ideal health state point is... Construct vectors respectively: Trajectory vector representing the true deviation from the trajectory Calculate its Euclidean modulus. .

[0145] The convergence vector representing the ideal recovery path: Its mold length .

[0146] Calculate the cosine similarity parameter of the angle between two vectors The vector dot product is 0.033 × 0.200 + 0.100 × 0.700 = 0.0766.

[0147] Substituting into the formula, we get This value is greater than zero and close to 1, indicating that the current evolution of gait fluctuations is highly consistent with the ideal rehabilitation path.

[0148] Substitute the experience penalty coefficient Calculate the compensatory attenuation index .

[0149] In contrast, if the patient's gait deteriorates during subsequent fatigue testing, the current state point will drift to... The new trajectory vector is (-0.050, 0.050), with a magnitude of... Recalculated The output of the compensatory decay index at this time is .

[0150] Numerical comparison mapping results show that although the absolute deviation distance of the deteriorating state (magnitude 0.071) is even smaller than that of the recovering state (magnitude 0.105), the exponential mapping function triggers an exponential nonlinear penalty due to the deviation of its evolution direction from the ideal node (i.e., the cosine of the angle changes from positive to negative). This amplifies the final output evaluation exponent by nearly 28 times (from 0.014 to 0.404). This mechanism improves the system's sensitivity to early detection of pathological compensatory deterioration trends.

[0151] To further verify the effectiveness and reliability of the present invention in real clinical scenarios, the following experimental verification and effect comparison analysis were conducted.

[0152] This experiment recruited 20 patients in the mid-recovery stage after unilateral lower limb orthopedic surgery (experimental group) and 20 healthy subjects with no history of lower limb disease (control group). The testing environment was a 100-meter barrier-free straight corridor in a temperature-controlled room. Subjects were asked to perform a long-distance back-and-forth walking fatigue test for 30 minutes at a comfortable pace.

[0153] The synchronous comparison method adopted the conventional kinematic analysis scheme that relies solely on the spatial step length asymmetry rate (affected side step length / healthy side step length) for judgment.

[0154] according to Figure 9 It can be seen that in the healthy control group, 30 minutes of continuous walking did not induce significant gait compensation, and the compensation attenuation index output by the system of this invention... The noise level remained consistently below 0.02, indicating an extremely low noise floor. This validates the effectiveness of the minimum bias introduced in this algorithm. Mechanisms such as the noise floor threshold successfully filtered out the risk of false alarms caused by normal physiological fluctuations.

[0155] In the postoperative experimental group, the impact cushioning capacity of the affected limb began to decline as walking time increased and lower limb muscle fatigue accumulated. Graph analysis showed that the traditional stride asymmetry index did not show significant statistical deterioration in the first 20 minutes of the test, and only crossed the abnormality threshold when the subject subjectively reported extreme fatigue and clinically visible slight limping (between the 22nd and 24th minutes of the test).

[0156] In contrast, the system of this invention, during the test interval of 14 to 16 minutes, triggered a nonlinear mapping penalty due to a slight increase in intralimb phase coupling delay and a negative angle between phase plane evolution vectors, resulting in a compensatory decay index. A steep upward step jump was observed. This data indicates that, by extracting kinematic features under dynamic weight constraints and fusing the evolutionary analysis of two-dimensional geometric vectors, this invention can advance the detection time window for the deterioration of compensatory fatigue in the lower limbs by approximately 6 to 8 minutes (equivalent to advancing 400 to 500 gait cycles) compared to traditional single-dimensional pure kinematic judgment methods. This advanced warning capability provides sufficient temporal margin for rehabilitation intervention, effectively avoiding secondary joint damage caused by excessive exercise, and fully demonstrating the technical advantages of the feature decoupling and multi-parameter fusion logic of this system.

[0157] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation, characterized in that, Includes the following steps: Simultaneously acquire the three-dimensional raw acceleration sequence, angular velocity sequence, and plantar pressure time sequence of the test subject's gait cycle; The original three-dimensional acceleration sequence is transformed into a spatial coordinate system and subjected to gravity stripping to generate a linear acceleration sequence. Simultaneously, the plantar pressure time series is used to calibrate plantar kinematic events to define the dual support phase interval, and a support phase impact gradient sequence is generated based on differential operations. The core integration window is defined based on the waveform characteristics of the supporting phase impact gradient sequence, and the data extraction buffer window is extended and truncated based on the set maximum phase delay constant. Based on the plantar kinematic events, drift correction is performed on the linear acceleration sequence to calculate the interlimb energy transfer ratio within the dual support phase interval; A dynamic weighted envelope is constructed using the impact gradient sequence of the support phase within the core integration window. Within the data extraction buffer window, the dynamic weighted envelope is used to perform dynamic weighted cross-correlation calculation on the angular velocity sequence to extract the intralimb phase coupling delay. A two-dimensional phase plane is constructed with the interlimb energy transfer ratio and the intralimb phase coupling delay as coordinate axes. The spatial evolution geometric relationship of the current test state point relative to the centroid of the historical baseline is analyzed in the two-dimensional phase plane, and the compensatory attenuation index is calculated and output.

2. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, Before synchronously acquiring the three-dimensional raw acceleration sequence, angular velocity sequence, and plantar pressure time sequence of the test subject's gait cycle, the process also includes: The inertial measurement unit was fixed to the center of the outer thigh, the center of the outer calf, and the instep or shoe upper of both lower limbs of the tester, respectively, and the plantar pressure sensor was placed inside the tester's shoes. A synchronization interrupt is triggered using a global common hardware clock source at a preset discrete time step to perform sampling frequency matching and time axis alignment between the inertial measurement unit and the plantar pressure sensor, and the aligned multimodal time series data is pushed onto the stack and written into a circular buffer.

3. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The method of using the plantar pressure time series to calibrate plantar kinematic events to define the dual support phase interval includes: The sole of the foot is divided into the heel region, arch region, forefoot region, and toe region, and the pressure time series of the corresponding region and the total plantar pressure time series are extracted respectively. When the pressure in the heel area exceeds the preset noise threshold from bottom to top, it is defined as the moment of heel contact with the ground; When the total plantar pressure time series crosses the full contact pressure threshold from bottom to top, and the plantar pressure in each region meets the preset plantar pressure distribution stability condition, it is defined as the start time of the full plantar contact event. When the pressure in the toe area drops to the noise threshold and the total plantar pressure time series decreases synchronously, it is defined as the moment when the toe leaves the ground; The double support period interval is formed by taking the heel-to-ground contact event of one lower limb as the starting boundary and the toe-to-ground departure event of the opposite lower limb as the ending boundary.

4. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The generation of the support phase impact gradient sequence based on differential operations includes: Extract the plantar pressure time series of the affected side and use the central difference algorithm to calculate the transient loading rate between adjacent sampling points; Boundary conditions are handled by forward difference and backward difference at the beginning and end of the sequence, respectively; The sequence after differentiation is smoothed using a moving average algorithm to generate the impact gradient sequence of the support phase.

5. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The step of defining a core integration window based on the waveform characteristics of the supporting phase impact gradient sequence, and extending and extracting a data extraction buffer window based on a set maximum phase delay constant, includes: The impact gradient sequence of the support phase is traversed backward from the moment of the heel strike on the affected side, and the discrete time point at which the gradient value crosses the gradient noise threshold from bottom to top is defined as the starting point of the physical window. Searching backward from the starting point of the physical window, the first discrete time point after the peak of the waveform to change from positive to negative is defined as the ending point of the physical window; The core integration window is formed by the physical window start point and the physical window end point; The maximum phase delay constant is extended to both sides of the time axis of the core integration window to form the data extraction buffer window.

6. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The step of performing drift correction on the linear acceleration sequence based on the plantar kinematic events to calculate the interlimb energy transfer ratio within the dual support phase interval includes: During the duration of the full-plantar landing event in the plantar kinematic event, the integral velocity state of the inertial measurement unit at the instep or upper position is updated to zero velocity, and the foot velocity closure error and vertical displacement closure error are obtained. The closure error is inversely distributed using a linear time decay function, and continuous error compensation is performed on the integral velocity and vertical displacement of the thigh and lower leg segments on the same side. The mechanical energy of each lower limb segment is calculated by combining the velocity modulus and vertical displacement of each segment after compensation. The total mechanical energy of the thigh segment, lower leg segment and foot segment is obtained by summing the mechanical energy of each segment. During the dual-support period, the absolute values ​​of the total mechanical energy change rates of the affected lower limb and the healthy lower limb are integrated over time and the ratio is calculated. The interlimb energy transfer ratio is then calculated and output.

7. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The construction of the dynamic weight envelope using the support phase impact gradient sequence within the core integral window includes: Extract the impact gradient sequence of the supporting phase within the core integral window and perform nonnegation processing; Obtain the local maximum absolute value of the gradient within the core integration window; The local maximum absolute value is used to perform a normalization operation on each nonnegative gradient scalar within the core integration window to generate the dynamic weight envelope with values ​​ranging from zero to one.

8. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The step of performing a dynamically weighted cross-correlation calculation on the angular velocity sequence using the dynamic weight envelope within the data extraction buffer window to extract the intralimb phase coupling delay includes: Extract the coronal axis rotation component orthogonal to the sagittal plane of the human body from the angular velocity sequence of the affected side to form the angular velocity scalar sequence of the thigh side and the angular velocity scalar sequence of the calf side; The product integral of the angular velocity scalar sequence of the thigh side, the angular velocity scalar sequence of the lower leg side after time translation, and the dynamic weight envelope at the corresponding time is calculated within the traversal range of the positive and negative maximum phase delay constants. Obtain the time shift term that causes the weighted cross-correlation coefficient to reach its maximum value, and use it as the intralimb phase coupling delay.

9. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The construction of a two-dimensional phase plane with the interlimb energy transfer ratio and the intralimb phase coupling delay as coordinate axes includes: Using the interlimb energy transfer ratio as the horizontal axis, and multiplying the intralimb phase coupling delay by a preset scaling constant as the vertical axis, a dimensionless two-dimensional space with a balanced scale is constructed as the two-dimensional phase plane. Extract the decoupling features of the current effective gait cycle and generate the current test state point in the two-dimensional phase plane; Extract the coordinates of valid historical state points within a fixed-length historical sliding window prior to the current gait cycle, calculate the arithmetic mean, and generate the historical baseline centroid.

10. The method for detecting, evaluating, and analyzing dynamic gait parameters in lower limb postoperative rehabilitation according to claim 1, characterized in that, The step of resolving the spatial evolution geometric relationship of the current test state point relative to the centroid of the historical baseline in the two-dimensional phase plane, and calculating and outputting the compensatory attenuation index includes: A trajectory vector is constructed by pointing from the centroid of the historical baseline to the current test state point; A convergence vector is constructed by pointing from the centroid of the historical baseline to a preset ideal health state point, wherein the horizontal axis coordinate of the ideal health state point is set to one and the vertical axis coordinate is set to zero. Calculate the Euclidean modulus parameter of the trajectory vector and the cosine similarity parameter of the angle between the trajectory vector and the convergence vector; The compensatory decay index is generated by multiplying the Euclidean modulus parameter with an exponential adjustment term constructed based on the cosine similarity parameter of the included angle, nonlinearly amplifying the weights that deviate from the ideal state evolution direction.

Images