JPWO2010095636A1 - 筋張力推定法及び装置 - Google Patents

Abstract

筋張力推定のための最適化計算の規模を縮小する。複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭｉ（ｉ＝１，２，．．．ｎ）に分類し、各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成し、同じサブグループに属する筋の筋活動度を同じとみなす。複数の筋グループＭｉの一部あるいは全部において、前記１つあるいは複数のサブグループのうちの１つは、筋電計が装着された１つの代表筋と当該代表筋と起始停止する骨が同じである筋とから形成される第１サブグループである。第１サブグループに属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得することで最適化計算における変数を削減し、あるいは、各サブグループを代表する筋活動度を最適化計算で推定することで、最適化計算における変数を削減する。

Description

本発明は、筋張力の推定及び推定された筋張力に基づく身体内部の活動情報の提示に関するものである。

人間の動作を解析するには、人体においてアクチュエータとして働く筋の筋張力パターンを知ることが重要である。例えば、筋肉の活動や筋肉を支配する脊髄神経束の活動を運動データから推定するためには、運動データから筋張力を取得することが必要である。
従来の筋張力推定手法としては、代表的に２つの手法が知られている。

一つの解法として、Hill-Stroeve筋モデル(非特許文献１、非特許文献２)のように筋の収縮特性をモデル化し、生理学や実験より求まるパラメータや筋長、計測したEMGから計算される筋活動度から筋張力を算出する生体力学的なアプローチが挙げられる。比較的簡便に測定可能で最も筋張力を良く表すパラメータである筋電位を用いていることから最も精度の高い推定値であると言える。また計算コストが小さく高速に計算できるという利点もある。しかし、筋張力を計測するためには EMGを測定する必要があり、筋張力を推定する筋すべてに筋電計を装着する必要がある。そのため筋電計の数の制限から筋張力の推定を行う筋の数には限界が生じる。仮に筋電計が無数に利用可能であったとしても、それらをすべての筋に装着することは不可能であり、被験者の行動を制限することにもなる。

もう一つの解法として、ロボティクスの分野で発達した逆動力学計算と最適化計算を用いた動力学的なアプローチがある。この解法では、逆運動学、逆動力学を基に、筋骨格モデルを用いて各関節トルクを算出することで、動力学的に整合性のとれた解を求めることができる。特許文献１、非特許文献３、非特許文献４では、モーションキャプチャにより得られた関節角度や床反力計のデータを用いて逆動力学計算を行って関節トルクを算出し、それらをもとに最適化計算によって筋張力を推定している。
このアプローチでは、モデルの自由度に対してそれらを駆動する筋の数が大きく逆動力学計算のみでは関節トルクから筋張力を一意に定めることはできず、最適化計算を用いて解を1つに絞る必要がある。すなわち、筋張力分配時の最適化計算において線形計画法や二次計画法などの計算が含まれるため計算コストが大きく多くの時間を要するという課題がある。
最適化計算の際の目的関数の中にHill-Stroeve筋モデル及び筋電計で計測する筋活動度から求まる筋張力との誤差項を考慮することで、より精度の高い全身の筋張力の算出が可能となるが（特許文献１、非特許文献３、非特許文献４）、逆動力学計算を用いて筋張力の最適化計算を行うことは、筋張力のリアルタイム推定を考えた場合、計算コストが問題となる。

このように様々な筋張力推定法が提案されているが、未だ完全な方法は確立しておらず、それぞれの方法は一長一短を有している。これら筋張力推定はトレーニング支援やリハビリテーションに応用されることが期待されており、リアルタイムで筋張力の推定を行い、それを理解しやすい形で提示することができれば効率のよい運動支援が行えるものと考えられる。しかし、生体力学的な方法は全身の筋張力を推定することが事実上不可能であり、最適化計算を用いる方法では計算に時間がかかるであるなどの問題を抱えている。そのため従来の方法を用いて全身の筋張力をリアルタイムで推定することは困難である。

国際公開番号ＷＯ２００５−１２２９００

A. Hill, the heat of shortening and thedynamic constants of muscle, Proceeding of the Royal Society of London, Vol.B126, pp. 136-195, 1938.

Sybert Stroeve, impedance characteristicsof a neuromusculoskeletal model of the human arm I: Posture control, Journal ofBiological Cybernetics, Vol. 81, pp. 475-494, 1999.

Y. Nakamura, K. Yamane, Y. Fujita, and I.Suzuki, somatosensory computation for man-machine interface from motion capturedata and musculoskeletal human model, IEEE Transactions on Robotics, Vol.21(1), pp. 5866, 2005.

Y. Nakamura, K. Yamane, and A. Murai,macroscopic modeling and identification of the human neuromuscular network,Proceedings of the 28th IEEE EMBS Annual International Conference, pp. 99-105,2006.

本発明は、筋張力推定のための最適化計算における変数を削減することで、最適化計算の規模を縮小することを目的とするものである。

本発明の他の目的は、運動時の筋張力をリアルタイムで推定すると共に、筋張力および／あるいは筋張力に基づいて取得される身体内部の活動状況をリアルタイムで提示することにある。

本発明が採用した第１の技術手段は、
被験者の運動時の各関節トルクを、筋骨格モデルを用いて逆動力学計算により算出し、該関節トルクを最適化計算により筋張力へ分配することで各筋の筋張力を推定する方法において、
被験者の複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類し、
複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）の各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成し、同じサブグループに属する筋の筋活動度を同じとみなし、
複数の筋グループＭ_ｉの一部あるいは全部において、前記１つあるいは複数のサブグループのうちの少なくとも１つは、筋電計が装着された１つの代表筋と当該代表筋と起始停止する骨が同じである筋とから形成される第１サブグループであり、
前記第１サブグループに属する筋の筋張力を、最適化計算を用いずに被験者の運動時に計測された前記代表筋の筋電位から取得し、最適化計算の対象から外すことで、最適化計算における変数を削減し、
あるいは、
前記１つあるいは複数のサブグループにおいて、各サブグループを代表する筋活動度を最適化計算で推定することで、最適化計算における変数を削減する、
筋張力の推定法、である。
第１の技術手段において、筋電計のチャンネル数と筋グループＭ_ｉの数ｎは必ずしも同じである必要はない。
第１の技術手段において、第１グループが形成されている筋グループでは、典型的には１つの筋グループにつき１つの第１グループが形成されるが、１つの筋グループにおいて２つ以上の代表筋を選択して２つ以上の第１グループを形成してもよい。

１つの態様では、前記第１サブグループに属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得し、
前記第１サブグループに属しない筋の筋張力を、被験者の運動を実現するのに必要な関節トルクにおいて、前記第１サブグループに属する筋により実現できない関節トルクを実現するように最適化することで推定する。

１つの態様では、各筋グループＭ_ｉにおいて、前記１つあるいは複数のサブグループは、前記第１サブグループと、前記第１サブグループに属しない筋について、起始停止する骨が同じ筋から分類された零個以上のサブグループと、を含み、
各サブグループを代表する筋活動度を最適化計算で推定する。
さらに、１つの態様では、前記第１グループに属する筋については、筋電位から取得した筋活動度と最適化計算により計算される筋活動度が一致すべきであるとして、
最適化計算において、計測された前記代表筋の筋活動度を参照値として用いる。

本明細書において、「筋の運動指向性に基づく筋のグルーピング」とは、言い換えると、「各関節の各運動方向に関わる筋のグルーピング」あるいは「各関節において筋が関与するトルクの向きに基づくグルーピング」ということができる。これらは、運動と筋のつき方の幾何学的な関係を表すものである。具体的には、例えば、「肘関節伸張に関わる筋」、「股関節屈曲に関わる筋」のような意味である。
すなわち、筋グループＭ_ｉへのグルーピングは、肘関節を伸ばす（伸展）ための筋、肘関節を曲げる（屈曲）ための筋といったように各関節の運動の方向に関わる筋の役目による分類である。筋グループＭ_ｉへのグルーピングを行うことによって、起始停止する骨のみによる筋のグルーピングで、関節を曲げるための筋と関節を伸ばすための筋とが同じグループに分類されてしまうことが防止される。
ある関節を屈曲させる方向に寄与する（寄与する関節が全く同じ）筋の筋グループＭ_ｉには、起始停止する骨が違う複数の筋が属しており、これをさらに起始停止する骨によってサブグループに分類する。筋グループＭ_ｉへのグルーピングを行うことによって、厳密に関与する関節に基づいてさらに下位のサブグループへの分類が可能となる。
また、筋グループＭ_ｉへのグルーピングは、多関節筋の影響を無視した運動指向性による筋の分類、サブグループへのグルーピングは多関節筋の影響も考慮し、筋が収縮した時に直接運動に影響を与える関節の一致する筋を同じグループに分類するということができる。
「異名筋促通に基づくグルーピング」は、同じ神経束で支配される筋同士をグループとして表したものである。
筋をグルーピングする際、力学的な面からと神経生理学的な面からの２つのアプローチがある。前者は、筋の走行から決定されるもの、後者は筋同士の神経による結合から決定されるものである。スポーツ科学の分野において、これら2つの間に関係があることは示されている。
サブグループへのグルーピングは、「運動学に基づいて協同筋に分類」しているということができる。「協同筋」とは、「ある関節について、主動筋と同じ方向へ曲げるのに働く筋」と定義することができる。「筋の運動指向性」、「協同筋」は、両方とも運動学の見地からの考えであり、上述のように「異名筋促通」は神経生理学の見地からの考えである。
１つの態様では、筋のグルーピングは、「筋の運動指向性」、「協同筋」に従って行うが、例えば、屈曲伸展に関しては協同筋でも、内転外転では拮抗筋になる場合等があるため、必要に応じて、「異名筋促通」でグルーピングを確認することが望ましい。
また、協同筋の筋グループM_iへの分類は体の姿勢や関節の変位によって変化することが知られている。近似的には分類は変化しない静的なものとして扱えるが、さらに精度を高めるためには分類を体の姿勢や関節の変位によって動的に変化させることもできる。

本発明が採用した第２の技術手段は、方法の発明としては、
被験者の運動時の各関節トルクを、筋骨格モデルを用いて逆動力学計算により算出し、該関節トルクを最適化計算により筋張力へ分配することで各筋の筋張力を推定する方法において、
被験者の複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類し、
複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）の各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成し、同じサブグループに属する筋の筋活動度を同じとみなし、
前記１つあるいは複数のサブグループにおいて、各サブグループを代表する筋活動度を最適化計算で推定することで、最適化計算における変数を削減する、
筋張力の推定法、である、

また、第２の技術手段は、装置の発明としては、
被験者の運動時の各関節トルクを、筋骨格モデルを用いて逆動力学計算により算出し、該関節トルクを最適化計算により筋張力へ分配することで各筋の筋張力を推定する筋張力取得手段と、
被験者の複数の筋を分類するグルーピング手段と、
を備え、
前記グルーピング手段は、
前記複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類する第１グルーピング手段と、
複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）の各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成する第２グルーピング手段と、からなり、
前記筋張力取得手段は、同じサブグループに属する筋の筋活動度を同じとみなし、前記１つあるいは複数のサブグループにおいて、各サブグループを代表する筋活動度を最適化計算で推定する、
筋張力の推定装置、である。

第２の技術手段において、１つの態様では、複数の筋グループＭ_ｉの一部あるいは全部において、前記１つあるいは複数のサブグループのうちの少なくとも１つは、筋電計が装着された１つの代表筋と当該代表筋と起始停止する骨が同じである筋とから形成される第１サブグループであり、
前記第１グループに属する筋については、筋電位から取得した筋活動度と最適化計算により計算される筋活動度が一致すべきであるとして、最適化計算において、計測された前記代表筋の筋活動度を参照値として用いる。

本発明が採用した第３の技術手段は、方法の発明としては、
被験者の全身あるいは身体の一部の複数の筋の少なくとも一部を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類し、各筋グループＭ_ｉから１つの代表筋を選択して当該代表筋に筋電計を装着し、
前記複数の筋を、
各筋グループＭ_ｉの前記代表筋からなる第１筋群Ｍ_ｉＥＭＧと、
各筋グループＭ_ｉにおいて、第１筋群Ｍ_ｉＥＭＧと起始停止する骨が同じである（寄与する関節が全く同じである）筋からなる第２筋群Ｍ_{ｉｈｉｇｈ}と、
前記第１筋群Ｍ_ｉＥＭＧ、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれない筋からなる第３筋群と、
に分け、
第１筋群Ｍ_ｉＥＭＧと第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得し、
前記第３筋群に属する筋の筋張力を、逆動力学計算により、計測した被験者の運動を実現するのに必要な関節トルクを計算し、前記関節トルクにおいて、前記第１筋群Ｍ_ｉＥＭＧ及び前記第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋により実現できない関節トルクを実現するように最適化することで推定する、
筋張力の推定法、である。

また、第３の技術手段は、装置の発明としては、
複数の筋電計と、
筋電位から筋張力を取得する第１筋張力取得手段と、
逆動力学計算により関節トルクを計算し、当該関節トルクを実現するように最適化計算を行うことで筋張力を推定する第２筋張力取得手段と、
被験者の複数の筋を分類するグルーピング手段と、
を備え、
前記グルーピング手段は、第１グルーピング手段と第２グルーピング手段とを備え、
前記第１グルーピング手段は、前記複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類するものであり、各筋電計は、各筋グループＭ_ｉから選択された１つの代表筋に装着されており、
前記第２グルーピング手段は、前記複数の筋を、
各筋グループＭ_ｉの前記代表筋からなる第１筋群Ｍ_ｉＥＭＧと、
各筋グループＭ_ｉにおいて、第１筋群Ｍ_ｉＥＭＧと起始停止する骨が同じである筋からなる第２筋群Ｍ_{ｉｈｉｇｈ}と、
前記第１筋群Ｍ_ｉＥＭＧ、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれない筋からなる第３筋群と、
に分けるものであり、
前記第１筋張力取得手段は、第１筋群Ｍ_ｉＥＭＧと第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得し、
前記第２筋張力取得手段は、前記第３筋群に属する筋の筋張力を、逆動力学計算により、計測した被験者の運動を実現するのに必要な関節トルクを計算し、前記関節トルクにおいて、前記第１筋群Ｍ_ｉＥＭＧ及び前記第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋により実現できない関節トルクを実現するように最適化することで推定する、
筋張力の推定装置、
である。

第３の技術手段には幾つかの態様がある。
１つの態様では、前記複数の筋において、いずれの筋グループＭ_ｉにも属しない筋を筋群Ｍ_{ｏｔｈｅｒｓ}とし、
前記第３筋群には、各筋グループＭ_ｉにおいて、前記第１筋群Ｍ_ｉＥＭＧ、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれない筋群Ｍ_ｉｌｏｗと、前記筋群Ｍ_{ｏｔｈｅｒｓ}と、が含まれる。
１つの態様では、前記複数の筋を、全ての筋がいずれかの筋グループＭ_ｉに属するように分類し、
前記第３筋群には、各筋グループＭ_ｉにおいて、前記第１筋群Ｍ_ｉＥＭＧ、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれない筋が含まれる。
１つの態様では、各筋グループＭ_ｉにおいて、筋群Ｍ_{ｉｈｉｇｈ}に属する複数の筋同士が拮抗しない。また、多関節筋に関しては、１つの態様では、主に関与する関節の方に含むようにグルーピングされる。
１つの態様では、前記第１筋群Ｍ_ｉＥＭＧと前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれる筋の筋張力は、筋電位データから得られる筋活動度、経験則により得られた筋のパラメータ、計測した運動データに基づく逆運動学計算により取得される筋長及び筋長の変化速度、から取得される。
１つの態様では、各筋グループＭ_ｉにおいて、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれる筋の筋電位は、前記代表筋の筋電位の関数として決定する。
典型的には、第２筋群Ｍ_{ｉｈｉｇｈ}に含まれる筋の筋電位は、前記代表筋の筋電位と同じと見なすが、第２筋群Ｍ_{ｉｈｉｇｈ}の筋電位は第１筋群Ｍ_ｉＥＭＧの筋電位と全く同じである必要はなく、これらの筋同士の幾何学的な位置、姿勢等から得られる関数により決定さ得る。
１つの態様では、筋電計のチャンネル数と筋グループＭ_ｉの数ｎが同じである。
筋電計のチャンネル数と筋グループＭ_ｉの数ｎが一致しない場合、例えば筋電位のチャンネル数が少ない場合には、筋グループＭ_ｉの任意のグループに対してＭ_ｉＥＭＧを設定し、筋電計を装着する。Ｍ_ｉＥＭＧが設定されなかった筋グループＭ_ｉに関しては、当該筋グループＭ_ｉに属する全ての筋がＭ_ｉlow、結果的に第３筋群に属する筋、となる。すなわち、EMGのチャンネル数などの制限から全てのM_iにおいてM_iＥＭＧが決定できるわけではなく、その場合はM_iの全ての筋は第３筋群に属することになる。

１つの態様では、
前記最適化計算は、
前記複数の筋群が骨格を駆動するとして、各筋群が出力すべき関節トルクを筋群毎に推定する第１ステップと、
各筋群において、前記第１ステップで推定された関節トルクを実現するように、各筋が出力する筋張力を推定する第２ステップと、
からなる。
最適化計算の問題の規模をEMG情報を用いて解決できる筋の数だけ減らすという考えは、一部の筋張力を求める場合にも有利である。このときにEMGが解決する筋を決めるためにグルーピングが使われるが、残る最適化計算で解決される筋の次元をさげることにもグルーピングが役立つ。
最適化手法には様々な方法があり、計算コストが大きく拘束条件を正確に満たすものと、計算コストが低く拘束条件を曖昧にしか満たさないものがある。ここでいう拘束条件とは、筋が伸張方向には力を出さないという不等式拘束条件を意味する。ただし後者の場合においても、最適化の際に考慮する筋同士の関係が拮抗するものを含まない場合、拘束条件を曖昧にしか満たさなくとも正確な解を得ることができる。
そこで、まず筋をグルーピングすることで、複数の筋グループが骨格を駆動するとして、前者の計算コストが大きい手法にて各筋グループが出力すべきトルクを求める。
そして、各筋グループにおいてそれらのトルクを実現するよう、後者の計算コストが小さい手法にて各筋が出力する筋張力を求めることができる。
筋の次元とは、大雑把に言うと、筋の本数を意味する。より具体的には、筋の次元とは、全身の筋の張力を決める上での独立な変数の数をさす。筋を一本一本独立に張力が決まるとする場合には、筋の総数が次元になり、この場合は大規模な最適化問題になる。筋をグループ化してグループ内での張力の分配規則を決めておくと、全身の筋張力を決める上で独立な変数は、グループの数になる。この場合の筋の次元は筋のグループの数になる。

１つの態様では、筋張力を、被験者の運動時に実時間で推定する。上記筋張力推定装置は、筋張力を、被験者の運動時に実時間で推定するリアルタイム筋張力推定装置である。

本発明が採用した第４の技術手段は、方法の発明としては、
被験者の撮影画像あるいは／および当該撮影画像に基づく合成画像を表示部に表示すると共に、表示された被験者の画像に筋骨格モデルをオーバーレイし、
上記推定法により取得した筋張力に基づく身体内部の活動情報を筋骨格モデルに反映させて視覚的に表示する、
身体内部の活動情報提示法、である。
また、第４の技術手段は、装置の発明としては、
リアルタイム筋張力推定装置と、
運動時の被験者を撮影する手段と、
被験者の撮影画像あるいは／および当該撮影画像に基づく合成画像を表示する表示手段と、
を備え、
前記表示手段に表示された被験者の画像に筋骨格モデルをオーバーレイし、前記リアルタイム筋張力推定装置により実時間で推定した筋張力に基づく身体内部の活動情報を前記筋骨格モデルに反映させて、視覚的に実時間表示するように構成されている、
身体内部の活動情報提示装置、である。
１つの態様では、身体内部活動情報を、被験者の運動時に実時間で表示する。

１つの態様では、前記身体内部の活動情報は、筋活動である。
１つの態様では、筋活動を、筋骨格モデルの筋の色あるいは／および形状の変化によって視覚的に表示する。

１つの態様では、前記身体内部の活動情報は、筋活動を、当該筋活動を支配する脊髄神経束の活動として表わしたものである。脊髄神経は複数の筋を支配しており、支配する全ての筋の活動をまとめたものを脊髄神経束の活動として計算している。
１つの態様では、脊髄神経束の活動は、筋骨格モデル上の各脊髄神経束の位置にシンボルを表示し、シンボルの色あるいは／および形状の変化によって視覚的に表示する。

本発明では、被験者の複数の筋をグルーピングすることによって、筋張力推定のための最適化計算における変数を削減し、最適化計算の規模を縮小する。
最適化計算の規模を縮小することによって、リアルタイムに筋張力を推定し、また、筋張力に基づく筋活動等をリアルタイムで可視化することがでる。
本発明によれば、EMG情報（Hill-Stroveモデルと共に利用）とグルーピングによって最適化計算の規模を縮小することができる。
さらに、限られた数の筋電計を用いるものでありながら、計算量が比較的小規模の最適化計算手法を組み合わせることで、全身の筋張力をリアルタイムで推定することができる。
したがって、スポーツトレーニングやリハビリテーション時に体性感覚情報を提示することも可能となる。筋張力や関節負荷を中心とする体性感覚情報をリアルタイムに提示することで、トレーナの代替となるシステムの構築が可能になる。
人間の筋肉の活動やそれを支配する脊髄神経束の活動を運動データから高速に推定し、それを実時間で本人の映像にオーバーレイするなどの方法で提示することによって、身体内部の活動状況を直感的に理解させることがきる。自分自身の運動とその時の身体内部の活動を直感的に見ることができ、運動の効果を確認できる。
スポーツトレーニング、リハビリテーション、医療診断、健康管理、エンターテイメントなどにおいて、運動の効果を確認させることができる。

筋長と最大等尺性筋力との関係を示す図である。 筋長の変化率と最大筋張力との関係を示す図である。 人体の筋骨格モデルを示す図である。 従来の最適化計算による筋張力推定の流れを示す図である。 リアルタイム筋張力可視化システムの概略図である。 リアルタイム筋張力可視化システムの概略フロー図である。 リアルタイム筋張力可視化システムにおける各スレッドとデータの受け渡しの概略図である。 被験者の運動時に、被験者の画像上に、筋骨格モデル及び推定された筋張力情報をリアルタイムでオーバーレイした図を示す。上図は、スクワット、下図は投球動作である。実際には、筋張力が大きくなると、筋の色が黄色から赤色に変化するようになっている。 図６のスクワットを複数フレームの時系列で示す図である。 図６の投球動作を複数フレームの時系列で示す図である。 表示された筋骨格モデルに脊髄神経束をシンボルで示す図である。画像中において各脊髄ごとに四方に配置されている球が脊髄神経束をシンボル化したものであり、各シンボルは該当脊髄神経束が活動すると赤色に変化する。 本発明の第１の実施形態に係る筋のグルーピングを示す図である。 本発明の第２の実施形態に係る筋のグルーピングを示す図である。

本発明の実施形態について詳細に説明する。先ず、本発明の背景となる概念および手法について説明する。これらの概念および手法は本発明の背景技術であると同時に、本発明の実施形態を実施する上でも用いられ得る技術である。次いで、本発明に係る筋張力推定の実施形態について説明する。尚、数式の番号についてはセクション毎で独立に付与している。

［Ａ］Hill-Stroeve筋モデル
HillとWilkieの筋モデルを定式化したStroeveの筋モデルでは、筋長と最大等尺性筋力の関係は図１のように表される。また、筋長の変化速度と最大筋張力の関係は図２のように表される。最大収縮速度ｖ_ｍａｘでは筋張力は０になり、筋長が変化しない時の最大筋張力が最大等尺性筋力に対応する。また、最大等尺性筋力より大きい力が加わった場合には筋は伸張する。

筋活動度ａと筋張力ｆ^＊の関係は式(1)で表される。
ここで、Ｆ_ｍａｘは最大筋張力、F_l(l)とF_v(l（ドット）)はそれぞれ正規化された筋張力と筋長、筋長の変化速度との関係を表す関数である。F_l(l)は図１に対応し、式(2)のガウス関数で近似する。
ここで、l₀は筋の自然長である。また、F_v(l（ドット）)は図２に対応し式(3)で近似する。
ここで、K_l，V_sh，V_shl，V_mlは定数であり、１つの態様では、Stroeveが示した値(表１)を用いる。また、これらの値をモーションキャプチャデータに基づいて同定してもよい。

Hill-Stroeve筋モデルを用いて筋張力推定をする場合に必要なデータとして、筋長、筋長変化速度、筋活動度がある。全ての筋の筋長、筋長変化速度、筋活動度が得られた場合、各筋の筋張力ｆ^＊は次の式で表される。
ここで、a_i，l_i，l（ドット）_i，F_maxは、各々i番目の筋の活動度、筋長、筋活動度、最大筋張力を表し、F_l，F_vはそれぞれ正規化された筋張力と筋長、筋長変化速度の関係を表す関数である。そして、Nmusは筋骨格モデルに含まれる筋の総数を表す。筋長l₁,……l_Nmus及び筋長変化速度l(ドット)₁,……l（ドット）_Nmusについては、全てモーションキャプチャから得られる運動データから算出できる。

筋活動度については、体幹、肘、肩、膝、足首等体の各部位において筋電位を用いた筋張力推定のモデルが考案されている。これらのモデルは、例えば、下記の文献に記載されている。

S.M.
McGill, A myoelectrically based dynamic three-dimensional model to predict
loads on lumbar spine tissues during lateral bending, Journal of Biomechanics,
Vol. 25, pp. 395-414, 1992.

T.S. Buchanan,
S.L. Delp, and J.A. Solbeck, muscular resistance to varus and valgus loads at
the elbow, Journal of Biomechanics, Vol. 120, pp. 634-639, 1998.

B.
Laursen, B. Jenson, G. Nemeth, and G. Sjogaard, A model predicting individual
shoulder muscle forces based on relationship between electromyographic and 3D
external forces in static position, Journal of Biomechanics, Vol. 31, pp.
731-739, 1998.

David G.
Lloyd, Thomas S. Buchanan, and Thor F. Besier, Neuromuscular Biomechanical Modeling
to Understand Knee Ligament Loading, Medicine & Science in Sports & Exercise,
Vol. 37, pp. 1939-1947, 2005.

A.L. Hof
and J.W. van den Berg, EMG-to-force processing. II. Estimation of parameters of
the Hill Muscle model for the human triceps surae by means of calf ergometer,
Journal of Biomechanics, Vol. 14, pp. 759770, 1981.

これらのモデルにおいて、各関節に対して複数の筋の筋電位を筋電計で計測し、その%MVC(Maximal Voluntary Contraction(MVC、最大随意努力時の筋活動)を100%として表した力の入れ具合を表す値)等をその筋の活動度として使用している。

［Ｂ］逆動力学・最適化計算を用いた筋張力算出
全身詳細筋骨格モデルにおける逆動力学計算の手法を示す（特許文献１、非特許文献３、非特許文献４）

［Ｂ―１］筋骨格モデル
本発明の実施形態で用いられる全身詳細筋骨格モデルについて述べる。図１に示すように、設計した詳細人体モデルは、適当な細かさでグループ分けされた骨格系剛体モデルと、骨格上に張られた筋・腱・靭帯系ワイヤモデルとからなる。骨格モデルは全身206個の骨からなる。そのうち頭蓋部、手部、足先部などは一つの剛体として扱い、計53個のリンクからなるモデルとなっている。各リンク間は、足根骨-足先部の回転1自由関節、第１胸椎-胸骨の6自由度関節を除いて全て球面3自由度関節となっている。骨格モデルは、全体の並進回転の6自由度を加えて、計155の自由度を持つ。

次に骨格モデルに筋、腱、靭帯を配置する。筋、腱、靭帯は各リンクに始点、終点及び経由点を通るワイヤとしてモデル化する。骨、筋、腱、靭帯はそれぞれ以下の性質を持つ。
骨：質量を持つ剛体リンク
筋:能動的に張力を発生するワイヤである。
腱:受動的に張力を発生するワイヤで、筋と接続し筋張力を骨へ伝達する。
靭帯:受動的に張力を発生するワイヤで、骨と骨とを接続し、それらの相対的な運動を拘束する。
また筋、腱、靭帯の機能の違いは、以下のようにモデル化する。
筋と腱の直列接続からなるような簡単な部位は、1本の筋ワイヤで代表する。
筋が骨の一部分に引っ掛かっている場合や腱鞘による腱の拘束をモデル化する場合には経由点を置く。
上腕二等筋など腱が分岐し、分岐した腱がそれぞれ別々の骨に接続するという配置になっている場合がある。ワイヤの始点、終点、経由点は全てリンクに固定されるため、この分岐点にヴァーチャルリンクを置く。ヴァーチャルリンクは質量を持たないが張力を伝達する。ヴァーチャルリンクは力、モーメントが0になるように自由に移動できる。
大胸筋や広背筋等の広い筋は、複数の並行な筋ワイヤで表現する。
このような筋骨格モデルについては、例えば特許文献１にも記載されており、この文献を参照することができる。
上述の筋骨格モデルは、例示に過ぎないものであり、本発明に適用される得る筋骨格モデルは、これらに限定されるものではない。

［Ｂ−２］筋骨格モデルを用いた筋張力の取得
筋骨格モデルを用いた筋張力の取得について説明する。一つの態様では、筋張力の取得装置は、マーカが付された被験者を撮影する複数の撮像手段（カメラ）と、床反力計測手段（フォースプレート）と、筋電位計手段（筋電位計）と、一つ又は複数のコンピュータ装置とを含み、コンピュータ装置は、各種計算を行う演算処理部、入力部、出力部、表示部、各種データを格納する記憶部を備えている。ここでは、モーションキャプチャデータ（運動データ）、筋電位、床反力を同時計測し、これを筋力の最適化において用いることで、力学的にも生理的にも妥当な筋力を得る。

全身詳細筋骨格モデルの筋張力計算について説明する。
特許文献１、非特許文献３、４に開示された方法では以下のように筋張力を計算する。
（１）モーションキャプチャシステムにより被験者の運動計測を行い、マーカの三次元位置の時系列データを得る。
（２）逆運動学計算によりマーカの三次元位置から関節角、関節角速度、関節角加速度を含む運動情報を計算する。
（３）ニュートンオイラ法などを用いた逆動力学計算により運動を実現するのに必要な関節トルクを計算する。
（４）関節角から得られる筋、腱、靭帯長変化と各関節角速度の関係を用いて（３）で求めた関節トルクを、床反力及び筋、腱、靭帯の張力に写像する。

逆動力学では、運動計測によって得られる運動データを元に、その運動を実現する筋・腱・靭帯の張力を求める。逆動力学の計算法の流れは、１．剛体リンク系の逆動力学による関節トルクの計算；２．ワイヤ長さの関節値に対するヤコビアンの計算；３．関節トルクのワイヤ張力への変換、となる。
剛体リンク系の逆動力学計算を用いると骨格モデルにおいて運動を実現するのに必要な関節トルクτ_ｇが計算できる。ダランベールの原理と仮想仕事の原理を用いるとτ_ｇと等価な筋、腱、靭帯張力ｆは、関節角θ_ｇに対する筋、腱、靭帯長ｌのヤコビアンＪを用いて、
と表される。
ヤコビアンＪの計算方法については、当業者によく知られているので、説明が煩雑になることを避ける目的で、ここでの詳述は省略する。ヤコビアンＪの計算方法については、例えば、特開２００３−３３９６７３号、あるいは、「D.E. Orin and W.W. Schrader. Efficient computation of the jacobian
for robot manipulators. Inter-national Journal of Robotics Research, Vol. 3,
No. 4, pp. 66.75, 1984」を参照することができる。

本実施形態において関節角ベクトルθ_Ｇは155次元であるのに対し、ワイヤ張力fは非特許文献３、非特許文献４のモデルでは989次元である。そのためτ_ｇからfが一意には定まらない冗長問題が生じる。ここで、筋骨格モデルの逆動力学計算において、運動を決定するパラメータに対して筋・腱・靭帯の要素数が非常に多く、力が一意に決まらないという未決定性問題が存在することは当業者に良く知られており、逆動力学計算により求められた関節モーメントを、最適化計算よって、各関節を駆動する筋の筋張力へ分配することが行なわれている。

ｆを決定するために、何らかの評価関数と拘束条件を設定し、数理計画法等による最適化を用いて解決する方法は、例えば、特許文献１、非特許文献３、４に開示されている。以下に最適化計算の具体例を示すが、筋張力計算に用いられる最適化計算としては幾つもの手法が提案されていることは当業者に知られており、本発明に適用され得る最適化計算は、本明細書に記載されたものに限定されないことは当業者に理解される。

全身詳細筋骨格モデルにおける逆動力学計算の１つの態様を説明する。以下の式を解くことにより、床との接触力及び全身の筋張力を算出する。
ここで、
τ_Ｇ：一般化力；
Ｊ：一般化座標からワイヤ長へのヤコビ行列；
ｆ：ワイヤ張力；
Ｊ_Ｃ：一般化座標から床との接触点へのヤコビ行列；
τ_Ｃ：床との接触力；
である。

式(4)を以下の流れで解く。
式(4)の内、腰関節の6DOFに対応する行のみを考慮して、床との接触力τ_Ｃを算出する。ここでは2次計画法を用いて最適化を行う。
式(4)からτ_Ｃを除き、下式を得る。
ここで線形計画法もしくは2次計画法を用いて筋張力を算出する。
床との接触力は、外界との接触力の典型例であり、床以外、例えば壁との接触力を用いることもできる。このような外界との接触力を差し引くことによる筋張力推定については、特許文献１、非特許文献３、非特許文献４に開示されている。

以下に、接触力τ_Ｃの最適化を考慮した関節トルクτ_Ｇ´の求め方について述べる。筋張力を求める複雑な最適化は2次計画法、または更に単純化して線形計画法で解くことができる。ここでは、高速な解法として線形計画法を用いた方法を考える。しかし、線形計画法特有の問題として、結果として得られる筋張力データが時間方向・空間方向に不連続になる。空間的不連続は、幾何学的に近い筋同士において筋張力が大きく異なることを示し、これは人体においては不自然である。

この問題を、協同筋グループ内における筋張力の平滑化によって解決する。以下の式を用いて、各要素が正であるベクトルａ_τ，ａ_ｆ，ａ_ｍに対して、次式である目的関数を最小化するδ_τ，δ_ｆ，δ_ｍ，ｆを求める最適化問題を線形計画法を用いて解く。
ここで、拘束条件は以下のように書ける。
式(6)の第3項は筋張力の平滑化のために付加している。以下式(6)から式(13)について説明する。

式(6)の第1項及び式(7)、式(8)は、式(5)における誤差の最小化を目的とし、動力学的な整合性を保障している。式(5)は等式の形で書けるが、式(5)が解を持たない場合を考慮して条件を緩和している。目的関数式(6)がａ_τ ^Ｔδ_τを含むことで、式(8)により正に拘束された最小のδ_τが得られる。一方、式(7)により、式(5)の誤差をδ_τより小さくしている。これらの拘束条件を考慮することで式(5)の誤差を最小化することができる。
式(6)の第2項及び式(9)、式(10)は、ｆを与えられた目標値ｆ^＊に近づける効果を持つ。例えば、適当な値のｆ^＊を与えることで、屈筋・伸筋間の筋力の関係を一意に定めることができる。バイオメカニカルなアプリケーションとして筋電計の測定値を用いるなどが考えられる（特許文献１、非特許文献３、非特許文献４）。ｆ^＊＝０とすると、最小の筋・腱・靱帯張力が得られる。

式(11)は筋・腱・靱帯張力の上限と下限の拘束であり、ｆ_max≧０は最大筋張力を表す。最大筋張力ｆ_maxは各筋ごとに独立に決めることができる。Hill-Stroeve筋モデルを用いて、筋長とその変化速度からｆ_maxを算出し、最大筋張力を与えることもできる。

最後に、式(6)の第3項及び式(12)、式(13)は協同筋群内において可能な限り筋張力を平滑化する効果を持つ。この効果は、2次計画法においては筋張力の2乗和の項を目的関数に付加することで実現できる。ｎ_ｍ本の筋からなるｍ番目の協同筋群Ｇ_ｍを考える。
この協同筋群内の筋張力の平均値は、
で算出される。ここでｆ_kはｋ番目の筋の筋張力を表す。ｋ(ｋ∈Ｇｍ)番目の筋の筋張力とそれが含まれる協同筋群における平均筋張力の差は、
で表すことができる。ここでＥ_Ｇｍｋは、i番目の要素が、
である行ベクトルである。全ての協同筋群についてのＥ_Ｇｍｋ(ｋ∈Ｇ_ｍ)を並べることで、式(12)に示すＥ_Ｇが得られる。

上記の例では、線形計画法を用いた最適化方法について説明したが、次に、二次計画法による筋張力の最適化方法の１つの態様を示す。不等式拘束条件つきの二次計画法「M. Renouf and P. Alart. Conjugate gradient type algorithms forfricional
multi-contact problems: Applications to granular materials. Vol. 194, pp. 2019.2041,
2005」に基づき、評価関数Ｚを、
として、Ｚを最小にするｆを求める。なお、Ｋ_１、Ｋ_２は重みである。τ_G´は、τ_GからJ_c ^Tτ_cが差し引かれた一般化力である。
これにより、筋張力を計測値に近づけることができる。また（τ_G´−Ｊ^Ｔｆ）も小さくなるので、力学的にも妥当な筋張力が計算できる。

［Ｃ］リアルタイム全身筋張力推定法
［Ｃ−１］リアルタイム筋張力可視化システム
図５に、リアルタイム筋張力可視化システムの概略図を示す。筋張力可視化システムは、筋張力推定手段と、推定された筋張力を用いて身体内部の活動情報を取得する手段と、運動時に撮影された被験者の画像及び推定された筋張力／取得された身体内部の活動情報を表示する手段と、を備えている。より具体的には、筋張力可視化システムは、身体の複数の所定部位に複数のマーカが付された被験者を撮影する複数の撮像手段（カメラ１）と、運動中の被験者を表示手段に表示するために撮影する撮影手段（DVカメラ２）と、床反力計測手段（フォースプレート３）と、筋電位計などの筋電位計手段（無線筋電位計４）と、一つ又は複数のコンピュータ装置５と、表示手段（スクリーン６）と、を含む。コンピュータ装置は、各種計算を行う演算処理部、入力部、出力部、表示部、各種データを格納する記憶部を備えている。モーションキャプチャデータ（運動データ）、筋電位、床反力を同時計測し、これを筋力の最適化において用いることで、力学的にも生理的にも妥当な筋張力を取得する。

図５Ａに、リアルタイム筋張力可視化システムのフロー図を示す。人体の運動データ、床反力は光学式モーションキャプチャ及びフォースプレートを用いて計測し、筋電位データは無線筋電計を用いて計測する。それぞれのデータはシステム制御PCに同期してリアルタイムで取得される。その後、運動データから筋骨格モデルの逆運動学計算（ＩＫ）を解くことで、関節角、筋長、筋長変化速度が得られる。このようにして、Hill-Stroeve筋モデルに基づく筋の筋張力の取得及び逆動力学計算（ＩＤ）及び最適化計算を用いた筋張力の推定に必要な情報を取得する。推定された筋張力は、筋骨格モデル上に配置された筋の色の変化（例えば、筋張力が大きくなるにしたがって黄色から赤色に変化させる）で可視化される。更に、同期してDVカメラを用いて実際の実験風景を撮影し、視点を合わせた上で上記筋骨格モデルに重ね合わせて可視化する。図６に実際に提示された映像のスクリーンショットを示す。

［Ｃ−２］筋張力推定の全体構成
本実施形態に係る筋張力取得は大きく分けて次の２つの工程を備えている。
先ず、EMGデータを用いて、EMG電極が装着された筋の筋張力、及び、この筋に密接に関連する筋の筋張力を求める。
次いで、筋張力ｆと関節トルクτ_Ｇ´との関係
を用いて他の筋の筋張力を推定する。また、未知数の数を低減することに加えて、EMGデータは、下記制約(3)を満足させる解の効率的な推定を可能とする。

表２に、筋骨格モデルの要素数及び自由度を示す。
左列は体幹の脊柱起立筋群等細部の筋まで全てモデル化した従来の解析用のモデルである。最適化計算の目的関数は、
であり、不等式拘束条件
を満たすように最適化計算が行われる（非特許文献３参照）。
表２において、右列のsimplified modelは、左列の詳細なモデルであるcomplex modelから重要度の低い要素を間引いたモデルである。表２に示すモデルは例示であって、本発明がこれらのモデルに限定されるものではない。また、後述する第２実施形態では、simplified modelにおける筋の本数を２７４から３１４に増加させている。
本発明の実施形態では、リアルタイム推定のために計算コスト低減を目的として、simplified modelをさらにグルーピングによって要素数を減らして低次元化することで、最適化計算の規模を縮小している。実際、本実施形態では、１６ｍｓ（従来の最適化計算に比べて１０倍以上の速度）で筋張力を推定することができる。

［Ｃ−３］筋のグルーピング
ここで神経生理学の分野における筋同士の神経結合について考える。神経結合は介在ニューロンを介した筋同士の結合であり、促通性と抑制性が考えられている。促通性の結合を持つ筋は協同筋として働き、抑制性の結合を持つ筋は拮抗筋として働く。この神経結合の機能的意義について考える場合、筋の作用ごとに分類し、例えば肘の屈筋群と伸筋群などの筋群(協同筋)に分けて論じられることが多い。これは、異名筋促通や拮抗筋抑制の考えを前提にすると、同一筋群の中では促通性の、拮抗作用を示す筋群の間では抑制性の結合が予想でき、筋群にまとめることで神経結合の機能をより単純化して考えることができるからである。

このことから、本実施形態では、全身の筋を異名筋促通のグループに分類し、グループ内で１つの筋を代表筋として選択して筋電計を設け、代表する筋の筋電位を計測し、その活動度をグループ内の全筋の活動度の代表値とする。本実施形態では、一般的な無線筋電計のチャンネル数である16chで全身の筋を計測する。

筋の分割については、四肢の関節の運動に注目し、肩関節(水平)屈曲/伸展、肘関節屈曲/伸展、股関節屈曲/伸展、足首関節背屈/底屈に分類する(膝関節については、股関節と連動する部分が多いため、まとめる)。各グループにおいて代表する筋を、三角筋前部、三角筋後部、上腕二頭筋長頭、上腕三頭筋外側頭、大腿直筋、大腿二頭筋長頭、前脛骨筋、腓腹筋外側頭とする。各筋の%MVCを計測し、その筋が含まれる異名筋促通のグループ全ての筋の筋活動度とし、モーションキャプチャから得られた筋長及び筋長変化速度から筋張力を求める。
本態様では、手足の運動に着目し、身体の左右のそれぞれの５つの関節（合計１０個）について考える。さらに、筋は表３に示すように、８つの群に分けられる。

各筋グループM_i(i＝1, 2, . . . , 8)のそれぞれの筋は、さらに、次の３つの組に分けられる。
１．Ｍ_ｉＥＭＧ:それらのEMG信号が測定される代表筋の筋群。
２．Ｍ_{ｉｈｉｇｈ}: 筋群Ｍ_ｉＥＭＧと寄与する関節が全く同じである筋からなる筋群。
３．Ｍ_ｉｌｏｗ: 筋群Ｍ_ｉＥＭＧと寄与する関節一部において同じである筋からなる筋群。
Ｍ_ＥＭＧ、Ｍ_ｈｉｇｈ、Ｍ_ｌｏｗは、それぞれ、
Ｍ_ＥＭＧ＝Ｍ_１ＥＭＧ∪Ｍ_２ＥＭＧ∪…∪Ｍ_８ＥＭＧのように規定する。
表３のいずれの筋グループにも属しない筋をＭ_{ｏｔｈｅｒｓ}とする。

これらの筋の筋張力を以下の流れで取得する。
まず、Ｍ_ＥＭＧ、Ｍ_ｈｉｇｈに含まれる筋については、筋電位及びHill-Stroeve筋モデルから筋張力を取得する。
そして、逆動力学計算及び最適化計算によって、残りの筋、すなわち、Ｍ_ｌｏｗ、Ｍ_{ｏｔｈｅｒｓ}に含まれる筋の筋張力を推定する。

［Ｃ−４］Hill-Stroeve筋モデルを用いた筋張力の取得
Hill-Stroeve筋モデルを用いて筋張力推定をする場合に必要なデータとして、筋長、筋長変化速度、筋活動度がある。全ての筋の筋長、筋長変化速度、筋活動度が得られた場合、各筋の筋張力ｆ^＊は次の式で表される。
ここで、ａ_ｉ、ｌ_ｉ、Ｆ_ｍａｘｉは、それぞれ、筋肉ｉの筋活動度、筋長、筋長の変化速度、最大等尺性筋力、であり、F_l(＊)、F_v(＊)は、筋張力と筋長、筋張力と筋長の変化速度、を表す関数である。関節角及び速度を用いた順動力学計算によって、l_i、l（ドット_）iが与えられる。

筋活動度a_iの進展を、１次微分方程式でモデル化する。
ここで、Ｔは時定数であり、u_iはMVCによって正規化されたEMG信号から計算される運動神経からの入力である。
ＥＭＧ信号からの筋活動度の算出には幾つかのやり方が当業者に知られており、例えば、以下の文献に記載された算出法を用いることができる。
S.Stroeve. Learning combined feedback and feedforward
control of a musculoskeletal system. Biological Cybernetics, Vol. 75, pp. 73.83,
1996.

各グループの代表筋の筋張力は、式(4)から直接計算することができる。
筋ｋ∈Ｍ_{ｉｈｉｇｈ}は、以下の式によって、同じグループの代表筋ｒ∈Ｍ_ｉＥＭＧの活動度から推定することができる。
ここで、E_{r→k（＊）}は、グループM_ｉに含まれる筋ｋの活動度a_iと、計測された代表筋の活動度ａ_ｒの関係を表す。E_{r→k（＊）}の関数としては、Georgepoulosらにより主張されるコサインチューニングに従う方法等が考えられるが、ここでは次式で定義する。

［Ｃ−５］逆動力学計算及び最適化計算による筋張力推定
ここまで、第１筋群、第２筋群の筋張力が取得され、したがって、最適化のための未知数の数を削減することができる。しかしながら、依然として、不等式拘束条件を備えた最適化計算の計算コストは大きい。以下に、不等式拘束条件を用いない効率的な筋張力推定について説明する。

先ず、式(1)における筋張力ｆ、ヤコビアンＪを各筋群に分配する。
ここで、J_EMG，J_high，J_low，J_othersは、関節角に対する各M_EMG，M_high，M_low，M_othersにおける筋長のヤコビ行列、f_EMG，f_high，f_low，f_othersは、各群における筋張力である。
τ_Ｇ´は、既に床反力τ_Ｃが差し引かれた一般化力である。
さらに、以下のように変形する。
ここで、
未知数の数は低減されているが、不等式拘束条件ｆ≦０を備え、したがって、反復計算を伴う最適化計算を行う必要がある。
ここで、singularity-robust
(SR) inverse [NAKAMURA, Y., AND HANAFUSA, H. 1986.“Inverse Kinematics Solutions
with Singularity Robustness for Robot Manipulator Control”. Journal of Dynamic Systems,
Measurement, and Control 108, 163.171.]を用いることで、不等式拘束条件を用いないで、反復計算を伴わない最適化計算を提案する。

先ず、式(6)を用いて、M_lowの筋張力の初期推定値を求める。
ここで、ｋ∈M_ilow、ｒ∈M_iEMGである。筋ｒと筋ｋは、同じ異名筋促通グループに属するので、式（１１）は初期値としては適当であると言える。

ついで、逆動力学計算によって得られた関節トルクを用いて、筋張力を補正していく。アルゴリズムは、各関節ｊについて、以下の工程を繰り返す。

ステップ１：関節ｊを駆動し、かつ、M_low及びM_othersに属する全ての筋を集めて、筋群M_Jjを形成する。また、式（９）における関節ｊ、筋群Ｍ_Ｊｊにおける筋に対応する行列を抽出して、
とする。

ステップ２：全ての筋ｋ∈M_Jjについて、以下の式によって、f_jk0の初期推定値を取得し、
全てのｋに対してｆ^＊ _jk0を補正することで、ｆ^＊ _j0を形成する。

ステップ３：式(12)からＪ^Ｔ _ｊｆ^＊ _ｊ０を差し引いて、次の式を得る。

ステップ４：Ｊ^Ｔ _ｊのSR-inverse、Ｊ^Ｔ＊ _ｊを用いて、Δf_j0を解き、次式によって更新されたｆ_ｊを得る。
SR-inverseは、
を最小化するものであり、ここで、ｔは正の重みである。Δf_j0の要素は、十分に小さく、f_j1が正とはならないことが予想される。
もし、f_j1≦０が保持されれば、ｆ_ｊ＝f_j1として終了する。
さもなければ、ステップ５に進む。

ステップ５：大きさがf_j1と同じで、要素がf_j1の最大値（正）であるベクトルf_j1max を形成する。第２の推定値f_＊j1＝f_j1―f_j1max≦0で計算する。

ステップ６：式(12)からＪ^Ｔ _ｊｆ^＊ _ｊ０を差し引いて、次の式を得る。

ステップ７：再び、J^T _jのSR-inverseを用いてΔf_j1を解き、次式によってf_jを更新する。
上記ステップ４における議論と同様に、f_j2の多くの要素は負であり、正であったとしても、少なくとも、小さいことが予想される。したがって、f_j2をf_jの近似として用いることができる。

上記アルゴリズムはJ^T _jのSR-inverseのみを用いるものであり、関節ｊを駆動する筋のみ考慮すればよいので、J_jのサイズは小さいものとなる。したがって、このアルゴリズムは、反復計算による最適化計算に比べて高速である。
本実施形態に係る筋張力推定法によれば、全身の筋張力推定に要する時間は16msであり、体性感覚情報の可視化に要した時間は68msであった（用いた計算機は、3.33 GHz Intel Xeon processor (3.25 GB RAM, NVIDIA Quadro FX3700)である）。結果として、15fpsのフレーム速度の視覚化システムが構築できた。

［Ｃ−６］体性感覚情報のリアルタイム提示
本発明の実施形態では、被験者の撮影画像あるいは／および当該撮影画像に基づく合成画像を表示部に表示すると共に、表示された被験者の画像に筋骨格モデルをオーバーレイし、上記推定法により取得した筋張力に基づく身体内部の活動情報（体性感覚情報）を筋骨格モデルに反映させて視覚的に表示する。
自分の動きに同期して動く自分の実写映像の上に筋活動の画像を被せることで、いかにも自分の体の中が透けて見えているような感覚を自然に感じさせる、すなわち、体内の体性感覚を透視できているという状況を実現することができる（図６参照）。

ここで、筋が活動する際の流れについて説明する。脊髄内のα運動ニューロンが、上位中枢からの運動指令信号あるいは筋の固有感覚受容器からの信号により興奮性に活動する。この信号は「脊髄神経束」を通って脊髄外へ出る。α運動ニューロンから各筋上の終板構造に向かい「脊髄神経束」が分化して信号が伝達される。終板構造において興奮性の信号がまとめられ、筋に活動電位を与える。これが「筋肉の活動」に繋がる。筋上を活動電位が伝達されていく際に「筋張力」が生じる。したがって、筋活動度や脊髄神経束の活動をリアルタイムで取得して提示するためには、筋張力をリアルタイムで取得することが重要である。

より具体的には、ビデオカメラにより撮られる映像の上に筋骨格モデルをオーバーレイして表示し、その筋骨格モデル上の筋の色や形状により筋活動度を表現する。筋張力/筋活動度は筋の色を、例えば、黄色→赤色に連続的に変化させることで表示する。したがって、オーバーレイにおいて、筋が活動している部分は、例えば赤色で表示される。また、筋の疲労を筋の太さで表現してもよい。１つの態様では、筋の疲労として、筋張力を時間積分した値を用いる。

脊髄神経束の活動は、各脊髄神経束の位置にシンボル、例えば球、を表示し、その色で活動を表現する（図９参照）。図９では、便宜上、表示部に表示された筋骨格モデルにおいて脊髄神経束の位置に球を表示しているが、これをさらに被験者の画像にオーバーレイすることができる。
脊髄神経束の活動の推定については、例えば、下記の文献を参照することができる。
Murai, A, Yamane, K, and Nakamura,
Y, "Modeling and Identification of Human Neuromusculoskeletal Network Based
on Biomechanical Property of Muscle," the 30^th IEEE EMBS Annual
International Conference, pages 3706-3709, Vancouver, August 2008.

［Ｄ］筋張力推定の第２実施形態
第２実施形態における装置の全体構成や全体の処理の流れは、筋のグルーピング及び最適化計算を除いて、第１実施形態と同じであり、図５、図５Ａに示す通りである。まず光学式モーションキャプチャと逆運動学計算により被験者の関節角度を得る。次に逆運動学から求まった関節角度、フォースプレートから得られる床反力から逆動力学計算を行って関節トルクを得る。最後にEMGや関節トルクから筋のグルーピングに基づいて最適化計算を行い筋張力の推定を行う。

生体における反射の一つとして、骨格筋を収縮させる体性反射がある。筋の固有受容器は筋長やその変化速度、筋張力の変化を感知する。筋長やその変化の速度は筋紡錘によって、また筋張力はGolgi腱器官によって感知される。固有受容器は中枢の運動調節に2つの側面から貢献している。1つは反射効果であり、もう1つは固有受容器からの情報が運動や姿勢の状態を上位脳に伝える事である。運動指令で運動が起こると筋の状態は変化し、固有受容器反射は必然的に誘発される。その反射効果は運動プログラムに反映され、運動パターンの形成と修正に寄与する。
筋紡錘に由来する反射には伸張反射、拮抗抑制、α−γ関連などが挙げられるが、伸張反射に注目する。伸張反射とはある主動筋の筋紡錘に由来する神経線維の発火は、その筋を支配する運動ニューロン
(同名筋運動ニューロン
)とその協同筋の運動ニューロンに単シナプス性の興奮が引き起こされることである。
伸張反射により、或る動作に必要な主動筋に上位中枢から信号が送られると、その筋のIa群神経線維から戻ってきた信号が協同筋の運動ニューロンに興奮効果を与え、協同筋に信号が送られる。その結果、関節トルクを出す為に協同筋の間で筋張力の分配が行われる。よって協同筋同士は同程度の筋活動度を持つことが期待される。そこで協同筋に注目し、協同筋同士は同程度の筋活動度を持つものとしグルーピングを行うことで計算の低次元化を図る。ただし、解剖学に基づきIa群神経線維と運動ニューロンの接続から協同筋は明らかにされていない。本実施形態では、解剖学ではなく運動学の面から協同筋について考え、協同筋のグルーピングを行う。

［Ｄ−１］筋のグルーピング
以下の記述において、i，jを筋のグループのインデックス、kを筋のインデックスとして用いる。
まず四肢と体幹の動きに注目し、全身の筋を異名筋促通のグループM_i(i=1,2,...,n_G)に分類する。n_GはM_iのグループ数である。各グループはさらに以下の2つのグループに分類される。

１．M_ihigh：グループを代表するEMGが計測されている筋、及び、EMGが計測されている筋と起始停止する骨が同じである筋。同じグループに属する筋は筋活動度が同程度であることが期待される。第２実施形態におけるM_highは、第１実施形態におけるＭ_ＥＭＧ＋Ｍ_ｈｉｇｈに対応する。
２．M_ilow：上記のM_ihighに属さない筋。さらに起始停止している骨によってM_i,1low,..., M_i,nilowに分類する。ただし、n_iはM_ilow内のM_i,jlowのグループ数である。

各グループM_iにはEMGを計測される筋は2つ以上存在しないものとする。さらに、M_high =M_1high∪M_2high∪...∪M_nGhigh，M_low=M_1low∪M_2low∪...∪M_nGlowによりM_high，M_lowを定義する。なお、筋電位を計測している筋の本数をn_EMG(≦n_G)としたとき、筋電位を計測している筋が含まれているグループはM₁,...,M_nEMGであり、M_(nEMG+1)high,..., M_nghighは空集合である。なお、本実施形態ではn_G=36，n_EMG=16とした。EMGを計測する筋及びその筋が寄与する関節の動きについては表４を参照することができる。表３と比べて各グループに属する筋本数が若干異なるのは、グルーピングの前提となる筋骨格モデル（simplified model）における筋本数が第１実施形態と第２実施形態とで異なるためである。

本実施形態で用いた筋のグループの詳細について以下に示す。M_highに属する筋のうち、グループを代表し筋電位が計測されている筋は太字で示した。大胸筋や広背筋など筋骨格モデルにおいては1つの筋が複数本のワイヤで表現されている場合もあるが、ここでは複数本のワイヤで表現されている場合も一つの筋として表記した。グループ１〜８は、M₁,...,M_nEMG（表４のグループ１〜８に対応する）に対応する。グループ９〜１８において、M_highは空集合である。尚、以下のグルーピングは一つの態様を例示するに過ぎないものであり、本発明に適用され得るグルーピングは、本明細書に記載されたものに限定されないことは当業者に理解される。

グループ１：肩関節の水平内転
グループ２：肩関節の水平外転
グループ３：肘関節の屈曲
グループ４：肘関節の伸展
グループ５：股関節の伸展と膝関節の屈曲
グループ６：股関節の屈曲と膝関節の伸展
グループ７：足首関節の背屈
グループ８：足首関節の底屈
グループ９：首関節の屈曲
グループ１０：肩関節の下制
グループ１１：肩関節の挙上
グループ１２：手首関節の屈曲
グループ１３：手首関節の伸展
グループ１４：股関節の外旋
グループ１５：足根中足関節の屈曲
グループ１６：体幹の屈曲
グループ１７：体幹の伸展
グループ１８：体幹の挙上

［Ｄ−２］２次計画法による筋張力推定
各筋グループ
に関しては、同じグループ
に属する筋は同じ関節の動きに寄与し異名筋促通により協同筋として働くと考えられるため、M_＊に属する筋は全て同程度の筋活動度
を持つことが期待される。Hill-Stroeve筋モデル（非特許文献１、２）により筋k(∈M_＊)の筋張力f_kは
のように表すことができる。
F_maxkは筋kの最大筋張力、F_l(l_k)は筋kの長さがl_kのときに筋の発揮しうる筋張力の最大張力に対する比、F_v(l(ドット)_k)は筋kの収縮速度がl（ドット）_kのときに発揮しうる筋張力の最大張力に対する比である。

ここで、式(1)中のF_l(l_k)，F_v(l(ドット)_k)，F_maxkは逆運動学計算により求まるため、筋kが発生させる関節トルクτ_Gk´は定数項と変数項を分離して、
とあらわすことができる。
ただし、J_k∈R^1×ndofは、関節角度に対する筋kの筋長のヤコビ行列、H_k∈R^1×ndofは、
である。さらに
を、
とおくと、式(2)より関節トルクτ_G´は、
となる。

以上を踏まえて各グループを代表する筋活動度ベクトルaを二次計画法により求める。二次計画法とはある目的関数、
を最小にするようなxを、xに関する線形な等式拘束条件や不等式拘束条件の下で求める手法である。ただしx，cはn次元ベクトル、Qはn×n行列である。

第１の実施形態では、M_highに属する筋の筋張力は、筋電位から求めた筋活動度をそのまま用いてHill-Stroeve筋モデルにより筋張力を求め、M_lowに属する筋の筋張力は逆運動学と二次計画法により求めている。
第２の実施形態では、EMGにより計算された筋活動度を参照値として与え、M_highに属する筋を含むすべての筋の筋張力を二次計画法により求める。
第２の実施形態では、最適化計算の計算量の減少は主に筋のグルーピングに依存する。また、ＩＫ、ＩＤ・筋張力推定、描画計算を並列処理することで計算の高速化を図ることができる。これらの各処理にそれぞれスレッドを割り当てることで全体として並列処理を行う（図５Ｂ参照）。複数の時刻のデータに対して、ＩＫ、ＩＤ・筋張力推定、描画計算が同時に行われるので、システムのスループットが向上する。並列計算により描画計算が律速されているため、二次計画法に少々時間がかかるようになっても出力画像のフレームレートが変わらない。

以下の4つの条件を考慮して目的関数及び不等式拘束条件を決定する。
条件１：逆運動学により求められた関節トルクと筋張力により発生する筋張力が等しい。つまり式(7)を満たす。
条件２：筋活動度の総和が最小となるような筋活動度のパターンが発生している。
条件３：筋電位が計測されている筋の筋活動度は筋電位から計算された筋活動度と等しい。
条件４：筋活動度a_＊は0から1の間である。

条件１から、逆運動学によって求まった関節トルクτ_G´と筋活動度より計算される関節トルクH^Taが一致して欲しいので、これらの差H^Ta−τ_G´の大きさの2乗を考える。これにより条件１を満たすような解を得るための目的関数
が得られる。
ただし、W_dyn∈R^ndofは重み行列である。

条件２より筋活動度の総和を最小とするために次の目的関数を考える。
ただし、W_tot∈R^ngは重み行列である。

条件３より筋活動度がHill-Stroeveモデルにより計算されている筋のグループに関しては二次計画法により計算される筋活動度と計測された筋電位より計算される筋活動度が一致すべきであるので、筋活動度がHill-Stroeveモデルにより計算されている筋のグループM_highのみの筋活動度を考える。
a_iEMG(i=1,…,n_EMG)を筋電位により計算されるM_ihighの筋活動度であるとし、a_EMGを以下のように定義する。
以上より、条件３の目的関数を重み行列W_EMG∈R^ngを用いて
とする。

以上の条件１、条件２、条件３より目的関数Zを、
とする。
ただし、k_tot、k_EMG(<0)は重み係数である。
式(13)は式(8)において示した二次計画法の目的関数の形となっており、
と考えることができる。

次に条件４より不等式拘束条件は、
となる。以上の目的関数Z及び不等式拘束条件(17)の下、二次計画法によりaを求めれば、筋k(∈M_＊)の筋張力は式(1)により得ることができる。

次に、目的関数中に含まれる重み行列について述べる。
W_dyn，W_tot，W_EMGは対角行列であるとし、それぞれ
のように表される。
なお、σ_idyn ²(i =1,...,n_dof )，σ_jtot ²(j =1,...,n_g)，σ_kEMG ²(k =1,...,n_EMG)
は、それぞれ、W_dyn=E_ndof×ndof，W_tot =E_ng×ng，W _EMG=E_nEMG×nEMGとおいた時に求まる(τ_G´−H^Ta)，a，(a_high−a_EMG)の各要素の分散である。これにより目的関数Zは無次元化される。
さらにユーザが必要に応じて重み係数k_tot，k_EMGを決定するものとする。

最後に、図１０、図１１に基づいて第１の実施形態の手法、第２の実施形態の手法を比較整理して説明する。図１０は、第１の実施形態の手法のグルーピングを示し、図１１は、第２の実施形態の手法のグルーピングを示す。図１０、図１１において、簡略のため、筋電計の数を３個とする。図１１において破線は空集合である。

図１０において、被験者の複数の筋の一部は、筋グループＭ_ｉ（ｉ＝１，２，３）に分類されている。筋グループＭ_ｉは、筋電計の付いている１つの筋M_iEMG(i=1〜３)と、M_iEMGとおなじ働きをする１つあるいは複数の協同筋をM_ihighとする。M_ihighの筋の活動度は一様にM_iEMGと同じとした上で、それ以外の残りの全ての筋(M_ilowとM_others)をまとめる。この残りの全ての筋を動的な力のつり合いから、SR-Inverse（計算が安定な逆行列のようなもの）を使って計算する。この残りの全ての筋（第１の実施形態では約３５０本）の各一本を変数として計算しており、二次計画法は使ってない。従前の最適化計算では、約１０００本の筋活動を計算する場合に二次計画法や線形計画法を使っていた。

図１１において、被験者の複数の全ての筋は、筋グループＭ_ｉ（ｉ＝１，２，３，４，５）に分類されている。筋グループＭ_ｉは、さらに、活動を同じくする協同筋群毎（起始停止する骨が共通する筋群）にまとめることで複数のサブグループに分類される。筋電計の付いている筋を含むサブグループ（M_ihighから構成される）は３つあり、これらのサブグループには、筋電計の付いている１つの筋とおなじ働きをする１つあるいは複数の協同筋も含まれる。すなわち、図１１におけるM_ihighは、図１０におけるM_iEMG＋M_ihighに相当する。
各筋グループＭ_ｉにおいて、筋M_ihighは以外の筋は全てM_ilowとされ、M_ilowを協同筋群毎にまとめることでM_ilowはさらにサブグループ（M_i,jlowから構成される）に分類される。M_i,jlowはM_ilowのなかで起始停止する骨が共通する筋群である。
第２の実施形態において、サブグループの総数は７０個程度になる。この各サブグループに属する筋の活動は一様として７０個程度の変数で表し、これを二次計画法の最適化計算を行う。このとき筋電が張られた筋M_iEMGの筋活動度は M_ihighの筋の参考値として使用される。

本発明は、スポーツトレーニング、リハビリテーション、医療診断、健康管理、エンターテイメント等の分野において利用することができる。

Claims (18)

1. 被験者の運動時の各関節トルクを、筋骨格モデルを用いて逆動力学計算により算出し、該関節トルクを最適化計算により筋張力へ分配することで各筋の筋張力を推定する方法において、
   被験者の複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類し、
   複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）の各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成し、同じサブグループに属する筋の筋活動度を同じとみなし、
   複数の筋グループＭ_ｉの一部あるいは全部において、前記１つあるいは複数のサブグループのうちの少なくとも１つは、筋電計が装着された１つの代表筋と当該代表筋と起始停止する骨が同じである筋とから形成される第１サブグループであり、
   前記第１サブグループに属する筋の筋張力を、最適化計算を用いずに被験者の運動時に計測された前記代表筋の筋電位から取得し、最適化計算の対象から外すことで、最適化計算における変数を削減し、
   あるいは、
   前記１つあるいは複数のサブグループにおいて、各サブグループを代表する筋活動度を最適化計算で推定することで、最適化計算における変数を削減する、
   筋張力の推定法。
2. 前記第１サブグループに属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得し、
   前記第１サブグループに属しない筋の筋張力を、被験者の運動を実現するのに必要な関節トルクにおいて、前記第１サブグループに属する筋により実現できない関節トルクを実現するように最適化することで推定する、
   請求項１に記載の筋張力の推定法。
3. 各筋グループＭ_ｉにおいて、前記１つあるいは複数のサブグループは、前記第１サブグループと、前記第１サブグループに属しない筋について、起始停止する骨が同じ筋から分類された零個以上のサブグループと、を含み、
   各サブグループを代表する筋活動度を最適化計算で推定する、
   請求項１に記載の筋張力の推定法。
4. 前記第１グループに属する筋については、筋電位から取得した筋活動度と最適化計算により計算される筋活動度が一致すべきであるとして、
   最適化計算において、計測された前記代表筋の筋活動度を参照値として用いる、
   請求項３に記載の筋張力の推定法。
5. 筋張力を、被験者の運動時に実時間で推定する、請求項１乃至４いずれかに記載の筋張力の推定法。
6. 被験者の撮影画像あるいは／および当該撮影画像に基づく合成画像を表示部に表示すると共に、表示された被験者の画像に筋骨格モデルをオーバーレイし、
   請求項１乃至５いずれかに記載の推定法により取得した筋張力に基づく身体内部の活動情報を筋骨格モデルに反映させて視覚的に表示する、
   身体内部の活動情報提示法。
7. 身体内部活動情報を、被験者の運動時に実時間で表示する、請求項６に記載の身体内部の活動情報提示法。
8. 前記身体内部の活動情報は、筋活動である、請求項６、７いずれかに記載の身体内部の活動情報提示法。
9. 筋活動を、筋骨格モデルの筋の色あるいは／および形状の変化によって視覚的に表示する、請求項８に記載の身体内部の活動情報提示法。
10. 前記身体内部の活動情報は、筋活動を、当該筋活動を支配する脊髄神経束の活動として表わしたものである、請求項６乃至９いずれかに記載の身体内部の活動情報提示法。
11. 脊髄神経束の活動は、筋骨格モデル上の各脊髄神経束の位置にシンボルを表示し、シンボルの色あるいは／および形状の変化によって視覚的に表示する、請求項１０に記載の身体内部の活動情報提示法。
12. 被験者の運動時の各関節トルクを、筋骨格モデルを用いて逆動力学計算により算出し、該関節トルクを最適化計算により筋張力へ分配することで各筋の筋張力を推定する方法において、
    被験者の複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類し、
    複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）の各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成し、同じサブグループに属する筋の筋活動度を同じとみなし、
    前記１つあるいは複数のサブグループにおいて、各サブグループを代表する筋活動度を最適化計算で推定することで、最適化計算における変数を削減する、
    筋張力の推定法。
13. 複数の筋グループＭ_ｉの一部あるいは全部において、前記１つあるいは複数のサブグループのうちの少なくとも１つは、筋電計が装着された１つの代表筋と当該代表筋と起始停止する骨が同じである筋とから形成される第１サブグループであり、
    前記第１グループに属する筋については、筋電位から取得した筋活動度と最適化計算により計算される筋活動度が一致すべきであるとして、最適化計算において、計測された前記代表筋の筋活動度を参照値として用いる、
    請求項１２に記載の筋張力の推定法。
14. 被験者の運動時の各関節トルクを、筋骨格モデルを用いて逆動力学計算により算出し、該関節トルクを最適化計算により筋張力へ分配することで各筋の筋張力を推定する筋張力取得手段と、
    被験者の複数の筋を分類するグルーピング手段と、
    を備え、
    前記グルーピング手段は、
    前記被験者の複数の筋を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類する第１グルーピング手段と、
    複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）の各筋グループにおいて、起始停止する骨が同じ筋から１つあるいは複数のサブグループを形成する第２グルーピング手段と、からなり、
    前記筋張力取得手段は、同じサブグループに属する筋の筋活動度を同じとみなし、前記１つあるいは複数のサブグループにおいて、各サブグループを代表する筋活動度を最適化計算で推定する、
    筋張力の推定装置。
15. 被験者の全身あるいは身体の一部の複数の筋の少なくとも一部を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類し、各筋グループＭ_ｉから１つの代表筋を選択して当該代表筋に筋電計を装着し、
    前記複数の筋を、
    各筋グループＭ_ｉの前記代表筋からなる第１筋群Ｍ_ｉＥＭＧと、
    各筋グループＭ_ｉにおいて、第１筋群Ｍ_ｉＥＭＧと起始停止する骨が同じである筋からなる第２筋群Ｍ_{ｉｈｉｇｈ}と、
    前記第１筋群Ｍ_ｉＥＭＧ、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれない筋からなる第３筋群と、
    に分け、
    第１筋群Ｍ_ｉＥＭＧと第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得し、
    前記第３筋群に属する筋の筋張力を、逆動力学計算により、計測した被験者の運動を実現するのに必要な関節トルクを計算し、前記関節トルクにおいて、前記第１筋群Ｍ_ｉＥＭＧ及び前記第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋により実現できない関節トルクを実現するように最適化することで推定する、
    筋張力の推定法。
16. 複数の筋電計と、
    筋電位から筋張力を取得する第１筋張力取得手段と、
    逆動力学計算により関節トルクを計算し、当該関節トルクを実現するように最適化計算を行うことで筋張力を推定する第２筋張力取得手段と、
    被験者の複数の筋を分類するグルーピング手段と、
    を備え、
    前記グルーピング手段は、第１グルーピング手段と第２グルーピング手段とを備え、
    前記第１グルーピング手段は、前記被験者の複数の筋の一部あるいは全部を、筋の運動指向性あるいは異名筋促通に基づいて複数の筋グループＭ_ｉ（ｉ＝１，２，．．．ｎ）に分類するものであり、各筋電計は、各筋グループＭ_ｉから選択された１つの代表筋に装着されており、
    前記第２グルーピング手段は、前記被験者の複数の筋を、
    各筋グループＭ_ｉの前記代表筋からなる第１筋群Ｍ_ｉＥＭＧと、
    各筋グループＭ_ｉにおいて、第１筋群Ｍ_ｉＥＭＧと起始停止する骨が同じである筋からなる第２筋群Ｍ_{ｉｈｉｇｈ}と、
    前記第１筋群Ｍ_ｉＥＭＧ、前記第２筋群Ｍ_{ｉｈｉｇｈ}に含まれない筋からなる第３筋群と、
    に分けるものであり、
    前記第１筋張力取得手段は、第１筋群Ｍ_ｉＥＭＧと第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋の筋張力を、被験者の運動時に計測された前記代表筋の筋電位から取得し、
    前記第２筋張力取得手段は、前記第３筋群に属する筋の筋張力を、逆動力学計算により、計測した被験者の運動を実現するのに必要な関節トルクを計算し、前記関節トルクにおいて、前記第１筋群Ｍ_ｉＥＭＧ及び前記第２筋群Ｍ_{ｉｈｉｇｈ}に属する筋により実現できない関節トルクを実現するように最適化することで推定する、
    筋張力の推定装置。
17. 請求項１４、１６いずれかに記載の筋張力の推定装置において、筋張力を、被験者の運動時に実時間で推定する、リアルタイム筋張力推定装置。
18. 請求項１７に記載のリアルタイム筋張力推定装置と、
    運動時の被験者を撮影する手段と、
    被験者の撮影画像あるいは／および当該撮影画像に基づく合成画像を表示する表示手段と、
    を備え、
    前記表示手段に表示された被験者の画像に筋骨格モデルをオーバーレイし、前記リアルタイム筋張力推定装置により実時間で推定した筋張力に基づく身体内部の活動情報を前記筋骨格モデルに反映させて、視覚的に実時間表示するように構成されている、
    身体内部の活動情報提示装置。