JP2010029340A

JP2010029340A - 筋張力データベースの構築方法、筋張力データベース、筋張力データベースを用いた筋張力計算方法及び装置

Info

Publication number: JP2010029340A
Application number: JP2008192924A
Authority: JP
Inventors: Yoshihiko Nakamura; 仁彦中村; Katsu Yamane; 克山根; Sadahiro Takaya; 定裕高屋; Akihiko Murai; 昭彦村井
Original assignee: University of Tokyo NUC
Current assignee: University of Tokyo NUC
Priority date: 2008-07-27
Filing date: 2008-07-27
Publication date: 2010-02-12
Anticipated expiration: 2028-07-27
Also published as: WO2010013631A1; JP5229796B2

Abstract

【課題】
非侵襲に計測できる運動データのみを用いて、生理的に妥当な筋張力を計算するためのデータベースを構築する。
【解決手段】
筋骨格モデルで表現される被験体が所定の運動をした時にフレーム毎に取得した運動データ及び筋骨格モデルの各筋の筋張力を用いた筋張力データベースの構築方法であって、前記運動データを、フレーム毎に筋骨格モデルの各関節の関節角データで特定し、筋骨格モデルの各関節を駆動する拮抗・協働関係にある筋を関節毎にグループ化して筋グループを形成し、前記筋張力を用いて、各筋グループにおける各筋の筋張力比を求めて、フレーム毎に筋張力比データを生成し、前記関節角データと前記筋張力比データとをフレーム毎に対応付けて格納することで、前記所定の運動についての筋張力データベースを得る。
【選択図】図３

Description

本発明は、筋張力の推定に関するものである。

医療、スポーツの分野において、筋・骨格をモデル化することにより、人体の運動を解析する研究が行われている。医療の分野では、筋骨格モデルを用いたリハビリテーションシステムの研究やリハビリテーションへの応用の研究が行われている。またスポーツの分野では、筋骨格モデルを用いたスポーツパフォーマンス向上のためのトレーニングの評価や人体の運動機能の向上の研究を行っている。

これらで用いられる筋骨格モデルの研究として、Delpらの筋骨格モデルによる逆動力学計算とシミュレーションの研究（非特許文献３）、Rasmussenらによる運動時の筋張力推定に関する研究がある（非特許文献４）。Nakamuraらは運動時に生じる筋・腱張力などの体性感覚情報のより精密な解析を行うために、全身詳細筋骨格モデルを構築し、体性感覚情報の高速動力学計算を行っている（非特許文献１）。ここでは、モーションキャプチャシステムにより計測された運動に逆運動学、逆動力学計算を用いて全身を駆動する筋の発生状態を推定している。この全身詳細筋骨格モデルは、骨格モデルと筋骨格ネットワークから構築され、骨格モデルを剛体リンク、筋、腱、靭帯を能動的もしくは受動的に張力を発生するワイヤとしてモデル化している。

ここで、運動データから筋張力を推定する際には、上述のような筋骨格モデルでは骨格モデルの自由度に対しそれを制御する筋の数が多いため、運動データから一意に筋張力が定まらないという冗長問題が生じる。これは、人体の関節が拮抗筋、協働筋といった複数の筋で駆動されるためである。この解決方法として、目的関数（評価関数）を設定し、最適化によって筋骨格系の制御を行う様々な研究がなされている。数理的な手法として、最大張力を最小化する最適化（非特許文献４、非特許文献５、非特許文献８）、筋活動度の二乗和、三乗和を最小化する最適化が提案されている（非特許文献６）。実験に基づき計測された表面筋電位、筋モデルにより筋張力を推定する研究も行われている（非特許文献７）。

このように、人間の運動における筋張力を推定しようとする場合、最適化計算による方法が用いられてきた。しかし、これらの方法では生理的に妥当な筋張力となる保証がない。ここで、「生理的に妥当」とは、実際に人間が発生している筋張力であることを表す。

Yamaneらは、モーションキャプチャシステムによる運動計測に加え、実際の運動時に発生する表面筋電位、床反力を計測し、筋張力の最適化に考慮することで実際の人体のもつ体性感覚に近い筋張力の計算を行っている（特許文献１、非特許文献２）。Yamaneらの表面筋電位計を用いた最適化はこの冗長問題を解決し、さらに実際の人体の持つ体性感覚情報に近い筋張力の計算が可能である。しかし、筋電データを用いる方法は、生理的に妥当な結果を得ることができるが、多数の電極を被験者に貼る必要があるため、計測に手間と時間を要するという不具合がある。
国際公開番号ＷＯ２００５／１２２９００ Y.Nakamura,K.Yamane, Y.Fujita, and I.Suzuki. So-matosensory computation for man-machineinterface from motion capture data and muscloskeltal human model. IEEETransactions on Robotics, Vol. 21, No. 1, 2004. K. Yamane, Y.Fujita, and Y. Nakamura. Estimation of physically and physiologically validsomatosensory information. In Proceedings of IEEE International Conference onRobotics and Automation, pp. 2635.2641, Barcelona,Spain, April 2005. S.L. Delp, F.C. Anderson, A.S. Arnold, P. Loan, A. Habib,C.T. John, E. Guendelman, and D.G. Thelen. OpenSim: Open-source software tocreate and analyze dynamic simulations of movement. IEEE Transactions onBiomedical Engineering, Vol. 54, No. 11, pp. 1940.1950, 2007. J. Rasmussen, M. Damsgaard, and M. Voigt. Musclerecruitment by the min/max criterion.a comparative study. Journal ofBiomechanics, Vol. 34, No. 3, pp. 409. 415, 2001. J. Rasmussen, M. Damsgaard, E. Surma, S.T. Christensen, M.de Zee, and V. Vondrak. AnyBody.a software system for ergonomic optimization.In Fifth World Congress on Structural and Multidisciplinary Optimization, 2003. R.D.Crowninshield and R.A. Brand. A physiologically based criterion of muscle forceprediction in locomotion. Journal of Biomechanics, Vol. 14, No. 11, pp.793.801, 1981. D.G. Lloydand T.F. Besier. An emg-driven musculoskeletal model to estimate muscle forcesand knee joint moments in vivo. Journal of Biomechanics, Vol. 36, No. 6, pp.765.776, 2003. 藤田，中村，山根，鈴木: "筋骨格人体モデルにおける筋張力計算の数理計画問題," 日本機械学会ロボティクス・メカトロニクス講演会'03講演論文集, 2P2-3F-B7, 2003.

本発明は、非侵襲に計測できる運動データのみを用いて、生理的に妥当な筋張力を計算する方法及び装置を提供すること、並びに、そのためのデータベース構築法及びデータベースを提供することを目的としている。

本発明が採用した第１の技術手段は、
筋骨格モデルで表現される被験体が所定の運動をした時にフレーム毎に取得した運動データ及び筋骨格モデルの各筋の筋張力を用いた筋張力データベースの構築方法であって、
前記運動データを、フレーム毎に筋骨格モデルの各関節の関節角データで特定し、
筋骨格モデルの各関節を駆動する拮抗・協働関係にある筋を関節毎にグループ化して筋グループを形成し、前記筋張力を用いて、各筋グループにおける各筋の筋張力比を求めて、フレーム毎に筋張力比データを生成し、
前記関節角データと前記筋張力比データとをフレーム毎に対応付けて格納することで、前記所定の運動についての筋張力データベースを得る、
筋張力データベースの構築方法、である。

１つの態様では、前記関節角データは、各関節の関節角、関節角速度、関節角加速度を含む。また、各関節の関節角、関節角速度のみを関節角データとしたり、あるいは、さらに高次の微分を関節角データに含めてもよい。
１つの態様では、前記各筋の筋張力比は、各筋グループの筋張力の総和で各筋の筋張力を割ったものである。また、筋張力比は、各筋グループの筋張力の総和で各筋の筋張力を割ったものに限定されるものではなく、要は、筋グループを構成する複数の筋に筋張力を分配できるような比であればよい。

１つの態様では、前記格納されている関節角データ、および／あるいは、前記格納されている筋張力比データは、複数の関節角データ、および／あるいは、複数の筋張力データの代表値である。例えば、別々の運動データに基づいて取得された複数の関節角データ及び筋張力データがある場合に、これらのデータを加工してデータベースに格納してもよい。代表値としては、平均、メディアン、トリム平均、重み付け平均等を用いることができる。

１つの態様では、前記関節角データは、所定の運動をしている被験体の運動データをモーションキャプチャで取得し、取得した運動データに基づく逆運動学計算により取得される。運動データを取得するためのモーションキャプチャシステムは一つの好適な例では光学式モーションキャプチャシステムであるが、運動データを取得するために用いるモーションキャプチャ方式は光学式には限定されない。
１つの態様では、前記筋張力は、
所定の運動をしている被験体の所定部位に装着した表面筋電計により筋電位データを取得し、
逆動力学計算により、計測した運動を実現するのに必要な関節トルクを計算し、
筋電位データ、関節トルクを用いて最適化計算を行うことで取得される。
１つの態様では、さらに、反力データを取得し、当該反力データを用いて環境から受ける接触力を最適化することを含み、
前記関節トルクは、一般化力から最適化された接触力を差し引いて決定される。
この場合、１つの態様では、筋張力の計算方法は、取得した反力データを用いて、環境から受ける接触力τ_Ｃを最適化する第1ステップと、取得した運動データ、取得した筋電位データ、最適化された接触力を用いて、筋張力ｆを最適化する第２ステップとを有する。
反力データは、一つの態様では、力センサによって取得される。用いられる反力データが床反力データの場合には、好ましい態様では、該床反力データはフォースプレートによって取得される。また、床反力データは、被験者の足の裏面側に装着した力センサによって取得してもよい。
また、一つの好ましい態様では、該運動データ、該床反力データ、該筋電位データは、同時計測されるものである。

本発明が採用した第２の技術手段は、
運動データと筋張力比データとを対応付ける筋張力データベースであって、
前記運動データは、フレーム毎に、筋骨格モデルの各関節の関節角データで特定されており、
前記筋張力比データは、筋骨格モデルの各関節を駆動する拮抗・協働関係にある筋を関節毎にグループ化して筋グループを形成し、フレーム毎に、各筋グループにおける各筋の筋張力比を計算することで生成され、
前記関節角データと前記筋張力比データとがフレーム毎に対応付けて格納されている、
筋張力データベース、である。

１つの態様では、前記各筋の筋張力比は、各筋グループの筋張力の総和で各筋の筋張力を割ったものである。

１つの態様では、前記格納されている関節角データ、および／あるいは、前記格納されている筋張力比データは、複数の関節角データ、および／あるいは、複数の筋張力データの代表値である。

本発明が採用した第３の技術手段は、
上記筋張力データベースを用いた筋張力取得方法であって、
筋骨格モデルで表現される被験体の運動データから、フレーム毎に筋骨格モデルの各関節の関節角データを求めるステップと、
１フレームないし複数フレームに対応する関節角データを、筋張力データベースに入力するステップと、
入力された関節角データに近い関節角データを探索し、探索された関節角データに対応付けられた１フレームないし複数フレームに対応する筋張力比データを筋張力データベースから出力するステップと、
出力された筋張力比データ、前記運動データを逆動力学計算することで算出した関節トルク、を用いて最適化計算を行なって筋張力を推定する、
筋張力の取得方法、である。

本発明が採用した第４の技術手段は、
上記筋張力データベースを備えた筋張力取得装置であって、
筋骨格モデルで表現される被験体の運動データから、フレーム毎に筋骨格モデルの各関節の関節角データを求める手段と、
１フレームないし複数フレームに対応する関節角データを、筋張力データベースに入力する手段と、
入力された関節角データに近い関節角データを探索し、探索された関節角データに対応付けられた１フレームないし複数フレームに対応する筋張力比データを筋張力データベースから出力する手段と、
出力された筋張力比データ、前記運動データを逆動力学計算することで算出した関節トルク、を用いて最適化計算を行なって筋張力を推定する手段と、
からなる筋張力の取得装置、である。

本発明の特徴は、関節を駆動する拮抗筋・協働筋の関係を、筋張力の比率として捉えることで、一つの関節を駆動する拮抗・協働関係にある筋をグループ化し、この筋グループ内の各筋の筋張力比を運動データに対応させた点にある。したがって、本発明は、筋骨格モデルの逆動力学計算における筋張力の最適化計算で筋張力比データを用いる点に特徴があるのであって、
「筋骨格モデルの逆動力学計算を行うことで筋張力を取得する方法であって、
に、反力データ、運動データ、当該運動データに対応付けられた筋張力比データを与えて、環境から受ける接触力τ_Ｃ、筋張力ｆを最適化することで筋張力を取得することを特徴とする筋張力の取得方法。
ここで、τ_Ｇは一般化力、Ｊは筋・腱・靭帯のヤコビアン、Ｊ_Ｃは接触点のヤコビアンである。」
と規定することもできる。

また、
「次式を用いて、筋骨格モデルの逆動力学計算を行うことで筋張力を取得する方法であって、
反力データを用いて、環境から受ける接触力τ_Ｃを最適化するステップと、
運動データ、当該運動データに対応付けられた筋張力比データ、最適化された接触力を用いて、筋張力ｆを最適化するステップと、
を有することを特徴とする筋張力の取得方法。
ここで、τ_Ｇは一般化力、Ｊは筋・腱・靭帯のヤコビアン、Ｊ_Ｃは接触点のヤコビアンである。」
と規定することもできる。

筋張力比を用いた筋張力の最適化計算（二次計画法）に用いる目的関数の一つは、
であり、ｚが最小となるｆを求めるものである。
ここで、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力ｆとの差；
Ｋ_１：関節トルクに対する重み；
Ｋ_２：筋張力比に対する重み；
である。
二次計画法に用いられる目的関数は、上記の式に限定されるものではなく、この式と実質的に等価な他の式が設定し得ることが当業者に理解される。

筋張力比を用いた筋張力の最適化計算（線形計画法）に用いる目的関数の一つは、
であり、であり、ｚが最小となるｆを求めるものである。
ここで、
であり、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力との差；
である。
線形計画法に用いられる目的関数は、上記の式に限定されるものではなく、この式と実質的に等価な他の式が設定し得ることが当業者に理解される。

本発明では、あらかじめ標準的な運動における人間の筋張力データを蓄積しておくことで、非侵襲な運動計測（モーションキャプチャから得られる運動データ）のみにより実際に人体で生じる筋張力に近い筋張力の計算を行うことができる。したがって、表面筋電位計を用いずに、表面筋電位を用いた時と同等な筋張力分配計算が可能となる。

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

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

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

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

逆動力学では、運動計測によって得られる運動データを元に、その運動を実現する筋・腱・靭帯の張力を求める。逆動力学の計算法の流れは、１．剛体リンク系の逆動力学による関節トルクの計算；２．ワイヤ長さの関節値に対するヤコビアンの計算；３．関節トルクのワイヤ張力への変換、となる。
以下に詳細に説明する。剛体リンク系の逆動力学計算を用いると骨格モデルにおいて運動を実現するのに必要な関節トルクτ_ｇが計算できる。ダランベールの原理と仮想仕事の原理を用いるとτ_ｇと等価な筋、腱、靭帯張力ｆは、関節角θ_ｇに対する筋、腱、靭帯長ｌのヤコビアンＪを用いて、
と表される。

Ｊの計算方法を示す。ワイヤiの長さl_iの関節値に対するヤコビアンＪ_Ｌｉはｌ_ｉの時間微分と、関節速度を次式のように関係付ける行列である。
このＪ_Ｌｉは以下の手順で計算する。ワイヤiがｍ_ｉ個の経由点（始点・終点を含む）からなるものとし、経由点ｊからｊ＋１までの距離をｌ_ｉ，ｊ（ｊ＝０，１，．．．，ｍ_ｉ）、ｌ_ｉ，ｊの関節速度に関するヤコビアンをＪ_Ｌｉ，ｊとすると、Ｊ_ＬｉはＪ_Ｌｉ，ｊの和、すなわち、
と表される。経由点ｊ位置をｐ_ｉ，ｊとすると、
よって、
となる。ここでＪ_ｐｉ，ｊ＝∂ｐ_ｉ，ｊ／∂θ_Ｇ、すなわち∂ｐ_ｉ，ｊのθ_Ｇに関するヤコビアンであるから、「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」の方法から計算できる。これによって得られたＪ_Ｌｉ，ｊの和からＪ_Ｌｉが計算できる。さらに、全ワイヤのＪ_Ｌｉを行方向に並べてJを得る。

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

fを決定するために、何らかの評価関数と拘束条件を設定し、数理計画法による最適化を用いて解決する方法は、例えば、非特許文献１、４、６、８に開示されている。また、特許文献１、非特許文献２では、この最適化において、実際の運動時に計測される表面筋電位を考慮した評価関数を用いて最適化を行う方法が提案されている（図２参照）。筋張力計算に用いられる最適化計算としては幾つもの手法が提案されていることは当業者に理解されることであり、本発明に適用され得る最適化計算は、本明細書に記載されたものに限定されるものではない。

表面筋電位を用いた最適化（特許文献１、非特許文献２）について述べる。
以下の流れで筋張力の最適化を行う
（１）運動を計測する際、被験者に表面筋電位計を装着すると共に、床反力計を設置し、モーションキャプチャシステムによる運動計測と同時に、表面筋電位、床反力を計測する。
（２）表面筋電位から生理学的筋モデルを用いて筋張力を計算する。
（３）（２）で計算された筋張力からの差と式2.1の誤差の和を評価関数として、これを最小化する筋張力の最適化を行う。

表面筋電位から筋張力を計算する方法について述べる。ここではHillとWilkieの筋モデルを定式化したStroeveの筋モデルを使う。筋電位計からの計測値を一定の時間幅で積分したIEMGは運動神経の活動度ｕを表す。ｕと筋活動度ａの関係は次式で表される。
ここでＴは時間遅れを表すパラメータである。また、ａと筋張力ｆの関係は次式で表される。
ここで、Ｆ_ｍａｘは最大筋張力、Ｆ_ｌ(l)とＦ_ｖ(l（ドット）)はそれぞれ正規化された筋力と筋長、筋長の変化率との関係を表す関数である。Ｆ_ｌ(l)を次式のガウス関数で近似する。
またＦ_ｖ(l（ドット）)は次式で近似する。
ここでK_l，V_sh，V_shl，V_mlは定数でStroeveが示した値を用いるか、モーションキャプチャデータに基づいて同定する。

数理計画法（線形計画法及び二次計画法）を用いた筋張力の最適化方法を示す。
先ず線形計画法を用いた最適化方法（非特許文献１参照）を示す。
式2.11及び前小節で求めた計測から計算される筋張力値ｆ^＊、計測値の対応をとる行列Ｋ_Ｆから以下のようにして最適化を行う。
として、
として、Ｚを最小にするΔτ，Δｆ_ｍａｘ，ｆを求める。
ここでａ^Ｔ _τ，ａ^Ｔ _ｍａｘａ^Ｔ _ｆは全要素が正の定ベクトルである。これにより、筋張力を計測値に近づけることができる。また、式2.12より式2.11の誤差も小さくなるので、力学的にも妥当な筋張力が計算できる。

次に、二次計画法による筋張力の最適化方法を示す。不等式拘束条件つきの二次計画法「M. Renouf and P. Alart.
Conjugate gradient type algorithms forfricional multi-contact problems:
Applications to granular materials. Vol. 194, pp. 2019.2041, 2005」に基づき、評価関数Ｚを、
として、Ｚを最小にするｆを求める。
これにより、筋張力を計測値に近づけることができる。また（τ_ｇ−Ｊ^Ｔｆ）も小さくなるので、力学的にも妥当な筋張力が計算できる。

［Ａ−３］床反力の利用
上記の２つの最適化計算では、説明の煩雑さを回避するため、床反力を用いた環境から受ける接触力τ_Ｃの最適化については触れていない。以下に、接触力τ_Ｃの最適化を考慮した関節トルクτ´_Ｇの求め方について述べる。

ヤコビアンJを用いると、一般化力と筋・腱・靭帯の張力ベクトルfと環境から受ける接触力の関係は次式のように表される。
ここで、τ_Ｃｉは接触点ｉにおける接触力、Ｊ_Ｃｉは接触点ｉのヤコビアンである。
接触点のヤコビアンJ_Ｃは次式で定義される。次式において、ｐｃは接触点の位置・姿勢を表すベクトルないしパラメータである。ヤコビアンの計算方法については上述の文献を参照することができる。
式（２）のうち接触力に関する項をまとめると
と書ける。

筋骨格モデルの逆動力学計算は、τ_Ｇから式（2.20）を満たすｆ、τ_Ｃを求める問題となる。既述のように、一般に筋骨格モデルでは筋の数が自由度数に比べて多いため、ｆは一意には決まらない。また、２つ以上のリンクが環境などと接触している場合には不整定問題となる。これらの値を決定するためには、何らかの評価関数によって最適化を行う必要がある。

まず接触力の計測値を考慮した評価関数と拘束条件を示す。前述のように、一般化力τ_Ｇは実際には駆動されていない腰リンクに作用する６軸力が含まれている。また、内力である筋力でこれに相当する力を発生することはできない。そこで、これを接触力で代用する。式（2.20）から腰リンクの６自由度に対応する行を取り出すと次式のようになる。
ここでＥ_ｈｉｐは必要な行を取り出すための行列である。接触しているリンクが1個の場合は一意な解が存在するが、２個以上の場合は不整定問題となる。

本実施形態では、フォースプレートによって計測された床反力データを使うことによってこれを解決する。そのため、次の項を最適化の評価関数に含める。
ここで、τ^＊ _Ｃは計測された接触力、Ｋ_Ｃは各リンクに作用する接触力を計測される接触力の座標系に変換するための行列である。例えば、複数のリンクが1つの力センサ上にある場合は、それらのリンクに作用する床反力の合計が計測されるので、Ｋ_Ｃによって変換する。

接触力は、法線方向の力、圧力中心点、摩擦力などに関する一定の条件を満たす必要がある。このうち、法線方向の力については次の不等式拘束条件を含むことによって対応する。
ここでＥ_ｖｅｒｔはτ_Ｃの法線方向成分を取り出す行列である。その他の条件に関しては，次式を評価関数に含むことで近似的に考慮する。

最適化計算について説明する。以上に示した式を用いて、最適化計算を以下の２段階で行う。
ステップ１：接触力の最適化（２次計画法）
ステップ２：筋張力の最適化

ステップ２については、既に説明したので、ここで、ステップ１について述べる。
まず、第1段階では不等式拘束条件
を満たし、評価関数
を最小化するτ_Ｃを計算する。ここでｗ_Ｈ、ｗ_Ｃは定数の重みである。なお、評価関数の第1項は式（2.21）の誤差を小さくする効果を持つ。ここではｗ_Ｈ＝１×１０^８およびｗ_Ｃ＝１０を用いた。

以上の計算により接触力が求まると、筋張力により発生しなければならない一般化力が次式で計算できる。
後は、式（2.16）、式（2.17）において、τＧは、τ´_Ｇに置き換えて考えればよい。

［Ｂ］運動と筋張力比の相関
筋電情報を使わずに筋張力配分を行うため、本発明では筋間の張力比と運動の関係に注目する。筋張力には伸張反射や拮抗抑制といった神経構造により一意に決まる部分がある。ある主働筋が働くと、それに応じて協働筋にも興奮効果をもつ神経信号が伝わり、協働筋間で筋張力が分配される。また拮抗筋は逆に抑制効果をもつ神経信号が伝わることで、拮抗筋が抑制される。以上から拮抗、協働関係にある筋の筋張力には一定の関係があると予想される。この筋張力の関係を同じ関節を駆動する拮抗、協働関係にある筋同士での筋張力の比率として考える。この筋張力比が一定になる運動の範囲が明らかになれば、その運動につては、運動データのみから冗長性の少ない筋張力を求めることができる。

運動と筋張力比の相関の検証を行なった。具体的には、同一動作において、異なる速度、異なる被験者間での運動と筋張力比の相関について計測に基づき得られた筋張力データを用いて検証する。

異なる被験者及び運動速度での同一動作の計測について説明する。
モーションキャプチャシステム、表面筋電位計を用いて3人の被験者の動作を計測した。計測パターンは、「遅く歩く×3」、「普通に歩く×3」、「速く歩く×3」、「ジョギング×3」である。

表面筋電位は、脚部の運動で主働筋、拮抗筋として働く以下の8つの筋を選んで計測した。

これらの計測データに対し、非特許文献１の二次計画法による最適化を用いて筋張力を計算した。なお、本実験では特に下肢の筋に着目する。

まず各関節の筋張力比を計算する。筋を左右の腰関節(hip)、膝関節(knee)、足関節(foot)を駆動する拮抗筋、協働筋を含むＮ_group(＝６)の個のグループに分類する。これらのグループには各関節を駆動するのに必要な筋全てが含まれる。二関節筋の場合は複数のグループに属する（すなわち、１つの筋が複数の筋グループに属してもよい。）。各筋が属する関節グループを表１、表２に示す。

計算された筋張力から、各グループにおいて各フレームtでの筋張力比を計算する。筋張力比の計算においては、先ず、各関節について、拮抗・協働関係にある筋をグループ化して関節毎に筋グループを形成する。ここでは、股関節、膝関節、足関節（左・右）の６つの関節に対応する６つの筋グループを形成して、筋グループ毎に各フレームｔの筋張力比を計算する。各筋の筋張力比は、各筋グループの筋張力の総和で各筋の筋張力を割ったものである（式（3.1）参照）。
ｄ_i,j[t](ｊ＝１，２，．．．，Ｎ_ｉ)は、グループｉにおける筋ｊの筋張力比；
ｆ_i,j[t]は、グループｉにおける筋ｊの筋張力；
Ｎ_ｉはグループｉに属する筋の数；である。
これから各グループごとにｄ_ｉ[t]のセットを作る。
これらを列方向に並べてＧ_ｉをつくる。
ここで、Ｎ_gruopはグループ数、Ｔは計測を行ったフレーム数である。

運動と筋張力比の相関を検証する際、以下の5種類のデータ間で検証を行った。
1．同一被験者による同一速度での歩行データ
2．同一被験者による異なる速度での歩行データ
3．異なる被験者の同一速度での歩行データ
4．異なる被験者の異なる速度での歩行データ
5．異なる動作（歩行データとジョギングデータ）

運動と筋張力比の相関を見るために各データにおいて動作の１周期分を取り出し、時間方向に正規化する。ここでは、実際の歩行、ジョギング動作を見て右足先部が地面を離れる瞬間から次に右足先が地面を離れるまでを１周期とした。歩行動作の1周期分は5[ms/frame]で200フレーム前後に相当する。速度に関係なく、1周期分の運動が200フレームになるようスケーリングを行なった。

同一動作において筋張力比にどれほどの相関があるかを大腿直筋及び中殿筋を選び、比較を行った。また大腿直筋は二関節筋であるため、データベース構築の際、二関節筋をどうグループ分けするべきかを検証するため、グループ(knee)及び(hip)における筋張力比を挙げた。

これらの筋張力比データを比較したところ、同一被験者間及び、同一速度間での筋張力比に相関があることが確認できた。また二関節筋である大腿直筋がどちらのグループにおいても筋張力比に相関が見られることから、筋のグループ分けの際、二関節筋を両方のグループに分類することができる。

相関を定量的に比較するため、相関係数を用いて比較した。二つの筋張力比データG1、G2に対し、各時刻t、筋j(j=1、2、3…Ni)における筋張力比d1、j[t]、d2、j[t]から各筋における相関係数rjを求める。
全筋においてrjを計算し、それらの平均r（バー）を求める。各関節において、r（バー）を求め比較を行った。結果を表４に示す。

これらを見ると同一速度間、同一被験者間では高い相関があることが確認できる。また速度、被験者ともに異なる場合は相関が低くなる。しかし、同一運動内での相関係数に比べて異なる運動間での相関係数は明らかに低い。これより同じ運動内での筋張力比には相関があることを確認できた。

［Ｃ］筋張力データベースの構築及び筋張力データベースを用いた最適化
上述の通り、同一種類の運動（例えば、歩行）では、異なる被験者間、異なる速度間で筋張力比に相関があることが確認された。従って、表面筋電位等を用いた詳細な筋張力比を、同一種類の新たなデータに対して適用することが可能である。そこで運動の種類ごとに筋張力比データと運動データを対応付けるデータベースを構築する。新たな運動データに対し、データベースから得られた筋張力比を筋張力の最適化に用いることで、運動データのみから実際の人体で生じる筋張力に近い筋張力を求めることが可能になる。以下にその方法を示す。

［Ｃ−１］筋張力データベースの構築
前章で述べた表面筋電位を用いた最適化により筋張力を計算する。
（１）モーションキャプチャ、表面筋電位計、床反力計を同時に用いて運動を計測する。
（２）逆運動学計算により運動データを計算する。
（３）逆動力学計算により、計測した運動を実現するのに必要な関節トルクを計算する。関節トルクを計算する際に、環境から受ける接触力τ_Ｃを考慮する（上述のτ´_Ｇの算出を参照）。
（４）計算された関節トルクに対し、表面筋電位を用いて最適化を行い、筋張力を計算する。最適化については、既述の説明を援用することができる。

各運動データｋに対し、ステップ（２）で得られた時刻tでの全ての関節角θ_ｋ[t]から、関節角速度、関節角加速度を求め、関節角、関節角速度、関節角加速度から関節角データｍ_ｋ[t]をつくる。
ここでΔtは1フレームあたりの計測時間である。
ｍ_ｔ[t]を時系列に並べて関節角データＭ_ｋを作る。(Ｍ_ｋ＝[ｍ_ｋ[t]，ｍ_ｋ[t]，．．．])

拮抗、協働関係にある筋をＮ_group(=6)個のグループに分け、筋張力比データＧ_ｉを式3.1、3.2、3.3より計算する。
データベースには、関節角データＭ_ｋと筋張力比データＧ_ｋｉ(ｉ=1，2，3，．．．，Ｎ_group)の組
を格納する。

すなわち、データベースには下記の対応が保存される。
ここで、ｄ_ｉｔは、フレームｔにおける関節ｉ（ｉ＝１，２，．．，Ｎ）に属する各筋ｊ（ｊ＝１，２，・・・，Ｎ_ｉ）の筋張力比を表している。
フレームｔにおいて、全ての関節には、当該関節を駆動する筋グループの筋張力比のセット（まとめて、ベクトルで表してもよい）が対応している。
一方、運動データは、フレームデータ（フレームｔにおける姿勢データ）を時系列に並べたものであり、フレームｔにおける姿勢データは、各関節の関節角データ（関節角度、関節角速度、関節角加速度）により特定される。
したがって、フレームｔにおいて、各関節の筋張力比データと、各関節の関節角データを対応させることができ、これらを組として保存する。フレームｔにおける関節角データが特定されれば、それと組である筋張力比データも特定される。

一つの態様では、図３、図４に示すように、筋張力データベースでは、運動毎（歩行、ジョギング等）にデータベースが作成されている。全ての関節の関節角データの時系列データＭ_１に対して、全ての関節（各筋グループ）に関する筋張力比データ（Ｇ_１１、Ｇ_１２、・・・）が対応して格納されている。

運動毎のデータベースを作成する場合、データベース構築用のデータを運動毎に分類する作業と、入力された運動データについて運動の種類を特定する作業が必要になる。この分類作業、特定作業は手作業で行うことができるが、手作業を行なわない方法の一つとして、運動データを自律的にセグメンテーション、クラスタリングし、運動の種類を認識する手法を適用することも可能である。このような手法は、例えば、「D. Kulic, W. Takano and Y. Nakamura: "Combining Automated
On-Line Segmentation and Incremental Clustering for Whole Body Motions,"
IEEE International Conference on Robotics and Automation, pp. 2591-2598, 2008.」に開示されている。この方法を用いることで、データベース構築時には自動的に、長い運動が適当な単位で分割され、似た運動がグループ化されるので、運動のグループ毎にデータベースを作成することができる。また、データベースを用いて検索を行なう場合には、観測された運動がどのグループに近いかが判るので、どのデータベースを使えばよいかがわかる。より具体的には、データベース中のある種類の運動を表す統計的モデル(隠れマルコフモデル)から、入力された運動が生成される尤度を計算することで、入力された運動に対応するデータベースを選択することができる。

関節角データと筋張力比データの対応を、さらに、具体的にどのような構造で保存するかについては、様々なやり方が取り得ることが当業者に理解される。例えば、最も単純な構造としては、多数の各フレームの関節角・速度・加速度・筋張力比のセットをフラットに格納しておき、データベースを使う際には関節角・速度・加速度が最も良く一致するフレームの筋張力比を採用することができる。この場合、各フレームの時間間隔が一定である必要はない。

関節角データと筋張力比データの対応を、木構造（典型的には、二分木構造）を備えた階層型データベース構造としてもよい。例えば、木構造の最上位層は、複数の運動データの全てのフレームが含まれている１つのノードであり、運動データの各フレームは、木構造の各階層のいずれかの１つのノードに含まれていると共に、各ノードには、上位層から下位層に行くにしたがってより近い状態量（関節角・速度・加速度）を備えたフレームが含まれている。各ノードに含まれた複数のフレームの状態量に基づく各ノードの代表値（平均、メディアン、トリム平均、重み付け平均等）を求めておく。そして、同一ノードに入るフレームの筋張力比の代表値（平均、メディアン、トリム平均、重み付け平均等）を求めておく。この場合、各フレームの時間間隔は一定であることが望ましい。

［Ｃ−２］筋張力データベースを用いた筋張力計算
前節で構築したデータベースを用いて以下の方法で、筋張力を計算する。全体の流れを図５に示す。
（１）モーションキャプチャによりキャプチャデータを得る。
（２）逆運動学計算により関節角、関節角速度、関節角加速度、を計算する。
（３）逆動力学計算から、計測した運動を実現するのに必要な関節トルクを計算する。
（４）関節角度、関節角速度、関節角加速度、をデータベースに入力する。データベースは入力された運動に最も近いデータを検索し、筋張力比Ｇ_ｉ′を出力する。
（５）データベースから得られた筋張力比Ｇ_ｉ′を用いて筋張力の最適化を行う。

新たな運動データに対し、前節と同様に関節角データＭ′を求める。新たな運動データの運動の種類と同じ種類の運動のＭ_ｋをデータベース内から探し、各時刻tにおいて|ｍ′[ｔ]−ｍ_ｋ[ｓ]|を最小にするｋ、ｓを求める。ここでｍ′[t]はＭ′の時刻tでの関節角データである。ｋは運動の種類に対応し、ｓはフレーム番号に対応する。

ｋ、ｓの探索は、用いられるデータベースの種類によって異なることは当業者に理解される。例えば、データベースが運動の種類毎に作成されている場合であって、かつ、用いるデータベースが予め決定されている場合には、入力された運動データのフレームｔの関節角データとの距離が最小となる関節角データを有するフレームｓを求める。データベースが複数の種類の運動データを含んでいる場合にも、入力された運動データのフレームｔの関節角データとの距離が最小となる関節角データを有するフレームｓを求める。例えば、データベースが運動の種類毎に作成されている場合であって、かつ、用いるデータベースが未知の場合には、データベース中のある種類の運動を表す統計的モデル(隠れマルコフモデル)から、入力された運動が生成される尤度を計算することで、入力された運動に対応する運動ｋを選択することができる。そして、運動ｋに含まれるフレームの中から入力された運動データのフレームｔの関節角データとの距離が最小となる関節角データを有するフレームｓを求める。
求めたｓに対応する筋張力比ｄ_ｉ[s]を時系列に並べて筋張力比Ｇ´_ｉ(ｉ＝1、2、3、．．．Ｎ_group)をつくる。

各フレームにおいてデータベースから得られた各関節の筋張力比ｄ_ｉ∈Ｒ^ＮｉからＤ_ｉ∈Ｒ^{Ｎｉ×Ｎｉ}を作る。Ｎ_ｉはグループｉに属する筋の数である。
ここでＤ_ｉｆ_ｉはグループｉの筋張力ｆ_ｉを筋張力比ｄ_ｉを用いてグループ内で分配した筋張力と筋張力ｆ_ｉとの差を表す。
ここで、Ｄ_ｉ（グループi内のみの差）を全てのグループについてまとめた行列Ｋ_ｆｄは、
となる。
そして、筋張力ｆを各関節の筋張力ｆ_ｉに写像する行列をＫ_ｇｉ∈Ｒ^{Ｎｉ×Ｎelement}とする。Ｎ_elementは全ての筋の数である。
Ｋ_ｇｆによってグループ毎に並べ替えられた全筋張力をＫ_ｆｄにかけると、全筋について理想的な筋張力比からの差Ｋ_ｆｄＫ_ｇｆが得られる。

一般化座標から筋、腱、靭帯長へ写像するヤコビアンＪ及びＫ_ｆｄからから以下のようにして二次計画法を解く。
ここで、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力ｆとの差；
Ｋ_１：関節トルクに対する重み；
Ｋ_２：筋張力比に対する重み；
である。
Ｚが最小になるようなｆを求める。これにより、筋張力比ｄに近づくような筋張力が求められる。

また線形計画法を用いた最適化方法では、
のＺを最小にするｆを求める。

本発明では、非接触・非侵襲な計測が可能な運動データのみにより生理的に妥当な筋張力の推定が可能となるので、筋骨格モデルの応用範囲が広がる。特に、ロボットの人間のコミュニケーションにおいてロボットが人間の内部状態を推定するためには、非接触な計測が必須である（従来の筋電データは、接触式電極が必要であった）。例えば、ロボットに搭載されたカメラ画像のみから計測した運動データを用いて、人間の筋張力を推定することができ、ロボットによる人間の作業支援に役立てることができる。

人の筋骨格モデルを表す図である。筋電位、床反力を用いた筋張力の最適化を示す図である。本発明に係る筋張力データベースを示す図である。図３と類似の図であり、本発明に係る筋張力データベースを示す図である。本発明に係る筋張力データベースを用いた筋張力計算を示す図である。

Claims

筋骨格モデルで表現される被験体が所定の運動をした時にフレーム毎に取得した運動データ及び筋骨格モデルの各筋の筋張力を用いた筋張力データベースの構築方法であって、
前記運動データを、フレーム毎に筋骨格モデルの各関節の関節角データで特定し、
筋骨格モデルの各関節を駆動する拮抗・協働関係にある筋を関節毎にグループ化して筋グループを形成し、前記筋張力を用いて、各筋グループにおける各筋の筋張力比を求めて、フレーム毎に筋張力比データを生成し、
前記関節角データと前記筋張力比データとをフレーム毎に対応付けて格納することで、前記所定の運動についての筋張力データベースを得る、
筋張力データベースの構築方法。
前記関節角データは、各関節の関節角、関節角速度、関節角加速度を含む、請求項１に記載の筋張力データベースの構築方法。
前記各筋の筋張力比は、各筋グループの筋張力の総和で各筋の筋張力を割ったものである、請求項１、２いずれかに記載の筋張力データベースの構築方法。
前記格納されている関節角データ、および／あるいは、前記格納されている筋張力比データは、複数の関節角データ、および／あるいは、複数の筋張力データの代表値である、請求項１乃至３いずれかに記載の筋張力データベースの構築方法。
前記関節角データは、所定の運動をしている被験体の運動データをモーションキャプチャで取得し、取得した運動データに基づく逆運動学計算により取得される、請求項１乃至４いずれかに記載の筋張力データベースの構築方法。
前記筋張力は、
所定の運動をしている被験体の所定部位に装着した表面筋電計により筋電位データを取得し、
逆動力学計算により、計測した運動を実現するのに必要な関節トルクを計算し、
筋電位データ、関節トルクを用いて最適化計算を行うことで取得される、
請求項１乃至５いずれかに記載の筋張力データベースの構築方法。
さらに、反力データを取得し、当該反力データを用いて環境から受ける接触力を最適化することを含み、
前記関節トルクは、一般化力から最適化された接触力を差し引いて決定される、請求項５に記載の筋張力データベースの構築方法。
運動データと筋張力比データとを対応付ける筋張力データベースであって、
前記運動データは、フレーム毎に、筋骨格モデルの各関節の関節角データで特定されており、
前記筋張力比データは、筋骨格モデルの各関節を駆動する拮抗・協働関係にある筋を関節毎にグループ化して筋グループを形成し、フレーム毎に、各筋グループにおける各筋の筋張力比を計算することで生成され、
前記関節角データと前記筋張力比データとがフレーム毎に対応付けて格納されている、
筋張力データベース。
前記関節角データは、各関節の関節角、関節角速度、関節角加速度を含む、請求項８に記載の筋張力データベース。
前記各筋の筋張力比は、各筋グループの筋張力の総和で各筋の筋張力を割ったものである、請求項８、９いずれかに記載の筋張力データベース。
前記格納されている関節角データ、および／あるいは、前記格納されている筋張力比データは、複数の関節角データ、および／あるいは、複数の筋張力データの代表値である、請求項８乃至１０いずれかに記載の筋張力データベース。
請求項８乃至１１いずれかに記載の筋張力データベースを用いた筋張力取得方法であって、
筋骨格モデルで表現される被験体の運動データから、フレーム毎に筋骨格モデルの各関節の関節角データを求めるステップと、
１フレームないし複数フレームに対応する関節角データを、筋張力データベースに入力するステップと、
入力された関節角データに近い関節角データを探索し、探索された関節角データに対応付けられた１フレームないし複数フレームに対応する筋張力比データを筋張力データベースから出力するステップと、
出力された筋張力比データ、前記運動データを逆動力学計算することで算出した関節トルク、を用いて最適化計算を行なって筋張力を推定する、
筋張力の取得方法。
前記関節角データは、各関節の関節角、関節角速度、関節角加速度を含む、請求項１２に記載の筋張力の取得方法。
さらに、反力データを取得し、当該反力データを用いて環境から受ける接触力を最適化することを含み、
前記関節トルクは、一般化力から最適化された接触力を差し引いて決定される、請求項１２、１３いずれかに記載の筋張力の取得方法。
最適化計算の目的関数は、
であり、ｚが最小となるｆを求めるものである、請求項１２乃至１４いずれかに記載の筋張力の取得方法。
ここで、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力ｆとの差；
Ｋ_１：関節トルクに対する重み；
Ｋ_２：筋張力比に対する重み；
である。
最適化計算の目的関数は、
であり、であり、ｚが最小となるｆを求めるものである、請求項１２乃至１４いずれかに記載の筋張力の取得方法。
ここで、
であり、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力ｆとの差；
である。
請求項８乃至１１いずれかに記載の筋張力データベースを備えた筋張力取得装置であって、
筋骨格モデルで表現される被験体の運動データから、フレーム毎に筋骨格モデルの各関節の関節角データを求める手段と、
１フレームないし複数フレームに対応する関節角データを、筋張力データベースに入力する手段と、
入力された関節角データに近い関節角データを探索し、探索された関節角データに対応付けられた１フレームないし複数フレームに対応する筋張力比データを筋張力データベースから出力する手段と、
出力された筋張力比データ、前記運動データを逆動力学計算することで算出した関節トルク、を用いて最適化計算を行なって筋張力を推定する手段と、
からなる筋張力の取得装置。
前記関節角データは、各関節の関節角、関節角速度、関節角加速度を含む、請求項１６に記載の筋張力の取得装置。
さらに、反力データを取得し、当該反力データを用いて環境から受ける接触力を最適化する手段を含み、
前記関節トルクは、一般化力から最適化された接触力を差し引いて決定される、請求項１６、１７いずれかに記載の筋張力の取得装置。
最適化計算の目的関数は、
であり、ｚが最小となるｆを求めるものである、請求項１７乃至１９いずれかに記載の筋張力の取得装置。
ここで、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力ｆとの差；
Ｋ_１：関節トルクに対する重み；
Ｋ_２：筋張力比に対する重み；
である。
最適化計算の目的関数は、
であり、であり、ｚが最小となるｆを求めるものである、請求項１７乃至１９いずれかに記載の筋張力の取得装置。
ここで、
であり、
τ_ｇ：関節トルク；
Ｊ：関節角から筋、腱、靭帯長へのヤコビアン；
ｆ：筋張力；
Ｋ_ｇ：筋張力ｆを各関節（各筋グループ）の筋張力へ写像する行列；
Ｋ_ｆｄｆ：筋張力ｆを、筋張力比を用いて筋グループ内で分配した筋張力と、筋張力ｆとの差；
である。