A wearable sensor-driven multi-dimensional musculoskeletal digital twin modeling method

Abstract

The application discloses a wearable sensor-driven multi-dimensional musculoskeletal digital twin modeling method, which belongs to the fields of signal and information processing, computer application technology and biomechanics and the like. The application can form an end-to-end multi-dimensional musculoskeletal digital twin modeling closed loop from kinematics, dynamics to muscle activation. The whole body joint angle and three-dimensional ground reaction force are taken as unified input, the joint net torque is solved by combining reverse dynamics, and the muscle activation is solved by static optimization under the torque balance constraint, so that multi-dimensional evaluation output containing joint kinematics, three-dimensional ground reaction force, joint torque and muscle activation is formed. Meanwhile, the closed loop framework can be applied to walking, running and stepping and the like, and can provide an interpretable, generalizable and unconstrained multi-dimensional musculoskeletal digital twin modeling for sports training, clinical rehabilitation and mobile health monitoring and the like.

Description

Technical Field

[0001] This invention belongs to multiple interdisciplinary fields, including signal and information processing, computer application technology, and biomechanics. Specifically, it relates to a wearable sensor-driven multidimensional musculoskeletal digital twin modeling method. This method can collect kinematic and dynamic data of human movement within a fully wearable sensor framework, enabling multidimensional digital twin modeling of human actions in multiple scenarios. Background Technology

[0002] Biomechanical assessment of human movement is crucial for understanding motor function, monitoring health status, and disease screening and intervention. Multidimensional biomechanical quantities such as joint angles, ground reaction forces, joint torques, and muscle activation are not only fundamental indicators describing human movement status but also important quantitative representations for assessing musculoskeletal and nervous system-related diseases. To achieve unified assessment of these multidimensional indicators, computational biomechanical methods, such as musculoskeletal models, typically require high-quality joint kinematics and three-dimensional ground reaction forces as key inputs.

[0003] Currently, optical motion capture systems and force tables in laboratory environments can provide high-precision kinematic and three-dimensional ground reaction force measurements. However, these devices are typically fixed in the laboratory, resulting in complex testing procedures and limited space and scenarios, making it difficult to meet the needs of continuous, multi-scenario, and multi-dimensional evaluation of human movement in real-world environments. As sports health monitoring and behavioral research gradually move towards natural settings, there is an urgent need for an evaluation method that can output highly reliable multi-dimensional biomechanical indicators under non-laboratory conditions.

[0004] In recent years, motion assessment based on wearable sensors has received widespread attention. Inertial measurement units (IMUs) can be used to estimate kinematic information such as joint angles; foot measurement devices such as pressure insoles can provide dynamic information such as vertical ground reaction force and pressure center trajectory. However, IMUs cannot acquire dynamic information such as ground reaction force and joint torque during human movement, and pressure insoles can only provide ground reaction force in the vertical direction of human gait, unable to acquire three-dimensional ground reaction force during gait. On the other hand, due to the flexibility of materials and multi-array sensor structure, pressure insoles are prone to measurement bias during prolonged wear or multi-tasking activities. Therefore, existing motion assessment methods using wearable sensors are limited in scope and cannot achieve multi-dimensional biomechanical parameter assessment and comprehensive evaluation of joint angles, three-dimensional ground reaction force, joint torque, and muscle activation during human movement.

[0005] Existing research to address the aforementioned problems typically includes data-driven approaches, physical modeling methods, and multimodal fusion methods. However, most data-driven methods rely heavily on large amounts of task- or population-related training data and lack sufficient consideration for biomechanical consistency, resulting in limited generalization ability and interpretability. Some methods remain focused on single-index or specific scenario evaluations, making it difficult to form a reusable multidimensional evaluation framework. On the other hand, while physical modeling-based ground reaction force estimation methods offer some interpretability, they often struggle to simultaneously achieve accuracy and robustness under complex tasks and varying environmental conditions. Furthermore, many methods fail to explicitly handle the statically indeterminate phase of dual-support systems or fail to establish a unified, physically consistent solution framework with wearable kinematic and dynamic information. Therefore, current technology still lacks a method capable of integrating multimodal wearable sensing with cyber-physical model constraints under unconstrained conditions to achieve multidimensional musculoskeletal digital twin modeling. The shortcomings of existing technologies in addressing the need for multidimensional musculoskeletal digital twin modeling under unconstrained conditions are mainly reflected in the following aspects: (1) Lack of a means to obtain kinematics and three-dimensional ground reaction force synchronously and completely: Existing wearable solutions can often only obtain partial kinematic or dynamic information, which is difficult to support the unified evaluation of multiple indicators such as joint angle, three-dimensional ground reaction force, joint torque and muscle activation.

[0006] (2) Foot insole sensors have drift and scale deviation, which will be amplified in subsequent physical solutions: Insoles are prone to deviation during long-term wear and multi-task exercise, and the accuracy of vertical ground reaction force directly constrains the feasible domain and solution stability of physical solutions. There is a lack of cross-modal, online self-calibration mechanism.

[0007] (3) The load distribution in the dual-support stage cannot be uniquely determined: When both feet are in contact with the ground at the same time, the system is statically unstable. Single-mode sensing makes it difficult to determine the load ratio of the left and right feet, which leads to instability or inconsistency between the three-dimensional ground reaction force estimation and the subsequent dynamic solution.

[0008] (4) Lack of end-to-end multidimensional musculoskeletal digital twin modeling methods: Some existing methods rely on big data training and ignore physical consistency. Some physical model methods are difficult to adapt to complex tasks and environmental changes, and can only evaluate motion parameters under a single working condition or a single dimension. Overall, there is still a lack of generalizable, interpretable and physically consistent multidimensional musculoskeletal digital twin modeling methods. Summary of the Invention

[0009] To address the aforementioned issues, this invention aims to develop a wearable, sensor-driven, multidimensional musculoskeletal digital twin modeling method. This method enables stable reconstruction and evaluation of multidimensional human motion states under unconstrained conditions, thereby meeting the application needs of scenarios such as sports training, clinical rehabilitation, and mobile health.

[0010] The technical solution of the present invention: A wearable, sensor-driven, multidimensional musculoskeletal digital twin modeling method includes the following steps: Step 1: Calculate joint angles based on the inertial measurement unit; Step 2: Cross-modal calibration of pressure insole sensors; Step 3: Solving for 3D ground reaction forces based on physical constraints; Step 4: Calculate joint torque and muscle activation based on the musculoskeletal model.

[0011] Furthermore, step 1 specifically includes: Step 1 is used to calculate the joint angles between adjacent segments of the human body based on the posture information collected by multiple inertial measurement units throughout the body. The whole-body joint angle vector obtained in Step 1 is used as the input for calculating the human segment motion state in Step 2, and also as the kinematic input for solving the three-dimensional ground reaction force in Step 3 and calculating joint torque and muscle activation in Step 4.

[0012] set up: Represents the global coordinate system; Indicates the first An inertial measurement unit coordinate system; Indicates the first The individual human body segment coordinate system has a one-to-one correspondence between the inertial measurement unit and the human body segment; Indicates time; Represents quaternion multiplication; Indicates quaternion conjugation; This indicates the matrix transpose.

[0013] (1.1) Establish a fixed attitude relationship between the inertial measurement unit and human body segments; The subject maintained a neutral, stationary posture after wearing the inertial measurement unit. During the neutral, stationary posture period... Inside, obtain the first The attitude quaternions of each inertial measurement unit relative to the global coordinate system are calculated, and the average attitude quaternion is calculated:

[0014] in, Indicates time No. The attitude quaternion of an inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; Indicates the first The average attitude quaternion of each inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; Indicates quaternion averaging; superscript It indicates a neutral and static posture.

[0015] Average attitude quaternion Convert to average direction cosine matrix:

[0016] in, Indicates the first The average direction cosine matrix of each inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; This represents the operation of converting a quaternion into a direction cosine matrix.

[0017] At the same time, according to the first The anatomical coordinate system of an individual body segment in a neutral, static posture yields the first... The cosine matrix of the reference direction of an individual body segment relative to the global coordinate system .Depend on and Calculate the first The inertial measurement unit to the first Fixed rotation matrix of individual body segment coordinate system :

[0018] Among them, superscript Indicates a reference state.

[0019] Fixed rotation matrix Convert to a fixed rotation quaternion:

[0020] in, Indicates the first The inertial measurement unit to the first Fixed rotational quaternions for an individual anthropometric coordinate system; This represents the operation of converting a rotation matrix into a quaternion. The purpose of the above process is to determine the fixed attitude deviation between the coordinate system of the inertial measurement unit (IMU) and the coordinate system of the human segment after the IMU is attached to the human segment.

[0021] (1.2) Calculate the dynamic global attitude of the inertial measurement unit; During the dynamic movement of the human body, the first... The angular velocity vector acquired by the first inertial measurement unit in its own coordinate system is used to calculate the angular velocity vector of the second inertial measurement unit. Rate of change of dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system Then, the dynamic attitude quaternion for the next time step is obtained through integration:

[0022] in, equal The first derivative with respect to time; Indicates time No. The angular velocity vector collected by each inertial measurement unit in its own coordinate system; Indicates time t No. The dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system; Indicates time No. The dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system; This represents the quaternion normalization operation, used to ensure that the updated dynamic attitude quaternions are still unit quaternions. This represents the change in time. The purpose of this step is to calculate the real-time attitude of the inertial measurement unit (IMU) during its motion, based on the angular velocity measured by the IMU.

[0023] (1.3) Obtain the human segment posture from the posture of the inertial measurement unit; Based on the fixed rotation quaternion obtained in step (1.1) Update the obtained step (1.2) Dynamic attitude quaternions of each inertial measurement unit Convert to the first Dynamic posture quaternions of individual body segments :

[0024] in, Indicates time No. The attitude quaternion of the individual body segment coordinate system relative to the global coordinate system; Represents a fixed rotation quaternion The conjugate quaternion, due to For a unit quaternion, its conjugate quaternion is... This is equivalent to its inverse quaternion. The purpose of this step is to convert the real-time attitude of the inertial measurement unit into the real-time attitude of the human body segments.

[0025] (1.4) Calculate the relative posture between adjacent human body segments; For any joint, let the proximal human segment of that joint be... The distal human body segment is Based on the quaternion of proximal human segment posture. and distal human segment posture quaternions Calculate the relative attitude quaternion of the joint. :

[0026] in, Indicates time The quaternion of joint relative posture of distal human segments relative to proximal human segments; The conjugate quaternion is used to transform the attitude relationships in the global coordinate system to the proximal human segment coordinate system. The purpose of this step is to represent joint motion using the attitude difference between two adjacent human segments.

[0027] (1.5) Convert the relative joint posture quaternions into joint angles; Let the joint relative pose quaternion obtained in step (1.4) be... for:

[0028] in, Indicates time Scalar components of the joint relative posture quaternion; , , Representing time respectively The three vector components of the joint relative posture quaternion.

[0029] according to Rotation sequence, joint relative pose quaternion Converted to joint internal / external rotation angles, joint adduction / abduction angles, and joint flexion / extension angles: (10) (11) (12) in, Indicates time The angles of internal and external rotation of the joint; Indicates time The angle of adduction and abduction of the joint; Indicates time The angle of joint flexion and extension. Represents the arctangent function. This represents the arcsine function. The purpose of this step is to convert the relative joint pose in quaternion form into joint angles in three directions.

[0030] Equations (9) to (12) are applied to all joints of the human body to be calculated. Let the human body have a total of... One joint to be calculated Indicates the joint number. For the first The three angles obtained from equations (10) to (12) for each joint are denoted as follows: , and .

[0031] (1.6) Output the angle vector of all joints in the body; The first The three angles of each joint form a joint angle subvector. Then, all joint angle subvectors are combined according to a preset joint order to form a whole-body joint angle vector. :

[0032] in, Indicates time No. Joint angle subvectors of each joint; , and Let each represent the result obtained in step (1.5). Each joint at any time The internal and external rotation angles, adduction and abduction angles, and flexion and extension angles; Indicates the number of joints to be calculated; superscript This represents the vector transpose. Step 1 outputs the whole-body joint angle vector. Used in step 2 to calculate the human segmental motion state and total vertical inertial reaction force, used in step 3 to solve the three-dimensional ground reaction force, and used in step 4 to calculate joint torque and muscle activation.

[0033] Furthermore, step 2 specifically includes: Step 2 is based on the whole-body joint angle vector obtained in Step 1. The musculoskeletal model is driven to calculate the vertical acceleration of the center of mass of the human segment, and the total vertical inertial reaction force of the human body is obtained according to the Newton-Euler method. Then, using the total vertical inertial reaction force of the human body as the physical reference, the vertical ground reaction forces of the left and right feet output by the pressure insole sensor are calibrated. The calibrated vertical ground reaction forces of the left and right feet and the vertical load ratio of the left and right feet obtained in step 2 are used in step 3.

[0034] (2.1) Calculate the total vertical inertial reaction force of the human body; Based on the whole-body joint angle vector obtained in step 1 , driving the musculoskeletal model to calculate the first Vertical acceleration of the center of mass of an individual body segment The mass of human body segments Vertical acceleration and gravitational acceleration constant Calculate the total vertical inertial reaction force of the human body. :

[0035] in, Indicates the total number of human body segments; subscript Indicates the vertical direction of the global coordinate system.

[0036] (2.2) Establish a vertical ground reaction force calibration model; The pressure insole sensor synchronously outputs the original vertical ground reaction force of the left foot. and the original vertical ground reaction force of the right foot A scale factor is introduced for each foot. , and zero drift item , The calibrated vertical ground reaction force of the left foot was obtained. and the reaction force perpendicular to the ground of the right foot :

[0037] Left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot That is, the three-dimensional ground reaction force of the left foot in the global coordinate system. The components in the direction and the three-dimensional ground reaction force of the right foot in the global coordinate system Components in direction; (2.3) Solve for the calibration parameters within the sliding window; Let time be The center sliding window is The time frame within the sliding window is The calibration parameters to be solved are denoted as a decision variable vector. And the change in scale factor between adjacent windows is denoted as and :

[0038] The total vertical inertial reaction force of the human body obtained in step (2.1) As a physical benchmark, establish the cross-modal calibration optimization function for the pressure insole sensor. : (17) in, Represents time frame The residual weights; Represents the regularization weights for zero-drift terms; This represents the scaling factor smoothing regularization weight. Equation (17) outputs the scaling factor. , and zero drift item , Used in equation (15), to make the sum of the vertical ground reaction forces of the calibrated left and right feet as close as possible to the total vertical inertial reaction force of the human body, while suppressing excessive zero drift and abrupt changes in scale factor.

[0039] Equation (17) satisfies the following constraints:

[0040] in, Represents time frame The left foot supporting state variable, This indicates that the left foot is in the support phase. This indicates that the left foot is in the swing phase; Represents time frame The right foot support state variable, This indicates that the right foot is in the support phase. This indicates that the right foot is in the swing phase; This represents the minimum permissible value of the scaling factor; This indicates the maximum permissible value of the scale factor; This represents the maximum permissible value of the absolute value of the zero-drift term.

[0041] (2.4) Output the vertical load ratio of the left and right feet; Based on the calibrated left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot Calculate the vertical load ratio of the left foot and the ratio of vertical load on the right foot :

[0042] Step 2 Output , , and This is used for solving the three-dimensional ground reaction force in step 3.

[0043] Furthermore, step 3 specifically includes: Step 3 is based on the whole-body joint angle vector output in Step 1. Step 2 output: calibrated left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot Left foot vertical load ratio and the ratio of vertical load on the right foot The three-dimensional ground reaction forces of the left and right feet are calculated by combining the pressure center trajectory output by the pressure insole sensor. The three-dimensional ground reaction forces of the left and right feet output in step 3 are used for the inverse dynamics calculation in step 4.

[0044] (3.1) Calculate the total inertial resultant force vector of the human body; Based on the whole-body joint angle vector obtained in step 1 , driving the musculoskeletal model to calculate the first Three-dimensional acceleration vector of the centroid of an individual segment The mass of human body segments and three-dimensional acceleration vector Calculate the total inertial resultant force vector of the human body. :

[0045] in , and These represent the total inertial force of the human body in the global coordinate system. , , Components in direction.

[0046] (3.2) Calculate the total inertial moment vector of the human body; Based on the centroid position vector of human segments Linear momentum vector Inertial tensor and angular velocity vector Calculate the total inertial moment vector of the human body. :

[0047] in, express The first derivative with respect to time; Indicates time No. The inertial tensor of an individual body segment about its center of mass; Indicates time No. The angular velocity vector of an individual body segment; express The first derivative with respect to time; This represents the cross product of vectors.

[0048] The vector of total inertial moment of the human body Decomposed into three directional components:

[0049] in, , and Representing time respectively Total inertial torque of the human body in the global coordinate system , , Components in direction.

[0050] (3.3) Calculate the torque arm components of the left and right feet; Pressure insole sensor output Foot pressure center coordinates The musculoskeletal model is based on the whole-body joint angle vectors from step 1. Calculate the coordinates of the human body's center of mass .Depend on and Calculate the first Components of the left and right foot torque arms from the foot pressure center to the body's center of mass :

[0051] in, Indicates foot separation, Indicates the left foot. Indicates the right foot; It represents the directions of the three coordinate axes of the global coordinate system.

[0052] (3.4) Calculate the resultant force and resultant moment of the three-dimensional ground reaction force during human movement; Let the three-dimensional ground reaction force vector of the left foot be... The three-dimensional ground reaction vector of the right foot is :

[0053] in, Indicates time Three-dimensional ground reaction force vector of the left foot; Indicates time Three-dimensional ground reaction force vector of the right foot; , , These represent the three-dimensional ground reaction forces of the left foot in the global coordinate system. , , Components in direction; , , These represent the three-dimensional ground reaction forces of the right foot in the global coordinate system. , , Components in direction.

[0054] The three-dimensional ground reaction forces of the left and right feet satisfy the resultant force equilibrium relationship:

[0055] The resultant moment components are calculated from the three-dimensional ground reaction forces of the left and right feet and the moment arms of the left and right feet. , and : (26) (3.5) Initialize the horizontal ground reaction forces of the left and right feet; During the dual-support phase, based on the left foot vertical load ratio obtained in step 2... and the ratio of vertical load on the right foot Initialize the horizontal ground reaction forces for the left and right feet:

[0056] in, , , and These represent the initial values ​​of the horizontal ground reaction forces of the left and right feet, respectively; superscript This indicates the initial values ​​for optimization. The purpose of this step is to provide an initial estimate of the distribution of horizontal ground reaction forces using the vertical load ratio.

[0057] (3.6) Establish friction cone constraints; To ensure the physical feasibility of foot-to-sole contact, a friction cone constraint is established for the three-dimensional ground reaction forces of the left and right feet:

[0058] in, This represents the coefficient of friction between the sole of the shoe and the ground. Equation (28) ensures that the horizontal ground reaction force obtained from the solution does not exceed the friction range that the sole of the foot can withstand.

[0059] (3.7) Solve for the three-dimensional ground reaction forces of the left and right feet; Using the horizontal ground reaction forces of the left and right feet as the unknowns, a three-dimensional decision variable vector for solving the ground reaction force is established. :

[0060] By minimizing the calculated total torque With Newton-Euler total moment of inertia The differences between them are used to establish an optimization function for solving the three-dimensional ground reaction force. :

[0061] Among them, the first two equality constraints in equation (30) ensure that the sum of the horizontal ground reaction forces of the left and right feet is consistent with the total inertial resultant force of the human body, and the latter two inequality constraints ensure that the three-dimensional ground reaction forces of the left and right feet satisfy the friction cone constraint.

[0062] By solving equation (30), the three-dimensional ground reaction vector of the left foot is obtained. and the three-dimensional ground reaction vector of the right foot :

[0063] in, and Together, they serve as the external load input for the inverse dynamics calculation in step 4.

[0064] Furthermore, step 4 specifically includes: Step 4 is based on the whole-body joint angle vector output in Step 1. and the three-dimensional ground reaction force of the left foot output in step 3 Three-dimensional ground reaction force of the right foot The joint torques are calculated using a musculoskeletal model, and muscle activation is solved through static optimization under the constraint of joint torque balance. The joint torques and muscle activation output in step 4, together with the outputs of steps 1 and 3, constitute the multidimensional biomechanical evaluation results.

[0065] (4.1) Calculate the joint moment vector using inverse dynamics; Vector of joint angles throughout the body Three-dimensional ground reaction force of the left foot and the three-dimensional ground reaction force of the right foot Input a musculoskeletal model. Under the condition of three-dimensional ground reaction force as an external load, calculate the joint moment vector according to the dynamic equations of the musculoskeletal model. :

[0066] in, This represents the first derivative of the vector of angles of all joints with respect to time. This represents the second derivative of the vector of angles of all joints with respect to time. Represents the musculoskeletal model in terms of the joint angle vectors throughout the body. The generalized mass matrix below; The Coriolis force and centrifugal force correlation matrix representing the musculoskeletal model; This represents the gravity term vector of the musculoskeletal model.

[0067] (4.2) Establish a balance between muscle force and joint torque; Let the first A muscle at all times Muscle activation level is The maximum isometric strength of this muscle is This muscle is relative to the first The lever arm of each joint degree of freedom is . No. One muscle against the first The torque contribution of each joint degree of freedom is The sum of the torque contributions of all relevant muscles equals the joint torque of that degree of freedom. :

[0068] in, This indicates the number of muscles involved in balancing joint torque.

[0069] (4.3) Solve for muscle activation through static optimization; Combining all muscle activation values ​​into a muscle activation vector and with neuromuscular control cost function With the goal of minimizing, solve for the muscle activation level at each time step:

[0070] in, This indicates the range of values ​​for muscle activation. Equation (34) aims to find the combination of muscle activations that minimizes the overall cost of muscle activation while satisfying joint torque balance.

[0071] (4.4) Output multidimensional biomechanical evaluation results; Joint torque vector and muscle activation vector Combined with the whole-body joint angle vector output in step 1 And the three-dimensional ground reaction vector of the left foot output in step 3. and the three-dimensional ground reaction vector of the right foot The multidimensional biomechanical evaluation results of the musculoskeletal digital twin system were obtained. :

[0072] Through steps 1 to 4 above, the musculoskeletal digital twin system forms a continuous data processing flow: the whole-body joint angle vector obtained in step 1 is used in steps 2, 3 and 4; the calibrated left and right foot vertical ground reaction force and left and right foot vertical load ratio obtained in step 2 are used in step 3; the left and right foot three-dimensional ground reaction force obtained in step 3 is used in step 4; step 4 finally outputs joint torque and muscle activation, which together with the whole-body joint angle vector and the left and right foot three-dimensional ground reaction force constitute a multi-dimensional biomechanical evaluation result.

[0073] The beneficial effects of this invention are: This invention enables the simultaneous acquisition of key inputs to whole-body kinematics and plantar dynamics without laboratory equipment. Through a wireless somatosensory network comprised of multiple IMUs and pressure-sensitive insoles for both feet, it achieves the simultaneous acquisition of whole-body segmental kinematic information and plantar dynamics-related information. This allows for the acquisition of precise kinematic and dynamic data of human movement even in outdoor settings without optical motion capture or force tables, overcoming the limitations of fixed laboratory environments.

[0074] This invention enables cross-modal self-calibration to ensure consistent and reliable scalability of vertical force inputs from the left and right feet. Using the inertial vertical resultant force derived from a whole-body IMU as the calibration benchmark, the vertical ground reaction forces output by the insole from the left and right feet are dynamically calibrated, resolving deviations caused by long-term insole wear and outputting stable and reliable ground reaction forces and weight relationships for the left and right feet. These weights are further used as initialization information for horizontal force distribution during the dual-support stage and subsequent constraint optimization, ensuring scalability and stability in subsequent calculations.

[0075] This invention enables the solution of complete three-dimensional ground reaction forces for both feet under information fusion and physical constraints. Addressing the issues of static indeterminate distribution of three-dimensional ground reaction forces during the dual-support stage and the difficulty in uniquely determining these forces using single-modal sensing, this invention embeds the multi-source observation consistency objective with Newton-Euler dynamics equations and friction cone feasible region constraints into a unified optimization framework. While ensuring physical feasibility, it minimizes the inconsistency in total force and total torque among multi-source data, thereby analytically obtaining the complete three-dimensional ground reaction force vectors for both feet and significantly improving the robustness, interpretability, and cross-task generalization ability of the solution.

[0076] This invention enables the formation of an end-to-end multidimensional biomechanical assessment closed loop encompassing kinematics, dynamics, and muscle activation. By using whole-body joint angles and three-dimensional ground reaction forces as unified inputs, and combining inverse dynamics to solve for joint net torques, and then statically optimizing muscle activation under torque balance constraints, a multidimensional assessment output is generated, including joint kinematics, three-dimensional ground reaction forces, joint torques, and muscle activation. Furthermore, this closed-loop framework can be applied to various movement tasks such as walking, running, and climbing stairs, providing interpretable and scalable unconstrained multidimensional assessment capabilities for sports training, clinical rehabilitation, and mobile health monitoring. Attached Figure Description

[0077] Figure 1 This is a schematic diagram of a wearable, sensor-driven, multidimensional musculoskeletal digital twin modeling method. Figure 2 This is a flowchart illustrating the overall technology of a wearable, sensor-driven, multidimensional musculoskeletal digital twin modeling method. Figure 3 This is a flowchart of a multidimensional musculoskeletal digital twin modeling technology based on wearable sensing and the fusion of kinematic and dynamic data. Detailed Implementation

[0078] The technical solution of the present invention will be further described below through specific embodiments.

[0079] A wearable sensor-driven multidimensional musculoskeletal digital twin modeling method, the overall schematic diagram of which is shown below. Figure 1 It includes the following steps ( Figure 2 ): Step 1: Calculate joint angles based on the inertial measurement unit; Step 2: Cross-modal calibration of pressure insole sensors; Step 3: Solving for 3D ground reaction forces based on physical constraints; Step 4: Calculate joint torque and muscle activation based on the musculoskeletal model.

[0080] Step 1 specifically includes: Step 1 calculates the joint angles between adjacent segments of the human body based on attitude information collected by multiple inertial measurement units (IMUs) throughout the body, including data from accelerometers, gyroscopes, and magnetometers. A total of 17 IMUs are deployed in the pelvis, trunk, both thighs, both calves, both feet, and relevant segments of the upper limbs to obtain the motion posture of the main segments of the body. During data preprocessing, outlier removal, coordinate orientation unification, low-pass filtering, and smoothing are performed on the acceleration, angular velocity, and attitude data output by the IMUs. The whole-body joint angle vectors obtained in Step 1 serve as inputs for calculating the segmental motion state in Step 2, and also as kinematic inputs for solving the three-dimensional ground reaction force in Step 3 and calculating joint torques and muscle activation in Step 4.

[0081] set up: Represents the global coordinate system; Indicates the first An inertial measurement unit coordinate system; Indicates the first The individual human body segment coordinate system has a one-to-one correspondence between the inertial measurement unit and the human body segment; Indicates time; Represents quaternion multiplication; Indicates quaternion conjugation; This indicates the matrix transpose.

[0082] (1.1) Establish a fixed attitude relationship between the inertial measurement unit and human body segments; After wearing the inertial measurement unit, the subject maintained a neutral, stationary posture for at least 5 seconds. During this period of neutral, stationary posture... Inside, obtain the first The attitude quaternions of each inertial measurement unit relative to the global coordinate system are calculated, and the average attitude quaternion is calculated:

[0083] in, Indicates time No. The attitude quaternion of an inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; Indicates the first The average attitude quaternion of each inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; Indicates quaternion averaging; superscript It indicates a neutral and static posture.

[0084] Average attitude quaternion Convert to average direction cosine matrix:

[0085] in, Indicates the first The average direction cosine matrix of each inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; This represents the operation of converting a quaternion into a direction cosine matrix.

[0086] At the same time, according to the first The anatomical coordinate system of an individual body segment in a neutral, static posture yields the first... The cosine matrix of the reference direction of an individual body segment relative to the global coordinate system .Depend on and Calculate the first The inertial measurement unit to the first Fixed rotation matrix of individual body segment coordinate system :

[0087] Among them, superscript Indicates a reference state.

[0088] Fixed rotation matrix Convert to a fixed rotation quaternion:

[0089] in, Indicates the first The inertial measurement unit to the first Fixed rotational quaternions for an individual anthropometric coordinate system; This represents the operation of converting a rotation matrix into a quaternion. The purpose of the above process is to determine the fixed attitude deviation between the coordinate system of the inertial measurement unit (IMU) and the coordinate system of the human segment after the IMU is attached to the human segment.

[0090] (1.2) Calculate the dynamic global attitude of the inertial measurement unit; During the dynamic movement of the human body, the first... The angular velocity vector acquired by the first inertial measurement unit in its own coordinate system is used to calculate the angular velocity vector of the second inertial measurement unit. Rate of change of dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system Then, the dynamic attitude quaternion for the next time step is obtained through integration:

[0091] in, equal The first derivative with respect to time; Indicates time No. The angular velocity vector collected by each inertial measurement unit in its own coordinate system; Indicates time t No. The dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system; Indicates time No. The dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system; This represents the quaternion normalization operation, used to ensure that the updated dynamic attitude quaternions are still unit quaternions. This represents the change in time. The purpose of this step is to calculate the real-time attitude of the inertial measurement unit (IMU) during its motion, based on the angular velocity measured by the IMU.

[0092] (1.3) Obtain the human segment posture from the posture of the inertial measurement unit; Based on the fixed rotation quaternion obtained in step (1.1) Update the obtained step (1.2) Dynamic attitude quaternions of each inertial measurement unit Convert to the first Dynamic posture quaternions of individual body segments :

[0093] in, Indicates time No. The attitude quaternion of the individual body segment coordinate system relative to the global coordinate system; Represents a fixed rotation quaternion The conjugate quaternion, due to For a unit quaternion, its conjugate quaternion is... This is equivalent to its inverse quaternion. The purpose of this step is to convert the real-time attitude of the inertial measurement unit into the real-time attitude of the human body segments.

[0094] (1.4) Calculate the relative posture between adjacent human body segments; For any joint, let the proximal human segment of that joint be... The distal human body segment is Based on the quaternion of proximal human segment posture. and distal human segment posture quaternions Calculate the relative attitude quaternion of the joint. :

[0095] in, Indicates time The quaternion of joint relative posture of distal human segments relative to proximal human segments; The conjugate quaternion is used to transform the attitude relationships in the global coordinate system to the proximal human segment coordinate system. The purpose of this step is to represent joint motion using the attitude difference between two adjacent human segments.

[0096] (1.5) Convert the relative joint posture quaternions into joint angles; Let the joint relative pose quaternion obtained in step (1.4) be... for:

[0097] in, Indicates time Scalar components of the joint relative posture quaternion; , , Representing time respectively The three vector components of the joint relative posture quaternion.

[0098] according to Rotation sequence, joint relative pose quaternion Converted to joint internal / external rotation angles, joint adduction / abduction angles, and joint flexion / extension angles: (10) (11) (12) in, Indicates time The angles of internal and external rotation of the joint; Indicates time The angle of adduction and abduction of the joint; Indicates time The angle of joint flexion and extension. Represents the arctangent function. This represents the arcsine function. The purpose of this step is to convert the relative joint pose in quaternion form into joint angles in three directions.

[0099] Equations (9) to (12) are applied to all joints of the human body to be calculated. Let the human body have a total of... One joint to be calculated Indicates the joint number. For the first The three angles obtained from equations (10) to (12) for each joint are denoted as follows: , and .

[0100] (1.6) Output the angle vector of all joints in the body; The first The three angles of each joint form a joint angle subvector. Then, all joint angle subvectors are combined according to a preset joint order to form a whole-body joint angle vector. :

[0101] in, Indicates time No. Joint angle subvectors of each joint; , and Let each represent the result obtained in step (1.5). Each joint at any time The internal and external rotation angles, adduction and abduction angles, and flexion and extension angles; Indicates the number of joints to be calculated; superscript This represents the vector transpose. Step 1 outputs the whole-body joint angle vector. Used in step 2 to calculate the human segmental motion state and total vertical inertial reaction force, used in step 3 to solve the three-dimensional ground reaction force, and used in step 4 to calculate joint torque and muscle activation.

[0102] Step 2 specifically refers to: Step 2 uses the total vertical inertial reaction force of the human body calculated by the inertial measurement unit and the musculoskeletal model as a physical reference to perform cross-modal calibration on the vertical load of the left and right feet output by the pressure insole sensor. This is based on the whole-body joint angle vector obtained in Step 1. The musculoskeletal model is driven to calculate the vertical acceleration of the center of mass of the human segment, and the total vertical inertial reaction force of the human body is obtained according to the Newton-Euler method. Then, using the total vertical inertial reaction force of the human body as the physical reference, the vertical ground reaction forces of the left and right feet output by the pressure insole sensor are calibrated. The calibrated vertical ground reaction forces of the left and right feet and the vertical load ratio of the left and right feet obtained in step 2 are used in step 3.

[0103] (2.1) Calculate the total vertical inertial reaction force of the human body; The total vertical inertial reaction force of the human body is determined by the mass of each segment and the vertical acceleration of the segment's center of mass. This is based on the whole-body joint angle vector obtained in step 1. , driving the musculoskeletal model to calculate the first Vertical acceleration of the center of mass of an individual body segment The mass of human body segments Vertical acceleration and gravitational acceleration constant Calculate the total vertical inertial reaction force of the human body. :

[0104] in, Indicates the total number of human body segments; subscript Indicates the vertical direction of the global coordinate system.

[0105] (2.2) Establish a vertical ground reaction force calibration model; The original vertical ground reaction forces of the left and right feet are corrected using a scaling factor and a zero-drift term. The scaling factor compensates for the overall underestimation or overestimation of the load amplitude by the pressure insole sensor, while the zero-drift term compensates for the bias output of the pressure insole sensor that still exists under no-load conditions. The pressure insole sensor synchronously outputs the original vertical ground reaction force of the left foot. and the original vertical ground reaction force of the right foot A scale factor is introduced for each foot. , and zero drift item , The calibrated vertical ground reaction force of the left foot was obtained. and the reaction force perpendicular to the ground of the right foot :

[0106] Left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot That is, the three-dimensional ground reaction force of the left foot in the global coordinate system. The components in the direction and the three-dimensional ground reaction force of the right foot in the global coordinate system Components in direction; (2.3) Solve for the calibration parameters within the sliding window; Using a sliding window to solve for calibration parameters allows the calibration process to simultaneously possess local adaptability and temporal continuity. The sliding window covers multiple sampling frames within a certain time range, preventing single-frame noise from having an excessive impact on the calibration parameters. Let time be... The center sliding window is The time frame within the sliding window is The calibration parameters to be solved are denoted as a decision variable vector. And the change in scale factor between adjacent windows is denoted as and :

[0107] The total vertical inertial reaction force of the human body obtained in step (2.1) As a physical benchmark, establish the cross-modal calibration optimization function for the pressure insole sensor. : (17) in, Represents time frame The residual weights of each time frame within the sliding window can be determined based on the foot contact state and signal stability. Represents the regularization weights for zero-drift terms; This represents the scaling factor smoothing regularization weight. Equation (17) outputs the scaling factor. , and zero drift item , Used in equation (15), to make the sum of the vertical ground reaction forces of the calibrated left and right feet as close as possible to the total vertical inertial reaction force of the human body, while suppressing excessive zero drift and abrupt changes in scale factor.

[0108] To avoid the pressure insole sensor outputting a false load when the foot leaves the ground, and to prevent the generation of a negative vertical ground reaction force that does not conform to physical meaning after calibration, it is necessary to supplement the optimization function with constraints, and to add upper and lower limits for the scale factor and the zero drift term to limit the variation range of the calibration parameters, so that the calibration results are kept within the reasonable physical range of the pressure insole sensor. Equation (17) needs to satisfy the following constraints:

[0109] in, Represents time frame The left foot supporting state variable, This indicates that the left foot is in the support phase. This indicates that the left foot is in the swing phase; Represents time frame The right foot support state variable, This indicates that the right foot is in the support phase. This indicates that the right foot is in the swing phase; This represents the minimum permissible value of the scaling factor; This indicates the maximum permissible value of the scale factor; This represents the maximum permissible value of the absolute value of the zero-drift term.

[0110] (2.4) Output the vertical load ratio of the left and right feet; The vertical load ratio between the left and right feet reflects the distribution of body weight and dynamic load between the two feet at any given moment, and also reflects the weight transfer process between the forefoot and hindfoot. This is based on the calibrated vertical ground reaction force of the left foot. and the reaction force perpendicular to the ground of the right foot Calculate the vertical load ratio of the left foot and the ratio of vertical load on the right foot :

[0111] Step 2 Output , , and This is used for solving the three-dimensional ground reaction force in step 3.

[0112] Step 3 specifically includes: The purpose of solving the three-dimensional ground reaction force in step 3 is to determine the three-dimensional ground reaction force of the left and right feet in the global coordinate system, without a fixed force table or relying on a laboratory force measurement platform, based on the overall human dynamics, the vertical load of the pressure insole sensor, and the trajectory of the pressure center. This is based on the whole-body joint angle vector output in step 1. Step 2 output: calibrated left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot Left foot vertical load ratio and the ratio of vertical load on the right foot The three-dimensional ground reaction forces of the left and right feet are calculated by combining the pressure center trajectory output by the pressure insole sensor. The three-dimensional ground reaction forces of the left and right feet output in step 3 are used for the inverse dynamics calculation in step 4.

[0113] (3.1) Calculate the total inertial resultant force vector of the human body; The total inertial force of the human body reflects the net external force required for the linear acceleration of each segment, while the total inertial torque reflects the overall torque requirement generated by the rotational inertia, angular velocity, and angular acceleration of each segment. Both provide overall dynamic constraints for solving the three-dimensional ground reaction forces of the left and right feet. Based on the whole-body joint angle vector obtained in step 1... , driving the musculoskeletal model to calculate the first Three-dimensional acceleration vector of the centroid of an individual segment The mass of human body segments and three-dimensional acceleration vector Calculate the total inertial force vector of the human body according to the Newton-Euler equations. :

[0114] in , and These represent the total inertial force of the human body in the global coordinate system. , , Components in direction.

[0115] (3.2) Calculate the total inertial moment vector of the human body; Based on the centroid position vector of human segments Linear momentum vector Inertial tensor and angular velocity vector Calculate the total inertial moment vector of the human body. :

[0116] in, express The first derivative with respect to time; Indicates time No. The inertial tensor of an individual body segment about its center of mass; Indicates time No. The angular velocity vector of an individual body segment; express The first derivative with respect to time; This represents the cross product of vectors.

[0117] The vector of total inertial moment of the human body Decomposed into three directional components:

[0118] in, , and Representing time respectively Total inertial torque of the human body in the global coordinate system , , Components in direction.

[0119] (3.3) Calculate the torque arm components of the left and right feet; Pressure insole sensor output Foot pressure center coordinates The musculoskeletal model is based on the whole-body joint angle vectors from step 1. Calculate the coordinates of the human body's center of mass .Depend on and Calculate the first Components of the left and right foot torque arms from the foot pressure center to the body's center of mass :

[0120] in, Indicates foot separation, Indicates the left foot. Indicates the right foot; It represents the directions of the three coordinate axes of the global coordinate system.

[0121] (3.4) Calculate the resultant force and resultant moment of the three-dimensional ground reaction force during human movement; Let the three-dimensional ground reaction force vector of the left foot be... The three-dimensional ground reaction vector of the right foot is :

[0122] in, Indicates time Three-dimensional ground reaction force vector of the left foot; Indicates time Three-dimensional ground reaction force vector of the right foot; , , These represent the three-dimensional ground reaction forces of the left foot in the global coordinate system. , , Components in direction; , , These represent the three-dimensional ground reaction forces of the right foot in the global coordinate system. , , Components in direction.

[0123] The three-dimensional ground reaction forces of the left and right feet satisfy the resultant force equilibrium relationship:

[0124] The resultant moment components are calculated from the three-dimensional ground reaction forces of the left and right feet and the moment arms of the left and right feet. , and : (26) The resultant force of the three-dimensional ground reaction forces of the left and right feet should be consistent with the resultant force of the total inertia of the human body, and the resultant torque generated by the three-dimensional ground reaction forces of the left and right feet about the reference point should be consistent with the total inertial torque of the human body. This equilibrium relationship reflects the dynamic correspondence between the external ground reaction force and the segmental inertial motion of the human body during human movement.

[0125] (3.5) Initialize the horizontal ground reaction forces of the left and right feet; During the dual-support phase, based on the left foot vertical load ratio obtained in step 2... and the ratio of vertical load on the right foot Initialize the horizontal ground reaction forces for the left and right feet:

[0126] in, , , and These represent the initial values ​​of the horizontal ground reaction forces of the left and right feet, respectively; superscript This represents the initial values ​​for optimization. The purpose of this step is to provide initial values ​​for the distribution of horizontal ground reaction forces using the vertical load ratio as an initial estimate for the optimization solution, thereby improving solution stability and convergence speed.

[0127] (3.6) Establish friction cone constraints; To ensure the physical feasibility of foot-to-sole contact, a friction cone constraint is established for the three-dimensional ground reaction forces of the left and right feet:

[0128] in, The coefficient of friction between the sole of the shoe and the ground is represented by . The purpose of equation (28) is to ensure that the horizontal ground reaction force obtained by the solution does not exceed the friction range that the sole of the foot can withstand, so as to avoid the situation where the horizontal ground reaction force is too large and cannot be achieved in the real physical environment.

[0129] (3.7) Solve for the three-dimensional ground reaction forces of the left and right feet; Using the horizontal ground reaction forces of the left and right feet as the unknowns, a three-dimensional decision variable vector for solving the ground reaction force is established. :

[0130] By minimizing the calculated total torque With Newton-Euler total moment of inertia The differences between them are used to establish an optimization function for solving the three-dimensional ground reaction force. :

[0131] Among them, the first two equality constraints in equation (30) ensure that the sum of the horizontal ground reaction forces of the left and right feet is consistent with the total inertial resultant force of the human body, and the latter two inequality constraints ensure that the three-dimensional ground reaction forces of the left and right feet satisfy the friction cone constraint.

[0132] By solving equation (30), the three-dimensional ground reaction vector of the left foot is obtained. and the three-dimensional ground reaction vector of the right foot :

[0133] in, and Together, they serve as the external load input for the inverse dynamics calculation in step 4.

[0134] Step 4 specifically includes: Step 4 is based on the whole-body joint angle vector output in Step 1. and the three-dimensional ground reaction force of the left foot output in step 3 Three-dimensional ground reaction force of the right foot Joint torques are calculated using a musculoskeletal model, and muscle activation is solved through static optimization under joint torque balance constraints. Figure 3 The joint torque and muscle activation output in step 4, together with the outputs of steps 1 and 3, constitute the multidimensional biomechanical evaluation results.

[0135] (4.1) Calculate the joint moment vector using inverse dynamics; Vector of joint angles throughout the body Three-dimensional ground reaction force of the left foot and the three-dimensional ground reaction force of the right foot Input a musculoskeletal model. Under the condition of three-dimensional ground reaction force as an external load, calculate the joint moment vector according to the dynamic equations of the musculoskeletal model. :

[0136] in, This represents the first derivative of the vector of angles of all joints with respect to time. This represents the second derivative of the vector of angles of all joints with respect to time. Represents the musculoskeletal model in terms of the joint angle vectors throughout the body. The generalized mass matrix below; The Coriolis force and centrifugal force correlation matrix representing the musculoskeletal model; The gravity term vector representing the musculoskeletal model is used to reduce the amplification of noise during the differentiation process. The joint angle sequence is first filtered and smoothed before calculating its first and second derivatives. The joint torques obtained through inverse dynamics calculations reflect the net joint torques required to maintain movement under the current motion state and external load conditions.

[0137] (4.2) Establish a balance between muscle force and joint torque; The balance between muscle force and joint torque describes the driving effect of muscles on joint movement. Different muscles contribute different torques to different degrees of freedom of the joint due to differences in their origin and insertion points, the number of joints they cross, and their lever arms. Let the ______ be the ______. A muscle at all times Muscle activation level is The maximum isometric strength of this muscle is This muscle is relative to the first The lever arm of each joint degree of freedom is . No. One muscle against the first The torque contribution of each joint degree of freedom is The sum of the torque contributions of all relevant muscles equals the joint torque of that degree of freedom. :

[0138] in, This indicates the number of muscles involved in balancing joint torque.

[0139] (4.3) Solve for muscle activation through static optimization; Combining all muscle activation values ​​into a muscle activation vector and with neuromuscular control cost function With the goal of minimizing, solve for the muscle activation level at each time step:

[0140] in, This indicates the range of values ​​for muscle activation. Equation (34) aims to find the combination of muscle activations that minimizes the overall cost of muscle activation while satisfying joint torque balance.

[0141] (4.4) Output multidimensional biomechanical evaluation results; Joint torque vector and muscle activation vector Then, combined with the whole-body joint angle vector output in step 1 And the three-dimensional ground reaction vector of the left foot output in step 3. and the three-dimensional ground reaction vector of the right foot The multidimensional biomechanical evaluation results of the musculoskeletal digital twin system were obtained. :

[0142] Through steps 1 to 4 above, the musculoskeletal digital twin system forms a continuous data processing flow: the whole-body joint angle vector obtained in step 1 is used in steps 2, 3 and 4; the calibrated left and right foot vertical ground reaction force and left and right foot vertical load ratio obtained in step 2 are used in step 3; the left and right foot three-dimensional ground reaction force obtained in step 3 is used in step 4; step 4 finally outputs joint torque and muscle activation, which together with the whole-body joint angle vector and the left and right foot three-dimensional ground reaction force constitute a multi-dimensional biomechanical evaluation result.

[0143] Compared with the prior art, the advantages of the present invention are: 1. This invention is based on a full-body sensory network composed of an IMU and a pressure insole. It can simultaneously acquire key inputs of whole-body kinematics and plantar dynamics without optical motion capture or force measurement platform. It breaks through the limitations of fixed laboratory scenarios and is suitable for continuous monitoring in outdoor and daily multi-scenario environments.

[0144] 2. This invention proposes a cross-modal self-calibration mechanism, which uses the total vertical reaction force of the system level derived from the IMU as a priori to perform dynamic calibration of the vertical force of the left and right feet of the insole through a sliding window, and introduces zero force during the swing period, non-negativity during the support period and smoothness constraints, which significantly suppress the drift and scale deviation caused by long-term wear and multi-tasking of the insole, and provide a stable and reliable priori vertical load for subsequent three-dimensional solution.

[0145] 3. In solving the three-dimensional ground reaction force, this invention unifies the consistency of multi-source information and the hard constraints of the cyber-physical model into the same optimization framework. It adopts system-level conservation equation constraints and friction cone feasible region constraints to solve the static indeterminate problem of horizontal component force distribution in the dual-support stage, and outputs physically feasible, interpretable and robust complete three-dimensional ground reaction forces of the left and right feet.

[0146] 4. This invention forms an end-to-end closed loop from kinematic reconstruction to reverse solving of multidimensional biomechanical parameters from the musculoskeletal model. Driven by wearable data, it can further obtain multidimensional indicators such as joint torque and muscle activation, realizing a comprehensive evaluation from external load to joint load and muscle activation. This avoids the limitations of existing methods that only focus on a single indicator or rely on a large amount of training data, which lack physical consistency and interpretability.

Claims

1. A wearable, sensor-driven, multidimensional musculoskeletal digital twin modeling method, characterized in that, Includes the following steps: Step 1: Calculate joint angles based on the inertial measurement unit; Step 2: Cross-modal calibration of pressure insole sensors; Step 3: Solving for 3D ground reaction forces based on physical constraints; Step 4: Calculate joint torque and muscle activation based on the musculoskeletal model.

2. The wearable sensor-driven multidimensional musculoskeletal digital twin modeling method according to claim 1, characterized in that, Step 1 specifically includes: Based on the attitude information collected by multiple inertial measurement units throughout the body, the joint angles between adjacent segments of the human body are calculated. set up: Represents the global coordinate system; Indicates the first An inertial measurement unit coordinate system; Indicates the first The individual human body segment coordinate system has a one-to-one correspondence between the inertial measurement unit and the human body segment; Indicates time; Represents quaternion multiplication; Indicates quaternion conjugation; Indicates matrix transpose; (1.1) Establish a fixed attitude relationship between the inertial measurement unit and human body segments; After wearing the inertial measurement unit, the subject maintained a neutral, stationary posture; during the period of neutral, stationary posture... Inside, obtain the first The attitude quaternions of each inertial measurement unit relative to the global coordinate system are calculated, and the average attitude quaternion is calculated: in, Indicates time No. The attitude quaternion of an inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; Indicates the first The average attitude quaternion of each inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; Indicates quaternion averaging; superscript Indicates a neutral and static posture; Average attitude quaternion Convert to average direction cosine matrix: in, Indicates the first The average direction cosine matrix of each inertial measurement unit relative to the global coordinate system in a neutral, stationary attitude; This represents the operation of converting a quaternion into a direction cosine matrix; At the same time, according to the first The anatomical coordinate system of an individual body segment in a neutral, static posture yields the first... The cosine matrix of the reference direction of an individual body segment relative to the global coordinate system ;Depend on and Calculate the first The inertial measurement unit to the first Fixed rotation matrix of individual body segment coordinate system : Among them, superscript Indicates a reference state; Fixed rotation matrix Convert to a fixed rotation quaternion: in, Indicates the first The inertial measurement unit to the first Fixed rotational quaternions for an individual anthropometric coordinate system; This represents the operation of converting a rotation matrix into a quaternion; (1.2) Calculate the dynamic global attitude of the inertial measurement unit; During the dynamic movement of the human body, the first... The angular velocity vector acquired by the first inertial measurement unit in its own coordinate system is used to calculate the angular velocity vector of the second inertial measurement unit. Rate of change of dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system Then, the dynamic attitude quaternion for the next time step is obtained through integration: in, equal The first derivative with respect to time; Indicates time No. The angular velocity vector collected by each inertial measurement unit in its own coordinate system; Indicates time t No. The dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system; Indicates time No. The dynamic attitude quaternion of each inertial measurement unit relative to the global coordinate system; This represents the quaternion normalization operation, used to ensure that the updated dynamic attitude quaternions are still unit quaternions. Indicates the amount of change over time; (1.3) Obtain the human segment posture from the posture of the inertial measurement unit; Based on the fixed rotation quaternion obtained in step (1.1) Update the obtained step (1.2) Dynamic attitude quaternions of each inertial measurement unit Convert to the first Dynamic posture quaternions of individual body segments : in, Indicates time No. The attitude quaternion of the individual body segment coordinate system relative to the global coordinate system; Represents a fixed rotation quaternion The conjugate quaternion, due to For a unit quaternion, its conjugate quaternion is... It is equivalent to its inverse quaternion; (1.4) Calculate the relative posture between adjacent human body segments; For any joint, let the proximal human segment of that joint be... The distal human body segment is Based on the quaternion of proximal human segment posture and distal human segment posture quaternions Calculate the relative attitude quaternion of the joint. : in, Indicates time The quaternion of joint relative posture of distal human segments relative to proximal human segments; The conjugate quaternion is used to transform the attitude relationship in the global coordinate system to the proximal human segment coordinate system; (1.5) Convert the relative joint posture quaternions into joint angles; Let the joint relative pose quaternion obtained in step (1.4) be... for: in, Indicates time Scalar components of the joint relative posture quaternion; , , Representing time respectively The three vector components of the joint relative pose quaternion; according to Rotation sequence, joint relative pose quaternion Converted to joint internal / external rotation angles, joint adduction / abduction angles, and joint flexion / extension angles: （10） （11） （12） in, Indicates time The angles of internal and external rotation of the joint; Indicates time The angle of adduction and abduction of the joint; Indicates time The angle of joint flexion and extension; Represents the arctangent function. Represents the arcsine function; Equations (9) to (12) are applied to all joints of the human body to be calculated; assuming the human body has a total of One joint to be calculated Indicates the joint number. For the first The three angles obtained from equations (10) to (12) for each joint are denoted as follows: , and ; (1.6) Output the angle vector of all joints in the body; The first The three angles of each joint form a joint angle subvector. Then, all joint angle subvectors are combined according to a preset joint order to form a whole-body joint angle vector. : in, Indicates time No. Joint angle subvectors of each joint; , and Let each represent the result obtained in step (1.5). Each joint at any time The internal and external rotation angles, adduction and abduction angles, and flexion and extension angles; Indicates the number of joints to be calculated; superscript This represents the transpose of a vector.

3. The wearable sensor-driven multidimensional musculoskeletal digital twin modeling method according to claim 2, characterized in that, Step 2 specifically includes: Step 2 is based on the whole-body joint angle vector obtained in Step 1. The musculoskeletal model is driven to calculate the vertical acceleration of the center of mass of human segments and obtain the total vertical inertial reaction force of the human body according to the Newton-Euler method. Then, the total vertical inertial reaction force of the human body is used as the physical reference to calibrate the vertical ground reaction force of the left and right feet output by the pressure insole sensor. (2.1) Calculate the total vertical inertial reaction force of the human body; Based on the whole-body joint angle vector obtained in step 1 , driving the musculoskeletal model to calculate the first Vertical acceleration of the center of mass of an individual body segment ;by human segment mass Vertical acceleration and gravitational acceleration constant Calculate the total vertical inertial reaction force of the human body. : in, Indicates the total number of human body segments; subscript Indicates the vertical direction of the global coordinate system; (2.2) Establish a vertical ground reaction force calibration model; The pressure insole sensor synchronously outputs the original vertical ground reaction force of the left foot. and the original vertical ground reaction force of the right foot ; A scale factor was introduced for the left and right feet respectively. , and zero drift item , The calibrated vertical ground reaction force of the left foot was obtained. and the reaction force perpendicular to the ground of the right foot : Left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot That is, the three-dimensional ground reaction force of the left foot in the global coordinate system. The components in the direction and the three-dimensional ground reaction force of the right foot in the global coordinate system Components in direction; (2.3) Solve for the calibration parameters within the sliding window; Let time be The center sliding window is The time frame within the sliding window is The calibration parameters to be solved are denoted as a decision variable vector. And the change in scale factor between adjacent windows is denoted as and : The total vertical inertial reaction force of the human body obtained in step (2.1) As a physical benchmark, establish the cross-modal calibration optimization function for the pressure insole sensor. : in, Represents time frame The residual weights; Represents the regularization weights for zero-drift terms; The scale factor represents the smoothing regularization weight; Equation (17) outputs the scale factor. , and zero drift item , Used in equation (15); Equation (17) satisfies the following constraints: in, Represents time frame The left foot supporting state variable, This indicates that the left foot is in the support phase. This indicates that the left foot is in the swing phase; Represents time frame The right foot support state variable, This indicates that the right foot is in the support phase. This indicates that the right foot is in the swing phase; This represents the minimum permissible value of the scaling factor; This indicates the maximum permissible value of the scale factor; This represents the maximum permissible value of the absolute value of the zero-drift term; (2.4) Output the vertical load ratio of the left and right feet; Based on the calibrated left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot Calculate the vertical load ratio of the left foot and the ratio of vertical load on the right foot : (19)。 4. The wearable sensor-driven multidimensional musculoskeletal digital twin modeling method according to claim 3, characterized in that, Step 3 specifically includes: Step 3 is based on the whole-body joint angle vector output in Step 1. Step 2 output: calibrated left foot vertical ground reaction force and the reaction force perpendicular to the ground of the right foot Left foot vertical load ratio and the ratio of vertical load on the right foot By combining the pressure center trajectory output by the pressure insole sensor, the three-dimensional ground reaction force of the left and right feet is calculated; (3.1) Calculate the total inertial resultant force vector of the human body; Based on the whole-body joint angle vector obtained in step 1 , driving the musculoskeletal model to calculate the first Three-dimensional acceleration vector of the centroid of an individual segment ;by human segment mass and three-dimensional acceleration vector Calculate the total inertial resultant force vector of the human body. : in , and These represent the total inertial force of the human body in the global coordinate system. , , Components in direction; (3.2) Calculate the total inertial moment vector of the human body; Based on the centroid position vector of human segments Linear momentum vector Inertial tensor and angular velocity vector Calculate the total inertial moment vector of the human body. : in, express The first derivative with respect to time; Indicates time No. The inertial tensor of an individual body segment about its center of mass; Indicates time No. The angular velocity vector of an individual body segment; express The first derivative with respect to time; Represents the cross product of vectors; The vector of total inertial moment of the human body Decomposed into three directional components: in, , and Representing time respectively Total inertial torque of the human body in the global coordinate system , , Components in direction; (3.3) Calculate the torque arm components of the left and right feet; Pressure insole sensor output Foot pressure center coordinates The musculoskeletal model is based on the whole-body joint angle vectors from step 1. Calculate the coordinates of the human body's center of mass ;Depend on and Calculate the first Components of the left and right foot torque arms from the foot pressure center to the body's center of mass : in, Indicates foot separation, Indicates the left foot. Indicates the right foot; Indicates the directions of the three coordinate axes of the global coordinate system; (3.4) Calculate the resultant force and resultant moment of the three-dimensional ground reaction force during human movement; Let the three-dimensional ground reaction force vector of the left foot be... The three-dimensional ground reaction vector of the right foot is : in, Indicates time Three-dimensional ground reaction force vector of the left foot; Indicates time The three-dimensional ground reaction force vector of the right foot; , , These represent the three-dimensional ground reaction forces of the left foot in the global coordinate system. , , Components in direction; , , These represent the three-dimensional ground reaction forces of the right foot in the global coordinate system. , , Components in direction; The three-dimensional ground reaction forces of the left and right feet satisfy the resultant force equilibrium relationship: The resultant moment components are calculated from the three-dimensional ground reaction forces of the left and right feet and the moment arms of the left and right feet. , and : （26） (3.5) Initialize the horizontal ground reaction forces of the left and right feet; During the dual-support phase, based on the left foot vertical load ratio obtained in step 2... and the ratio of vertical load on the right foot Initialize the horizontal ground reaction forces for the left and right feet: in, , , and These represent the initial values ​​of the horizontal ground reaction forces of the left and right feet, respectively; superscript Indicates the initial value to be optimized; (3.6) Establish friction cone constraints; Establish friction cone constraints for the three-dimensional ground reaction forces of the left and right feet: in, This indicates the coefficient of friction between the shoe sole and the ground; (3.7) Solve for the three-dimensional ground reaction forces of the left and right feet; Using the horizontal ground reaction forces of the left and right feet as the unknowns, a three-dimensional decision variable vector for solving the ground reaction force is established. : By minimizing the calculated total torque With Newton-Euler total moment of inertia The differences between them are used to establish an optimization function for solving the three-dimensional ground reaction force. : By solving equation (30), the three-dimensional ground reaction vector of the left foot is obtained. and the three-dimensional ground reaction vector of the right foot : (31)。 5. The wearable sensor-driven multidimensional musculoskeletal digital twin modeling method according to claim 4, characterized in that, Step 4 specifically includes: Step 4 is based on the whole-body joint angle vector output in Step 1. and the three-dimensional ground reaction force of the left foot output in step 3 Three-dimensional ground reaction force of the right foot The joint torque is calculated using a musculoskeletal model, and muscle activation is solved by static optimization under the constraint of joint torque balance. (4.1) Calculate the joint moment vector using inverse dynamics; Vector of joint angles throughout the body Three-dimensional ground reaction force of the left foot and the three-dimensional ground reaction force of the right foot Input a musculoskeletal model; under the condition of three-dimensional ground reaction force as an external load, calculate the joint moment vector according to the dynamic equations of the musculoskeletal model. : in, This represents the first derivative of the vector of angles of all joints with respect to time. This represents the second derivative of the vector of angles of all joints with respect to time. Represents the musculoskeletal model in terms of the joint angle vectors throughout the body. The generalized mass matrix below; The Coriolis force and centrifugal force correlation matrix representing the musculoskeletal model; Represents the gravity term vector of the musculoskeletal model; (4.2) Establish a balance between muscle force and joint torque; Let the first A muscle at all times Muscle activation level is The maximum isometric strength of this muscle is This muscle is relative to the first The lever arm of each joint degree of freedom is ;No. One muscle against the first The torque contribution of each joint degree of freedom is The sum of the torque contributions of all relevant muscles equals the joint torque of that degree of freedom. : in, Indicates the number of muscles involved in balancing joint torques; (4.3) Solve for muscle activation through static optimization; Combining all muscle activation values ​​into a muscle activation vector and with neuromuscular control cost function With the goal of minimizing, solve for the muscle activation level at each time step: in, This indicates the range of values ​​for muscle activation. (4.4) Output multidimensional biomechanical evaluation results; Joint torque vector and muscle activation vector Combined with the whole-body joint angle vector output in step 1 And the three-dimensional ground reaction vector of the left foot output in step 3. and the three-dimensional ground reaction vector of the right foot The multidimensional biomechanical evaluation results of the musculoskeletal digital twin system were obtained. : (35)。

Images