Muscular strength acquiring method and device based on musculoskeletal model

[0001] The present invention relates to a method capable of acquiring a physically and physiologically valid muscular strength from motion data based on a musculoskeletal model. “Physically valid” means herein that a muscular tension satisfies a dynamic equation. “Physiologically valid” means that a muscular tension is actually generated by a human.

[0002] The somatic sensation is a sensation perceived by the receptors in the superficial layer and deeper tissue of the body. It is divided into tactile cutaneous sensation and deep sensation generated by motor organs such as muscles and tendons. Deep somatic sensation accompanying exercises plays an important part in medical fields such as rehabilitation and surgery, and also sports science. Additionally, it is considered that machines should have the capability of estimating the somatic sensation of the human body for smooth interaction between people and machines.

[0003] The first step in estimating the somatic sensation from motions is to estimate the muscle force from motion data. Conventionally, there have been many studies to estimate muscular force from motion data. However, in these studies the muscular force was calculated using the geometrical nature of the motion such as the physiological model of muscles and percentage change of muscular length. They did not consider whether the motion could be realized with the obtained muscular force. A method for calculating the muscular force considering dynamic validness includes the AnyBody Project (J. Rasmussen, M. Damsgaard, E. Suruma, S. T. Christensen, M. de Zee, and V. Vondrak: “AnyBody-a software system for ergonomic optimization,” Fifth World Congress on Structural and Multidisciplinary Optimization, 2003) and studies of the inventors (Y. Nakamura, K. Yamane, I. Suzuki, and Y. Fujita: “Dynamic Computation of Musculoskeletal Human Model Based on Efficient Algorithm for Closed Kinematic Chains,” Proceedings of International Symposium on Adaptive Motion of Animals and Machines, SaP-I-2, 2003). However, they did not consider the physiological validness. The number of muscular tensions physically valid for a single motion is innumerable. However, they are not always physiologically valid. Conventionally, there have been methods for calculating muscular tensions that are satisfactory in either a physical manner or physiological manner, but there has been no method for calculating the muscular tension satisfying the physical and physiological characteristics at the same time.

Patent document: Japanese published unexamined application 2003-339673

[0004] The present invention provides a method and a device capable of acquiring a physically and physiologically valid muscular strength from motion data based on a musculoskeletal model.

[0005] A technical means adopted by the present invention to solve the problems shown above is a method and a device for obtaining muscular tension by performing inverse dynamics calculation of a musculoskeletal model using the equation shown below.
τ_G=J^Tf+J_C^Tτ_C
where τ_G is generalized force, J is Jacobian of muscles, tendons and ligaments and J_c is Jacobian of contact point.

[0006] The present invention obtains muscular tension by providing reaction force data, motion data and myogenic potential data to the above equation and optimizing a contact force τ_C and muscle tension f.

[0007] In a preferred mode, a method for calculating muscle tension comprises a first step for optimizing a contact force τ_C received from an environment using acquired reaction force data and a second step for optimizing the muscle tension f using acquired motion data, acquired myogenic potential data, and optimized contact force.

[0008] In a preferred mode, the motion data are obtained by a motion capture system. In a preferred example, a motion capture system for acquiring motion data is an optical motion capture system. A motion capture method used for acquiring motion data is not limited to an optical type. In one mode, reaction force data are acquired with a force sensor. If used reaction force data are ground reaction force data, in a preferred mode, the ground reaction force data are acquired with a force plate. Ground reaction force data may be acquired with a force sensor attached to the soles of a subject's feet. In a preferred mode, the myogenic potential data are acquired with an electromyograph. In a preferred mode, the motion data, ground reaction force data and myogenic potential data are measured at the same time.

[0009] In a preferred mode, contact force is optimized by quadratic programming using a certain objective function, and inequality constraint conditions. Muscular tension is optimized by linear programming. In the above equation, generalized force τ_G and Jacobian J of length of muscle, tendon and ligament concerning generalized coordinate can be calculated by acquired motion data. Jacobian J_C of a contact point can be calculated by acquired reaction force data.

[0010] In a preferred mode, the present invention comprises means for storing the obtained motion data and calculated muscular tension by linking them with each other.

[0011] In addition to the method and device, the present invention can also provide a computer program for causing a computer to operate to acquire muscular tension by performing inverse dynamics calculation of musculoskeletal model.

[0012] In the present invention, physically and physiologically valid muscular tension can be calculated by performing inverse dynamics calculation and optimization calculation using motion data, measurement data of myogenic potential data and ground reaction force. If muscular tension is calculated for many tasks and subjects using this method, and similarity of antagonistic muscles is found for each task, muscular tension can be estimated to a certain degree even if there is no myogenic potential and ground reaction force data.

[0019] The present invention provides a method for calculating a physically and physiologically valid muscular strength from motion data obtained by motion capture and information about myogenic potential and ground reaction force measured at the same time. Human musculoskeletal model consists of many muscles and is generally a redundant drive system. Therefore, muscle force cannot be determined uniquely. If standard usages of muscles in various motions are compiled in a database by using this method, muscular force may be estimated even if myogenic potential and ground reaction information are not acquired.

[0020] In embodiments of the present invention, a device for acquiring muscle tension comprises a plurality image capturing means (a camera) to capture subjects with marks affixed, ground reaction force measuring means (a force plate), myogenic potential measuring means (an electromyograph) and one or more computers. A computer comprises a processing unit, an input device, an output device, a display device and a memory device storing various data.

[0021] In this embodiment, a behavior capture system measures electromyogram (EMG) and ground reaction force, as well as motion data, at the same time. Physically and physiologically valid muscular strength is acquired by considering them in optimization of muscular strength. The behavior capture system is disclosed in Seri, Yamane, Nakamura: “Simultaneous and real-time measurement of conscious actions with behavior capture system”, Lecture meeting on robotics and mechatronics, The Japan Society of Mechanical Engineers, proceedings 01, 1P1-H7, 2001.

[0022] The behavior capture system is a system for measuring interactive information of human behavior in addition to a motion capture system that measures only human whole-body motions. The behavior capture system can measure and process information input to humans (point of gaze), information output from humans (behavior, ground reaction force, brain waves and internal human information) in simultaneous and in a real-time manner. The information can be reviewed in a 3D manner on a screen. In the present invention, the behavior capture system measures the motion capture data (motion data), electromyogram (EMG), and ground reaction force at the same time. A physically and physiologically valid muscular strength is acquired by using the data in optimization of muscular strength.

[0023] Musculoskeletal model and inverse dynamics of musculoskeletal, which are prerequisites for calculating muscular strength, will now be explained. As a human musculoskeletal model, a model developed by a group of inventors consisting of a total 547 of muscles, tendons and ligaments, and a bone structure with a degree of freedom of 155 is used (FIG. 1). In this model; muscles, tendons and ligaments are expressed as a wire having no thickness. Kinetic properties such as via-points and divergence points are based on anatomical documents. Details of the musculoskeletal model are disclosed in Suzuki, Yamane, Taya, Kurihara, Nakamura: “Human dynamic computation based on detailed human body model using high-speed kinetics calculations”, collected manuscripts, the 20th academic lecture of The Robotics Society of Japan, 3C14, 2002, and Japanese published unexamined applications 2003-339673, and 2004-13474. Data of musculoskeletal model are stored in the memory device of the device of the present invention.

[0024] A designed detailed human model consists of bones grouped with appropriate detail, and network among muscles, tendons and ligaments covering areas between bones. Each element has the following characteristics. A bone is a rigid-body link having mass. A muscle is a wire that generates tension positively and is relatively elastic. A tendon is a wire that generates tension passively and is relatively inelastic. Additionally, a tendon transmits the muscle tension to distant parts by connecting to muscles. A ligament is a wire that generates tension passively and is relatively inelastic. A ligament connects bones, and constrains relative motion among bones.

[0025] Muscles, tendons and ligaments are modeled as wires that pass through start points, end points and zero or more routing points. A routing point is established in cases where a muscle gets caught in a part of a bone, and a model where a tendon is constrained by a tendon sheath. In general, the start point, end point and routing point are fixed on bones. A virtual link is established at a bifurcation point to describe diverging muscles and ligaments. In this way, a wire starts from a start point, passes zero or more routing points and ends at an end point.

[0026] Inverse dynamics will now be explained. In inverse dynamics, the tension of muscles, tendons and ligaments to realize a motion is calculated based on motion data acquired from the motion measurements or motion generation such as choreography. The flow of inverse dynamics calculation includes: 1. Calculation of joint torque by inverse dynamics of a rigid-body link system; 2. Calculation of Jacobian for joint value of wire length; and 3. Conversion of joint torque to wire tension.

[0027] The inverse dynamics of this model can be formulated as follows. First, generalized force TG required for a given motion is calculated. If generalized coordinate 0 is a joint angle of a bone structure, the generalized force corresponds to joint torque. Therefore, the generalized force (joint torque) can be calculated using Newton-Euler Formulation. In this case, assuming that the six degrees of freedom of the hip link are also driven, τ_G includes six axial forces acting on the hip link. Then, Jacobian J of the length of muscles, tendons and ligaments relating to generalized coordinates is calculated. In other words, if the length of muscles, tendons and ligaments is made aligned vector 1, J is expressed by the following equation.

[0028] $\begin{array}{cc}J\stackrel{△}{=}\frac{\partial l}{\partial {\theta }_G}& \left(1\right)\end{array}$
Detailed methods for calculating Jacobian J are disclosed in Suzuki, Yamane, Taya, Kurihara, Nakamura: “Human dynamic computation based on detailed human body model using high-speed kinetics calculations”, collected manuscripts, the 20th academic lecture of The Robotics Society of Japan, 3C14, 2002, in Japanese published unexamined application 2003-339673, and in D. E. Orin and W. W. Schrader: “Efficient Computation of the Jacobian for Robot Manipulators,” International Journal of Robotics Research, vol. 3, no. 4, pp. 66-75, 1984.

[0029] Using Jacobian J, the relationship between generalized force, tension vectors f of muscles, tendons and ligament, and contact force received from environment, is described as

[0030] $\begin{array}{cc}{\tau }_G=J^{T}f+\sum J_{\mathrm{Ci}}^T{\tau }_{\mathrm{Ci}}& \left(2\right)\end{array}$
where τ_Ci is a contact force at the contact point i, J_Ci is Jacobian at the contact point i.

[0031] Jacobian J_c of a contact point is defined by the following equation. In the following equation, pc is a vector or a parameter expressing the position and posture of a contact point. The method for calculating Jacobian is disclosed in the documents shown above.

[0032] $\begin{array}{cc}{J}_c\stackrel{△}{=}\frac{\partial p_c}{\partial {\theta }_G}& \left(1^{\prime }\right)\end{array}$

[0033] Terms concerning contact force in equation (2) are rearranged as follows:
τ_G=J^Tf+J_C^Tτ_C  (3)

[0034] The inverse dynamics calculation of the musculoskeletal model is to obtain f and τ_C satisfying equation (3) from τ_G. Generally, in the musculoskeletal model, the number of muscles is larger compared to number of degrees of freedom. Therefore, f is not determined uniquely. If two or more links contact the environment, the contact forces become indeterminate. It is required to perform optimization by some criteria to determine these values.

[0035] Optimization will now be explained. First, objective function and restraint conditions are mentioned considering measured values of contact force and electromyogram (EMG), followed by description of formulation of optimization problem and solution.

[0036] The contact force will now be explained. As mentioned above, the generalized force τ_G includes six axial forces acting on the hip link that is not actually driven. The muscular strength that is an internal force cannot generate a corresponding force. Therefore, it is substituted by the contact force. The rows corresponding to six degrees of freedom of the hip link extracted from equation (3) is as follows:
E_hipτ_G=E_hipJ_C^Tτ_C  (4)
where E_hip is a matrix that extracts required rows. If the number of contacted links is one, unique solution exists. However, if the number is two or more, it will become indeterminate.

[0037] In the present invention, this problem is solved by using ground reaction force data measured with a force plate. Therefore, the following terms are included in the optimization objective function.
|K_Cτ_C−τ*_C|²  (5)
where τ*_C is a measured contact force, and K_C is a matrix that maps the contact forces acting on each link into a coordinate system of measured contact forces. For example, if a plurality of links exist on a force sensor, the sum of the ground reaction forces acting on those links is measured, and so is converted by K_C.

[0038] The contact force shall satisfy a certain condition regarding the vertical force, center of pressure, friction force, etc. Among them, the vertical force is handled by incorporating the following inequality constraint condition.
E_vertτ_C≧0  (6)
where E_vert is a matrix that extracts the vertical components from τ_C. Other conditions are considered approximately by incorporating the following equation in the objective function.
|τ_C|²  (7)

[0039] Electromyogram (EMG) will now be explained. In the present invention, the Stroeve muscle model that formulates the muscle model of Hill and Wilkie is used. In this model, the relationship between muscle length and maximum isometric muscular strength is expressed as FIG. 2. The relationship between change rate of muscle strength and maximum muscle tension is shown in FIG. 3. At the maximum contraction velocity v_max, the muscular strength is zero. The maximum muscle tension when the muscle length does not change corresponds to the maximum isometric muscular strength. A muscle extends when a force larger than the maximum isometric muscular strength is applied.

[0040] IEMG (the integrated EMG), which is a value measured by an electromyograph integrated by a certain time interval, expresses the activity factor u of a motor nerve. The relationship between u and muscle activity factor a is expressed by the following equation.

[0041] $\begin{array}{cc}\stackrel{.}{a}=\frac{u-a}{T}& \left(8\right)\end{array}$
where T is a parameter expressing time delay. The relationship between a and muscle tension f* is expressed by the following equation.
f*=aF_l(l)F_v({dot over (l)})F_max  (9)
where F_max is a function to express the maximum muscle tension. F_l(l) and F_v(l with dot) are functions expressing the relationship between normalized muscular strength and muscular length, and change rate of muscular length, respectively.

[0042] F_l(l) corresponds to FIG. 2 and is approximated by the Gaussian function in the following equation.

[0043] $\begin{array}{cc}{F}_l\left(l\right)=\exp \left\{-{\left(\frac{l-l_0}{K_l}\right)}^2\right\}& \left(10\right)\end{array}$
F_v(l with dot) corresponds to FIG. 3 and is approximated by the following equation.

[0044] $\begin{array}{cc}{F}_v\left(\stackrel{.}{l}\right)=\{\begin{array}{cc}0& \left(i\le -v_{\max }\right)\\ \frac{V_{\mathrm{sh}}\left(v_{\max }+i\right)}{V_{\mathrm{sh}}{v}_{\max }-i}& \left(-v_{\max }\le i\le 0\right)\\ \frac{V_{\mathrm{sh}}{V}_{\mathrm{shl}}{v}_{\max }+V_{\mathrm{ml}}i}{V_{\mathrm{sh}}{V}_{\mathrm{shl}}{v}_{\max }+i}& \left(0\le i\right)\end{array}& \left(11\right)\end{array}$
where K_l, V_sh, V_shl, and V_ml are constants. Values mentioned by Stroeve are used here.

[0045] n_ch muscular strength values f* are obtained by using an electromyograph with n_ch channels. As for the surface electromyograph used in the present invention, an electromyogram (EMG) measured in a channel does not always correspond to a muscle in the musculoskeletal model. In some cases, the sum of potential of several muscles is measured. Considering this case, the measured value of electromyogram (EMG) shall be considered by incorporating the following equation in the objective function with the matrix K_F that makes f and measured value correspond to each other.
|K_Ff−f*|  (12)
A muscle, tendon and ligament can only generate force in the direction of tension. Therefore, the following inequality constraint is applied.
f≦0  (13)

[0046] Optimization calculation will now be explained. Optimization calculation is performed in the following two phases by using equations shown above.

[0047] Step 1: Optimization of contact force using quadratic programming

[0048] Step 2: Optimization of muscle tension by using linear programming.

[0049] First, in the step 1, τ_C that satisfies inequality constraint condition
E_vertτc≧0  (14)
and minimizes the objective function
Z=w_H|E_hipτ_G−E_hipJ_C^Tτ_C|²
+w_C|K_Cτ_c−τ_C*|²+|τ_C|²  (15)
is calculated, where w_H and w_C are weighting parameters. The first term in the objective function has the effect of reducing the error of equation (4). w_H=1×10⁸ and w_c=10 are used here.

[0050] When a contact force is acquired by the calculations shown above, the generalized force that should be generated by the muscle tension can be calculated by the following equation.
τ′_G=τ_G−J_C^Tτ_C=J^Tf.  (16)
In the step 2 of the optimization, δ_f, δ_τ and f that satisfy the inequality constraint condition
−δ_f≦K_Ff−f*≦δ_f  (17)
E_mtlf≦0  (18)
δ_f≧0  (19)
−δ_τ≦τ′_G−J^Tf≦δ_τ  (20)
δ_τ≧0  (21)
and minimize the objective function
Z=a_f^Tδ_f+a_τ^Tδ_τ  (22)
are obtained by using τ′_G, where a_f and a_τ are constant vectors with positive components. The muscle tension can be approximated to the measured value by minimizing equation (22) based on the equation (17). Equation (20) reduces the error of equation (16) induced from the dynamics equation. Therefore, dynamic consistency is also considered.

[0051] Modifications to muscular strength optimization calculation will now be explained. The deviation from average muscle tension of a muscle group is added to the objective function for optimizing the muscular strength. In other words, the foregoing objective function
Z=a_f^Tδ_f+a_τ^Tδ_τ  (22)
is replaced by
Z=a_f^Tδ_f+a_τ^Tδ_τ+a_m^Tδ_m  (23)
where a_m is a constant vector with positive components. δ_m is a vector that satisfies the following restraints.
−δ_m≦E_Gf≦δ_m  (24)
0≦δ_m  (25)
E_G is a matrix for calculating the difference between the average value of muscle tensions belonging to each of pre-defined muscle group and each muscle tension. They are calculated in a method shown below.

[0052] Assuming that a class including the number of muscles belonging to muscle group k is G_k, the average f_k (with bar) of muscular tension belonging to the muscle group k is obtained from

[0053] $\begin{array}{cc}\stackrel{\_}{f_k}=\frac{1}{n_k}\sum _{i\in G_k}^{}{f}_i& \left(26\right)\end{array}$
where n_k is a number of muscles belonging to the muscle group k. Consequently, the difference between the tension of the muscle n(nεG_k) and average of muscles belonging to muscle group k is calculated by
Δf_nk= f_k−f_n=E_Gknf  (27)
where E_Gkn is a horizontal vector with the same number of elements as the number of muscles. The i-th element is (1−n_k/n_k) for I=n, 1/n_k for i≠n and iεG_k and 0 in other cases. Finally, E_G is obtained by vertically accumulating E_Gkn of muscles contained in all muscle groups.

[0054] In the conventional calculating method using linear programming, muscles that mainly generate tensions in the muscles (cooperative muscles) that have the similar role (elbow joint extension, etc) change discontinuously. As a result, the method has a problem that the tension of each muscle changes rapidly. Muscle tension is averaged among cooperation muscles and smoother tension change is obtained by changing the objective function as shown above.

[0055] Experimental examples will now be explained. Motion capture markers and electrodes for 16 channels were affixed to both legs of a subject. EMG, ground reaction force, and motion capture data for squat exercise were obtained at the same time. Two force plates were used for measuring the ground reaction force. The force and moment applied to left and right leg were measured independently. As a result, EMG and ground reaction force data per 1 ms and motion capture data per 33 ms were acquired.

[0056] The upper row of FIG. 4 shows the result of optimization for the right lateral great muscle without using the data of EMG and ground reaction force assuming τ*_C=0 and f*=0. Initial rise of muscle tension is delayed compared to the EMG data in the middle row of FIG. 4, indicating that the muscle started activity before starting actual motion. The pattern similar to the measured values is obtained if actual EMG data is used for calculation (lower row of FIG. 4). The computation time required for optimizing the muscle tension is 0.81 s on average.

[0057] If optimization is performed without using the ground reaction force data, the ground reaction force as shown in the upper row of FIG. 5 is acquired because the force is distributed evenly to left and right. On the other hand, with the measured value of elements of the ground reaction force in the vertical direction, the force distribution for left and right legs is uneven after starting activity as shown in the middle row of FIG. 5. Distribution similar to measured value as shown in the lower row of FIG. 5 is also acquired by considering measured data. When both legs were in contact with the ground, the computation time required for optimizing the ground reaction force was 8.8 ms on average.

[0058] FIG. 6 shows the visualized results of the case where no EMG and ground reaction force data was used (upper row) and the case where EMG and ground reaction force data were used (lower row). The larger muscle tension is indicated in black not white. It indicates that unbalance between left and right reaction force is reflected in the muscle tension if ground reaction force data are considered.

[0059] In the present invention, a physically and physiologically valid muscular strength can be calculated by performing inverse dynamics and optimized calculation by using motion data, and measured data of EMG and ground reaction force. If muscular tension is calculated for many tasks and subjects, and similarity of antagonistic muscles is found for each task, muscular tension can be estimated to a certain degree even if there is no EMG and ground reaction force data.