JP2015112453A

JP2015112453A - 運動解析装置、運動解析方法及び運動解析プログラム

Info

Publication number: JP2015112453A
Application number: JP2013259093A
Authority: JP
Inventors: 文夫宮崎; Fumio Miyazaki; 平井　宏明; Hiroaki Hirai; 宏明平井; 充典植村; Mitsunori Uemura; かんな宇野; Kanna Uno; 貴紀奥; Takanori Oku
Original assignee: Osaka University NUC
Current assignee: Osaka University NUC
Priority date: 2013-12-16
Filing date: 2013-12-16
Publication date: 2015-06-22
Anticipated expiration: 2033-12-16
Also published as: JP6518932B2; US10631751B2; US20160324436A1; WO2015093224A1

Abstract

【課題】所望の運動を実現するために中枢神経系が選択する運動指令を、筋拮抗比、筋拮抗和の概念の下、筋骨格モデルに基づいて３つの特徴量で推定する装置を提供する。。【解決手段】運動解析装置は、運動者の筋電位を測定する筋電位測定部２０と身体運動を測定する運動測定部３０と、測定部２０，３０で得られた測定情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を、筋骨格モデルに基づき算出する、剛性楕円算出部１０８、平衡点算出部１０７、及び平衡点を記述する基底ベクトルである筋シナジー算出部１０６とを備えた。【選択図】図４

Description

本発明は、身体適所、例えば肢の先端等の動作点の運動制御に関わる特徴量を算出する運動解析技術に関する。

例えば運動技能の習得において効果的な練習を行うには、学習者の習熟度を適切に判断し、質の高い情報を提供することが望ましい。運動実現に直接関与する動作点には運動評価及び運動学習支援のための豊富な情報が含まれており、従来、運動対象のインピーダンス特性や平衡点軌道についての研究がなされている。現在、これらの特徴量を推定する方法として、対象運動に摂動を与え、その変動を身体モデルに基づき解析する方法（非特許文献１，２）、対象運動時の筋群活動、運動学及び手先力を計測し、これらの関係を記述するモデルを同定することで上記２つの特徴量を推定する方法（特許文献１）が提案されている。また、動作点のインピーダンス、平衡点の実現に関わる筋群活動を統計手法を用いて解析し、運動実現のための機能単位（筋シナジー）に分解する方法（特許文献２）も提案されている。

特開２００４−７３３８６号公報 WO２０１１／０３０７８１号公報

Hiroaki Gomi and Mitsuo Kawato, ``Equilibrium-Point Control Hypothesis Examined byMeasured Arm Stiffness During Multijoint Movement," Science, vol. 272 , no. 5258, pp. 117-120, 1996. 五味裕章、川人光男、``水平面内における多関節運動中の人腕機械インピーダンスの計測" 計測自動制御学会論文集、vol. 32、no. 3、pp. 369-378、1996.

非特許文献１，２には、動作点のインピーダンスと平衡点を推定する方法が記載されているが、かかる方法では、高精度なロボットと運動計測装置が必要となる。従って、装置が大掛かりになったり、適用可能な運動タスクが限られたり、また摂動が反射系に影響を与えるといった技術的課題が存在する。また、筋シナジーに関する記載はない。

特許文献１では、モデル作成のために膨大な数のデータ収集と解析が必要となり、さらに神経回路モデルを用いて推定する場合、その物理的意味が不明瞭であるといった問題がある。また、筋シナジーに関しても記載されていない。一方、特許文献２には、インピーダンス、平衡点及び筋シナジーに関する記載はあるものの、抽出される筋シナジーの物理的意味に関しては統計データに基づく示唆に留まっている。また、剛性楕円に関する記載はない。

本発明は、上記に鑑みてなされたもので、身体運動を筋骨格構造の機械系と捉え、筋拮抗比及び筋拮抗和の概念を導入することによって、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトル（筋シナジー）の３つの特徴量を、より簡単な計算によって解析的に求め得る運動解析装置、運動解析方法及び運動解析プログラムを提供するものである。

本発明に係る運動解析装置は、運動者の筋活動及び身体運動を周期的に測定する測定部と、前記測定部で得られた測定情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を筋骨格モデルに基づき算出する演算手段とを備えたものである。

また、本発明に係る運動解析方法は、測定部によって運動者の筋活動及び身体運動を周期的に測定し、前記測定部で得られた測定情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を筋骨格モデルに基づき算出するものである。

また、本発明に係る運動解析プログラムは、前記測定部で周期的に測定された運動者の筋活動及び身体運動の情報を取り込む測定処理手段、前記測定処理手段で取得された情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を筋骨格モデルに基づき算出する演算手段、としてコンピュータを機能させるものである。

これらの発明によれば、測定部によって運動者の筋活動及び身体運動が周期的に測定され、得られた測定情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの３つの特徴量が筋骨格モデルに基づいて算出される。３つの特徴量の算出は、筋活動、身体運動及び筋骨格モデルの情報を用いてリアルタイムで代数計算処理が可能である。算出結果は、一旦保管する態様としてもよいし、必要に応じてリアルタイムで取り出されて利用する態様とすることも可能である。

従って、特徴量の推定に際しては、物理モデルに基づき解析的に推定が行え、また煩雑なパラメータ同定を必要とせず、大掛かりな装置を要しない代数計算によって推定でき、さらに筋拮抗比、筋拮抗和の概念の下、前記特徴量を統一的に扱うことが可能となる。

また、推定される筋シナジーを含む３つの特徴量は、物理モデルに基づき推定される量であるため、神経回路モデルや統計手法から推定されるものに比べ、その物理的意味が明確化される。

また、本発明は、前記測定部は、前記筋活動として筋電位を測定する筋電位検出部と、前記身体運動として関節の動きを測定する運動検出部とを備え、前記演算手段は、筋拮抗比及び筋拮抗和を算出する第１の算出手段と、前記筋骨格モデルに基づいて身体運動の測定部位を前記動作点に写像するヤコビ行列を算出する第２の算出手段とを備え、算出された筋拮抗比、筋拮抗和及びヤコビ行列を用いて前記特徴量を算出することを特徴とする。この構成によれば、剛性楕円及び平衡点の推定に際し、摂動法では技術的に推定が困難もしくは不可能な複雑な運動タスクに対しても適用可能である。また、摂動に伴う反射系の影響もない。さらに、剛性楕円及び平衡点の推定に際し、神経回路モデル等を用いて筋電図を解析する方法のように膨大な数量のデータ収集や解析が不要である。

また、本発明は、算出結果を出力する出力部を備えたものである。この構成によれば、所望の運動を実現するために中枢神経系が選択する運動解（運動指令）を３つの特徴量に基づき、モニタ等の出力部で時々刻々と観測することができる。

また、本発明は、前記身体運動を観察する撮像部を備え、前記出力部は、画像を表示する表示部を備え、前記表示部に、前記特徴量のうち、少なくとも前記平衡点及び前記剛性楕円の表示と、前記撮像部で撮像された身体運動の画像とを表示することを特徴とする。この構成によれば、算出結果が運動中の身体画像と対応付けて表示されることで、特徴量に対する視認性が高まる。

また、本発明は、前記出力部は、前記平衡点を前記身体運動の画像に対応させて表示することを特徴とする。この構成によれば、平衡点の位置を適正に視認することができる。

また、本発明は、前記出力部は、前記剛性楕円を前記身体運動の画像の動作点に対応させて表示することを特徴とする。この構成によれば、剛性楕円の状態を適正に視認することができる。

また、本発明は、前記出力部は、前記特徴量に基づき支援対象に物理的支援を施す駆動部であることを特徴とする。この構成によれば、リハビリテーション分野やインタフェース評価装置への適用が可能となる。

本発明によれば、所望の運動を実現するために中枢神経系が選択する運動指令を、筋拮抗比、筋拮抗和の概念の下、筋骨格モデルに基づいて３つの特徴量で推定することができる。

本発明に係る運動解析方法の概要及び適用態様を説明する図である。本発明の運動解析方法の一例を説明する図で、（Ａ）は上肢による身体運動の一例を示す図、（Ｂ）は上肢の筋骨格モデルである筋と座標系との定義を示す図である。筋拮抗比と筋拮抗和及びその機能を示す図である。本発明に係る運動解析装置の一実施形態を示すブロック構成図である。平衡点の算出方法を幾何学的に説明する図である。剛性楕円の算出方法を幾何学的に説明する図である。運動解析処理の手順の一例を示すフローチャートである。表示処理の手順の一例を示すフローチャートである。被験者Ａの実験から得られた筋拮抗和と筋拮抗比の推移を示す図で、（Ａ）は筋拮抗和を示し、（Ｂ）は筋拮抗比を示す。動径、偏角、Ｎｕｌｌ方向のシナジーベクトルと、それぞれのシナジースコアを示す特性図で、（Ａ）は被験者Ａのデータであり、（Ｂ）は被験者Ｂのデータである。摂動法により算出された関節剛性を示す特性図で、（Ａ）は被験者Ａのデータであり、（Ｂ）は被験者Ｂのデータである。摂動法より求めた関節剛性を手先剛性に変換したものと、筋拮抗和およびＮｕｌｌ方向のシナジースコアから推定した手先剛性を示す特性図で、（Ａ）は摂動法より求めたデータであり、（Ｂ）は今回の推定法から求めたデータであり、（Ｃ）はＮｕｌｌ方向のシナジースコアを足し合わせずに求めたデータである。手先剛性より算出された手先の平衡点軌道と、筋シナジースコアから推定した手先の平衡点軌道、およびそれらの接線速度の大きさを示す図で、（Ａ）は摂動法より求めたデータであり、（Ｂ）は今回の推定法から求めたデータである。

図１は、本発明に係る運動解析方法の概要及び適用態様を説明する図である。図２は、本発明の運動解析方法の一例を説明する図で、（Ａ）は上肢による身体運動の一例を示す図、（Ｂ）は上肢の筋骨格モデルである筋と座標系との定義を示す図である。図３は、筋拮抗比と筋拮抗和及びその機能を示す図である。図４は、本発明に係る運動解析装置の一実施形態を示すブロック構成図である。

図１において、運動者１は中枢神経系からの運動指令に従って運動を行う。ここでは、運動者１の運動は筋骨格構造の機械系（筋骨格モデル）として扱われる。筋骨格モデルの情報は、好ましくは装置に予め記憶されている。運動者１の運動は、筋電位及び身体運動（運動学）として測定される。その結果、運動者１の情報として、筋骨格モデル、筋電位及び運動学が得られる。これらの情報を用いて３つの特徴量が推定演算され、モニタ１５２に表示（可視化）される。

３つの特徴量とは、剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーである。運動者１の運動状況が撮影されており、その撮影映像がモニタ１５２に表示されている。剛性楕円は、動作点を中心にして（対応して）表示され、平衡点は、撮影映像の座標系に対応して例えば点状の所定のマークで表示されている。剛性楕円は、動作点の動き難さ（易さ）を規定するもので、中心からの方向と径の長さとによって表される。平衡点は、動作点に対する運動指令を示す位置を表している。筋シナジーは、本実施形態では、偏角方向、動径方向、及びそれらに直交するＮｕｌｌ方向（零方向）の各成分で記述される。

モニタ１５２に３つの特徴量を表示することで、運動者１の運動状況に対する運動評価に供することが可能となる。また、運動者１の運動状況をより高めるべく運動学習支援用に供することも可能となる。

図２（Ａ）において、運動者１（被験者）は、作業机２上で右肩がマニピュレータ３の正面に来るように着席し、椅子の背もたれに両肩を固定されている。被験筋は、本例では、水平面内の上肢運動で主要な役割を担う３対６筋（うち１対は二関節筋）であり、被験筋と座標系の定義を図２（Ｂ）に示す。前腕と上腕は、天井から吊り下げ部材６によって吊り下げられ、あるいは机上でスライド自在なように下方から支持されることで、重力の影響を取り除いている。また、本例では、プラスチックのギブスによって右腕の手首関節は固定され、下方に向けられたレーザーポインタ４とマニピュレータ３の先端が手先（動作点）に固定されている。作業机２の上には、運動方向と距離を示す紙５を置き、レーザーポインタ４で被験者が手先位置をリアルタイムに視認できるようにした。上部には動画を撮影する撮像部としてのカメラ４０が備えられ、少なくとも運動部位が撮影される。モーションキャプチャー等の光学トラッキングシステムのためのマーカ８が、所定位置、例えば左肩、右肩、右肘の関節中心と、右手先のレーザーポインタ４の直上に固定される。各マーカ８は、後述するように運動測定部３０（図４参照）で検知される。

図２（Ｂ）は、右肩を基準にしたｘ−ｙ座標系を示し、ρは右肩からの距離、φは角度である。また、運動学データであるθ１，θ２は、光学トラッキングシステムで計測されるもので、肩関節の関節角、及び肘関節の関節角を表わしている。

図２（Ｂ）に示すように、上肢の各筋は、三角筋後部繊維ｍ１、三角筋前部繊維ｍ２、上腕三頭筋長頭ｍ３、上腕二頭筋ｍ４、上腕三頭筋外側頭ｍ５、腕橈骨筋ｍ６であり、ｍ１とｍ２，ｍ３とｍ４、ｍ５とｍ６とは拮抗筋対である。また、図３に示すように、ｍｉ（ｉ＝１〜６）を筋電位として表すものとしたとき、拮抗筋比、拮抗筋和は、数１中の、式（１）、式（２）に示すようにして算出される。

なお、筋電位は、ＥＭＧ（筋活動を直接計測する筋電図：electromyogram）計測で得られ、筋協調の最小単位と考える拮抗する筋同士の活動の比（筋拮抗比）の概念と、拮抗する筋同士の活動の和（筋拮抗和）の概念とを導入している。また、筋電位は、筋間での正規化をとるために最大等尺性収縮時の筋電位に対する百分率％ＭＶＣで表している。

図４において、運動解析装置は、運動者１の筋電位を検出する筋電位測定部２０と、運動者１に取り付けたマーカ８を検出して運動学データである関節角θ１，θ２を測定する運動測定部３０と、測定した筋電位、関節角に基づいて所定の信号処理を施すマイクロコンピュータを代表とするコンピュータ等で構成される信号処理部１０とを備える。

筋電位測定部２０は、所要数の電極２１１と、各電極２１１に生じる電気信号を検出する筋電位検出回路２１と、検出された電気信号に対して所定の前処理を施す前処理回路２２とを備える。電極２１１は、本例では、筋ｍ１〜ｍ６の筋電位を検出するためのものである。体表で得られる筋電位信号は、レベルが数十μＶ〜数百μＶ、周波数が５Hz〜５００Hz程度の交流信号である。そこで、前処理回路２２は、筋電位を処理可能なレベル（数千倍）まで増幅するアンプ、筋電位の主要周波数帯の信号のみを通過させるバンドパスフィルタ、及び全波整流回路を備えている。また、前処理回路２２は、筋電位信号をデジタルで処理するべく、出力側にＡＤ変換部を備えている。

運動測定部３０は、前述のモーションキャプチャー等の光学トラッキングシステムのための、例えば特定色を発する発光素子あるいは特定色に反応する表面処理されたマーカ８を検出する所要数の撮像素子３１１と、各撮像素子３１１からの撮像信号を検出する位置検出回路３１と、検出された撮像信号に対して所定の前処理を施す前処理回路３２とを備える。撮像素子３１１は所要領域においてマーカ８の位置変化を検出する。運動測定部３０は、光学的な位置測定の他、磁気発生器と磁気センサとから構成される３次元上の位置と向きが検出可能な公知の測定器を採用してもよい。

信号処理部１０は、ＣＰＵを有するマイクロコンピュータで構成されており、本発明に係る信号処理のための処理プログラムを記憶するＲＯＭ１１、処理途中のデータを一時的に記憶するＲＡＭ１２、所要の指示を行うためのテンキーやマウス等で構成される操作部１３、及び人体筋骨格モデルのデータその他の所要のデータを記憶する外部記憶部１４が接続されている。なお、人体筋骨格モデルのデータは、ＲＯＭ１１又はＲＡＭ１２内に格納される態様でもよい。また、信号処理部１０には表示部１５が接続されている。表示部１５は、表示処理部１５１とモニタ１５２とを備える。モニタ１５２は、ここでは操作部１３からの入力情報の確認や処理結果を表示すると共に、カメラ４０による撮影映像を表示するものである。表示処理部１５１は、信号処理部１０での処理結果とカメラ４０からの映像とを対応付けて表示する処理を行う。なお、表示処理部１５１は信号処理部１０における表示処理プログラムで実行される機能部とハードウエアの構成部とを備える態様でもよい。

信号処理部１０のＣＰＵは、ＲＯＭ１１からＲＡＭ１２に読み出した処理プログラムを実行することによって、筋電位の測定処理を周期的に実行する筋電位測定処理部１０１、マーカ８の測定処理を周期的に実行する運動測定処理部１０２、筋電位測定処理部１０１及び運動測定処理部１０２から得られた信号から前準備処理を実行する、筋拮抗和算出部１０３、筋拮抗比算出部１０４及びヤコビ行列算出部１０５として機能する。また、信号処理部１０のＣＰＵは、前準備処理の結果を利用して、筋シナジー算出部１０６、平衡点算出部１０７及び剛性楕円算出部１０８として機能する。

筋電位測定処理部１０１は、後述するように%ＭＶＣ化されたＥＭＧを算出する。運動測定処理部１０２は、マーカ８の検出位置からマーカ間の距離や前記した関節角を算出する。筋拮抗和算出部１０３は、前記数１の式（１）を実行し、筋拮抗比算出部１０４は、前記数１の式（２）を実行する。また、ヤコビ行列算出部１０５は、前記数１の式（３）を実行する。なお、式（３）において、ヤコビ行列Ｊ_ｘ（θ）は、手先速度と関節速度とを関係づけるもの（写像）である。また、ベクトルｘは手先位置を示す。

筋シナジー算出部１０６は、数２の式（４）を実行する。

すなわち、筋骨格モデルに基づいた筋シナジー抽出法（宇野かんな、富永健太、奥貴紀、松居和寛、植村充典、平井宏明、宮崎文夫、``多様なヒト上肢運動間で共通する筋協調パターンの抽出” 第３１回日本ロボット学会学術講演会、2F2-03、首都大学東京、2013年9月4日）によると、図２（Ｂ）に示す３対６筋モデルにおいて、手先平衡点の極座標速度ｐ^・＝（ρ^・，φ^・）^Ｔは、筋拮抗比ｒ_ｉ(ｉ＝１，２，３）と筋拮抗和ｓ_ｉ(ｉ＝１，２，３）を用いて算出される。なお、ｕ_ρ（ｓ）,ｕ_φ（ｓ）は動径方向、偏角方向のシナジーベクトルであり、ｕ_ρ×φ（ｓ）はＮｕｌｌ方向のシナジーベクトルである。

平衡点算出部１０７は、数３の式（５）を実行する。

式（５）において、シナジーベクトルｕ_ρ（ｓ）、ｕ_φ（ｓ）、ｕ_ρ×φ（ｓ）と、筋拮抗比ｒの平均値ｒ⁻からの変化量ベクトルｒ−ｒ⁻との内積をシナジースコアと定義する。これらの関係から、極座標系での手先平衡点と筋拮抗比は、式（５）の関係を持つ。なお、Ｎｕｌｌ方向のシナジースコアに関しては，その大きさのみを評価する。図５に示す点ＥＰは、運動前後の静止時間において、ｕ_ρ（ｓ）・（ｒ−ｒ⁻）およびｕ_φ（ｓ）・（ｒ−ｒ⁻）がそれぞれ実測値の動径及び偏角と一致するようにレンジを調節し、ρ及びφを算出したのち、ｘ＝ρcosφ，ｙ＝ρsinφとして、ｘ−ｙ空間に変換したものである。

剛性楕円算出部１０８は、数４の式（６）で示される剛性楕円推定法１、あるいは式（８）で示される剛性楕円推定法２を実行する。

なお、式（６）の剛性楕円推定法１は、平衡点位置のＮｕｌｌ成分に関わるｗ_ρ×φ（ｓ, r）が関与していないが、式（８）の剛性楕円推定法２では、式（７）のように平衡点位置のＮｕｌｌ成分に関わるｗ_ρ×φ（ｓ, r）が関与している。図６に示す剛性楕円Ｄは、推定結果を動作点に対応付けて重畳表示したもので、λmaxは長径方向及び長さを示し、λminは短径方向及び長さを示している。径の長さは関節の動かし難さ（手先を各方向へ一定量だけ変位させるのに必要な外力の大きさ）を示している。

駆動ロボット５０は、例えば、機構体の所要部位を運動させるための複数のアクチュエータを有するもので、電気信号に応じて空圧を調整し、空気圧ゴムチューブを伸張、収縮可能にするMcKibben型人工筋肉等で構成されて、前記特徴量に応じた信号を受信し、動作を行う。アクチュエータはMcKibben型人工筋肉に限定されず、電気信号で、あるいは電気信号を変換した別の力によって機械的運動に変える電気機械変換素子、例えば電磁ソレノイド、圧電素子、モータ等であってもよい。制御対象が、人工筋アームを備えた筋骨格ロボット、例えば上肢筋骨格ロボットである場合、人体の場合と同様に目標手先運動を実現する人工筋活動コマンドを生成するメカニズムが必要となる。また、人体であるマスタ側からの出力信号に対応して、スレーブ機側であるロボットを制御する場合（いわゆるマスタースレーブシステム）には、ヒトとロボットとでは筋特性や構造が異なるため、ヒトの筋電位計測結果からリアルタイムで得られる特徴量をロボットの筋拮抗比や筋拮抗和に変換する。

図７は、運動解析処理の手順を示すフローチャートである。まず、筋電位測定部２０、運動測定部３０から筋電位の検出信号、運動位置の信号が取り込まれ（ステップＳ１１）、これらの信号から筋電位の測定処理、運動位置等の測定処理が実行される（ステップＳ１３）。次いで、測定された筋電位信号、運動位置等の信号から、前処理が実行される（ステップＳ１５）。すなわち、式（１）〜式（３）に従って、筋拮抗和、筋拮抗比の各値、及びヤコビ行列が算出される。

続いて、前準備で得られた各情報を用いて、３つの特徴量の算出が行われる（ステップＳ１７）。すなわち、筋シナジー、平衡点、及び剛性楕円の値が算出される。そして、これらの３つの特徴量は、表示処理部１５１に出力される（ステップＳ１９）。次いで、処理が終了か否かが判断され（ステップＳ２１）、終了であれば、本フローを抜け、一方、終了でなければ、ステップＳ１１に戻って同様の処理が繰り返される。

図８は、表示処理の手順を示すフローチャートである。まず、カメラ映像の受信、及び受信映像のモニタ１５２への表示が実行される（ステップＳ３１）。次いで、３つの特徴量の受信が行われると（ステップＳ３３）、筋シナジーの値が所定の図表に表現されてモニタ１５２の所定領域に表示される（ステップＳ３５）。また、平衡点及び剛性楕円の値が所定の図形（図１、図６参照）に変換されてモニタ１５２の所定領域に表示される（ステップＳ３７）。次いで、処理が終了か否かが判断され（ステップＳ３９）、終了であれば、本フローを抜け、一方、終了でなければ、ステップＳ３１に戻って同様の処理が繰り返される。なお、図８の表示処理は、信号処理部１０のコンピュータで処理する態様でもよい。

次に、実験例について説明する。

ヒトは、冗長な筋骨格系を適切に制御し、環境に柔軟に対応することで様々な運動を実現している。その制御メカニズムの解明は、これまで多くの研究者が取り組んできた課題である。特に、剛性の調節はヒトの運動戦略を理解する上で重要な要素の一つである。しかしながら、これらの制御則と剛性の調節メカニズムを統合する知見は未だ得られていない。この問題に対して、本発明者らはこれまで関節周りに拮抗的に配置される筋対（拮抗筋対）に着目し、拮抗筋対のＥＭＧの比を筋拮抗比、拮抗筋対のＥＭＧの和を筋拮抗和と定義した（特許文献２参照。ただし、本文献では筋拮抗和を筋活性度と呼称している。）。筋拮抗比が関節平衡点の調整に、また筋拮抗和が関節剛性の調整にそれぞれ寄与することから、筋拮抗比と筋拮抗和は、拮抗筋対の協調性・関節平衡点・関節剛性を関連付ける概念であるといえる。

そこで、以下では、筋シナジー仮説の観点から、運動中の手先剛性の調節メカニズム解明を試みる。実験は、リーチング運動中のＥＭＧと手先剛性を測定し、筋拮抗比と筋拮抗和の概念及び筋骨格モデルをもとに、筋シナジーを抽出した後、（１）抽出した筋シナジーを用いて推定した手先剛性と、従来技術（非特許文献１、２参照。）を用いて算出した手先剛性とを比較する。さらに、（２）抽出した筋シナジーを用いて推定した手先の平衡点軌道と、前記従来技術を用いて算出した手先の平衡点軌道との類似性を明らかにする。

＜実験例＞
（１）方法
（１．１）実験内容
被験者２名（Ａ，Ｂ）は、いずれも２０代男性、右利きである。被験者は右肩がマニピュレータ３の正面に来るように着席し、図２（Ａ）に示すように椅子の背もたれに両肩を固定された。被験筋は、水平面内の上肢運動で主要な役割を担う３対６筋（うち１対は二関節筋）であり、被験筋と座標系の定義を図２（Ｂ）に示す。

前腕と上腕は天井からつりさげられて重力の影響が除外された。また、プラスチックのギブスによって右腕の手首関節は固定され、レーザーポインタ４とマニピュレータ３（PHANTOM Premium 3.0/6DOF, SensAbleTechnologies 製，サンプリングレート１００[Hz]）先端が手先に固定された。作業台の上には、運動方向と距離を示す紙５を置き、レーザーポインタ４で被験者が手先位置をリアルタイムに視認できるようにした。光学トラッキングシステム（OptiTrack,NaturalPoint 社,サンプリングレート１００[Hz]）のためのマーカ８が左肩、右肩、右肘の関節中心と、右手先のレーザーポインタ直上に固定された。筋電位取得装置及び生体アンプ（WEB-5000, 日本光電社）とＡＤ変換器（PowerLab, AD Instruments Inc.) を用いて、ＥＭＧは１０００[Hz]で計測された。また、トラッキングシステムより得られる運動学データとマニピュレータ３が発生する力学データ、および筋電位取得装置から得られるＥＭＧデータは、動作点の平衡点や剛性の調整メカニズムの解明を試みるために、すべて同期して計測された。

水平面内のリーチングでは、（−０．２，０．４５）[m]をスタート地点とし、（０．２，０．４５）[m]をゴール地点とした。スタート地点とゴール地点とは被験者にとって自然な運動が可能な範囲に収めた。リーチング運動前は、スタート地点で１秒間静止し、その後１秒間でリーチング運動し、ゴール地点で１秒間静止した。これを１試行とした。“運動開始１[s]前”、“運動開始０．５[s]前”、“運動開始”、“運動終了”、および“運動終了後１[s]後”の５回のタイミングでビープ音を発し、運動の合図とした。また、試行と試行の間は、２±０．５[s]の間隔をランダムに設けた。

剛性計測のための摂動は、リーチング運動直前（最初のビープ音から０．９[s]後）、手先が行程の１／４，２／４，３／４を通過した瞬間、およびリーチング運動終了直前（最初のビープ音から１．９[s]後）の計５か所のうち１か所で与え、場所の選び方はランダムに設定した。摂動を与える時間は、短い時間（０．２[s]）とした。また、摂動はｘ軸方向から４５[deg]間隔の８方向に与え、どの方向に摂動を加えるかはランダムな順番に設定した。１セットの計測は４０試行（５箇所，８方向）行われ，これを８セット計測した。摂動の大きさは、被験者の体格にあわせて６[N]、または８[N]に設定した。また、摂動を与えない状態でのリーチング運動を４０試行行った。

（１．２）データ解析
（１．２．１）ＥＭＧデータ
ＥＭＧデータは１０〜４５０[Hz]のバンドパスフィルタ、整流化、および５[Hz]ローパスフィルタによる平滑化を施した後、最大随意筋収縮時（Maximum Voluntary Contraction, MVC）の値を用いて%ＭＶＣ化し、平均をとったものをｍ_ｉ(ｉ＝１，２，…，６）とした。ｍ_ｉを用いて定義した筋拮抗比と筋拮抗和およびその機能は、図３にまとめている。

（１．３）筋シナジーの抽出と平衡点軌道の推定
本発明者が提案している、筋骨格モデルに基づいた筋シナジー抽出法によると、図２（Ｂ）に示す３対６筋モデルにおいて、極座標系で表される手先平衡点の速度ベクトル
ｐ^・＝（ρ^・，φ^・）^Ｔは、筋拮抗比の速度ベクトルｒ^・＝（ｒ^・ _１，ｒ^・ _２，ｒ^・ _３）^Ｔと筋拮抗和ｓ_ｉ(ｉ＝１，２，３）とを用いて、式（９）で与えられる。なお、本実験では、簡単のため上腕と前腕のリンク長を同じとした。

であり、

である。ここで、ｃは筋の特性を表す定数、Ｌはリンク長、θ_２は肘関節角度である。これらの式は、手先平衡点の速度が、筋拮抗比速度ベクトルｒ^・＝（ｒ^・ _１，ｒ^・ _２，ｒ^・ _３）^Ｔの、（筋拮抗和の結合によって張られる）ベクトル空間ｑ₂×(ｑ₁+ｑ₂／2)への射影によって推定できることを示している。これをもとに、動径方向、偏角方向、動径と偏角に直交する方向の基底ベクトルを求めると、式（１１）に示すように、

となる
この基底ベクトルは筋拮抗比ベクトルの動径・偏角方向への分配を示すものであるから、ｕ_ρ（ｓ）、ｕ_φ（ｓ）をそれぞれ動径、偏角方向のシナジーベクトル、ｕ_ρ×φ（ｓ）をＮｕｌｌ方向のシナジーベクトルとした。これらのシナジーベクトルと、筋拮抗比ｒの平均値ｒ⁻からの変化量ベクトルｒ−ｒ⁻との内積をシナジースコアと定義した。これらの関係から、極座標系での手先平衡点と筋拮抗比は、式（１２）で示すような関係を持つことがわかった。

ここでは、運動前後の静止時間においてｕ_ρ・（ｒ−ｒ⁻）およびｕ_φ・（ｒ−ｒ⁻）がそれぞれ実測値の動径および偏角と一致するようにレンジを調節し、ρおよびφを算出したのち、
ｘ_ｅｑ ^ｅｍｇ＝（ρcosφ，ρsinφ）としてｘ−ｙ空間に変換した。

なお、Ｎｕｌｌ方向のシナジースコアに関しては、その大きさのみを評価した。

（１．３．１）ＥＭＧデータを用いた手先剛性の推定
筋骨格モデルをもとに、筋拮抗和ｓ_i（ｉ＝１，２，３）を用いて関節剛性行列Ｋ_θ ^{ｅｍｇ＿ｓ}（ｓ）を表すと、式（１３）で示すような、

となることが過去の研究（前出の論文）に示されている。ただし、ｋ_ｓは無次元数の筋拮抗和を関節剛性に変換するための定数である。ここで、特に動的な状況下では、Ｎｕｌｌ方向のシナジーベクトルが剛性の大きさと向きを調整しているという仮説を検証するため、関節剛性行列Ｋ_θ ^{ｅｍｇ＿ｓ}（ｓ）の各成分に、Ｎｕｌｌ方向のシナジースコアｗ_ρ×φ（ｓ，r）＝ｕ_ρ×φ（ｓ）・（ｒ−ｒ⁻）を足し合わせたものを、式（１４）で示すように、ＥＭＧデータより推定した関節剛性行列Ｋ_θ ^{ｅｍｇ＿ｓ＋ｎ}（ｓ，r）とした。なお、ｋ_ｎの値は、実験的に２００[Ｎｍ／ｒａｄ]とした。

関節空間の剛性行列Ｋ_θ ^{ｅｍｇ＿ｓ＋ｎ}（ｓ，r）を、式（１５）によりタスク空間の、
剛性Ｋ_ｘ ^{ｅｍｇ＿ｓ＋ｎ}（θ，ｓ，r）に変換したものを、ＥＭＧデータから推定した手先剛性行列とする。

なおＪ_ｘ（θ）は手先速度と関節速度を関係づけるヤコビ行列である。

（１．３．２）摂動法による手先剛性と平衡点軌道の推定
手先の運動は、リーチング運動（被験者が本来行おうとした運動）と、摂動による変位量の２つが混在した状態で計測された。リーチング運動からマニピュレータ３による摂動の方向、摂動を与えるタイミング、および摂動の有無の影響を取り除くため、全データで平均を取った関節角度ベクトルθ（ｔ）＝（θ_１（ｔ），θ_２（ｔ））^Ｔをリーチング運動とした。

一方、摂動による変位量、δθ（ｔ）＝（δθ_１（ｔ），δθ_２（ｔ））^Ｔは、式（１６）に示すように、摂動開始時からの変位量ベクトルθ_ｓ，ｉ（ｔ）から、摂動開始時からの変位量の摂動方向平均ベクトルθ_ｓａｖ（ｔ）を減じたものとした。

ただし，ｉ＝１，２，…，８は摂動の方向を表し、ｔ_ｓ≦ｔ≦ｔ_ｓ＋Ｔとした。ここに、ｔ_ｓは摂動開始の時刻、Ｔは剛性推定計算に利用する区間の移動に要する時間で０．４[s]とした。摂動を与える際、マニピュレータ３には一定の手先力を発生するように指令したため、前記区間の移動時間中は一定の外力が発生したとみなした。この関係を式（１７）で示した。

ｆの大きさは被験者Ａで８[N]、被験者Ｂで６[N]であった。この外力を、関節空間に変換したものを摂動トルクδτ（ｔ）＝Ｊ^Ｔ δｆ（ｔ）とした。

θ（ｔ）、δθ（ｔ）、およびδτ（ｔ）から、慣性行列Ｉ（θ）、
関節粘性行列Ｄ＝[[Ｄ_ｓｓ，Ｄ_ｓｅ]^Ｔ，[Ｄ_ｅｓ，Ｄ_ｅｅ]^Ｔ]、
および関節剛性行列Ｋ_θ ^ｐｔｂ＝[[Ｋ_ｓｓ，Ｋ_ｓｅ]^Ｔ，[Ｋ_ｅｓ，Ｋ_ｅｅ]^Ｔ]を動力学方程式により算出した（Hiroaki Gomi, Mitsuo Kawato,“Human arm stiffness and equilibrium-point trajectory during multijointmovement,” Biological Cybernetics, vol.76,no.3, pp.163-171, 1997.参照）。なお、手先剛性行列、
Ｋ_ｘ ^ｐｔｂ（θ）は、関節剛性行列Ｋ_θ ^ｐｔｂとヤコビ行列Ｊ_ｘ（θ）および前記動力学方程式より算出される物理パラメータから推定した腕が手先に発生する力Ｆ_ｉｎとを用いて、式（１８）により求めた。

算出した関節剛性行列Ｋ_θ ^ｐｔｂおよび関節粘性行列Ｄを用いて、手先の平衡点軌道
ｘ_ｅｑ ^ｐｔｂを、式（１９）に従って算出した。

ここで，Φは関節角座標を手先のｘ−ｙ座標に変換する順運動学関数であり、Ｉ（θ），Ｈ（θ^・, θ）は先の方法で求めた物理パラメータと、θ＝（θ_１，θ_２）^Ｔより算出した。また，ＤとＫ_θ ^ｐｔｂは、手先位置５か所ごとに求めた値を線形補間し、それぞれ２０[Hz]のデータに直して計算した。

（２）結果
（２．１）筋拮抗比と筋拮抗和、および筋シナジー抽出
被験者Ａの実験から得られた筋拮抗和ｓ_ｉと筋拮抗比ｒ_ｉの推移を、図９に示す。図９（Ａ）は筋拮抗和ｓ_ｉを示し、図９（Ｂ）は筋拮抗比ｒ_ｉを示す。これらの筋拮抗和ｓ_ｉと筋拮抗比ｒ_ｉをもとに得られた動径、偏角、Ｎｕｌｌ方向のシナジーベクトルｕ_ρ（ｓ），ｕ_φ（ｓ），ｕ_ρ×φ（ｓ）と、それぞれのシナジースコアｕ_ρ（ｓ）・（ｒ−ｒ⁻），ｕ_φ（ｓ）・（ｒ−ｒ⁻），｜ｕ_ρ×φ（ｓ）・（ｒ−ｒ⁻）｜を、図１０に示す。図１０（Ａ）は被験者Ａのデータであり、図１０（Ｂ）は被験者Ｂのデータである。なお、図中、左側の各シナジーベクトルを示す棒グラフは、左側から、
ｒ_１、ｒ_２、ｒ_３に対応する。被験者Ａ，Ｂ共に類似した結果が見られる。

（２．２）関節剛性と手先剛性
摂動法により算出された関節剛性を、図１１に示す。図１１（Ａ）は被験者Ａのデータであり、図１１（Ｂ）は被験者Ｂのデータである。次に、摂動法より求めた関節剛性を手先剛性に変換したものと、筋拮抗和およびＮｕｌｌ方向のシナジースコアから推定した手先剛性を、図１２に示す。図１２（Ａ）は摂動法より求めた特性データであり、図１２（Ｂ）は今回の推定法から求めた特性データである。また、比較のため、Ｎｕｌｌ方向のシナジースコアを足し合わせずに求めた手先剛性（数１６の式（２０）に示す）

を、図１２（Ｃ）の特性図に示す。

（２．３）手先の平衡点軌道と速度分布
関節剛性行列Ｋ_θ ^ｐｔｂ等を用いて算出された手先の平衡点軌道ｘ_ｅｑ ^ｐｔｂと、筋シナジースコアから推定した手先の平衡点軌道ｘ_ｅｑ ^ｅｍｇ、およびそれらの接線速度の大きさ｜ｘ^・ _ｅｑ ^ｐｔｂ｜、｜ｘ^・ _ｅｑ ^ｅｍｇ｜を、図１３に示す。図１３（Ａ）は摂動法より求めた特性データであり、図１３（Ｂ）は今回の推定法から求めた特性データである。

（３）考察
（３．１）抽出された筋シナジーについて
本解析法で抽出された筋シナジーベクトルは、図１０の左側の図に示す通り、２名の被験者間でそれぞれ高い一致を示した（内積０．９４±０．０２）。これにより、この筋シナジー抽出法には被験者間の個人差がほとんど反映されないと考えられ、換言すれば一般性が高いことがわかった。

また、図１０の左側の図より、動径方向のシナジーベクトルｕ_ρ（ｓ）の働きに着目すると、ヒトが動径の正方向に手先を運動させるとき、二関節筋対は肘の伸展をサポートする役割を担っていることがわかった。一方、偏角方向のシナジーベクトルｕ_φ（ｓ）の働きに着目すると、ヒトが偏角の正方向の手先運動を考えるとき、肩の単関節筋対と肩肘の二関節筋対が同調して肩関節の運動に作用しているのに対して、肘の単関節筋対は作用していないことがわかった。

それに対し，Ｎｕｌｌ方向のシナジーベクトルｕ_ρ×φ（ｓ）は、肩と肘の単関節筋対と、肩肘の二関節筋対が拮抗している状態を表している。筋拮抗比ベクトルｒがＮｕｌｌ方向に変化しても、手先平衡点は変化しないと考えられる。しかし、図１０の右側の図に示すように、Ｎｕｌｌ方向のシナジースコアはリーチング運動の開始直後と終了の直前に大きく変動し、またこの傾向は被験者間で共通して観察された。したがって、Ｎｕｌｌ方向のシナジーベクトルは、手先平衡点の調節以外の役割を担った筋シナジーであることがわかった。

（３．２）手先剛性と筋拮抗和およびＮｕｌｌ方向のシナジースコア
摂動法により求めた手先剛性は、図１２（Ａ）に示すように、運動の開始直後と終了直前に大きさが増加し、運動の中盤では大きさが減少した。

また、運動中盤の剛性楕円は、肩と手先を結んだ方向からやや傾いており、手先剛性の方向が運動中に変化していることが分かった。また、運動中における手先剛性の方向の変化特性は、いずれの図においても共通している。また、運動中における、手先剛性を示す楕円形状の大小変化特性も、図１２（Ａ）（Ｂ）間で共通している。

また、筋拮抗和とＮｕｌｌ方向のシナジースコアに着目すると、運動開始直後と運動終了直前の２か所でピークを持つ点で、手先剛性の大きさの推移と特徴が類似している。図１２（Ｂ）に示すように、筋拮抗和とＮｕｌｌ方向のシナジースコアとを用いて推定した手先剛性は、図１２（Ｃ）の筋拮抗和のみを用いて推定した手先剛性楕円に対して、より実測値に近い形状の推移をたどっている。したがって、手先剛性の大きさと方向の調節に、筋拮抗和と筋拮抗比のＮｕｌｌ方向のシナジースコアが共に寄与していると考えられる。

（３．３）手先の平衡点軌道および速度分布と動径・偏角方向のシナジースコア
実際の速度分布がベル型であるのに対し、関節剛性を用いて求めた手先の平衡点速度は、図１３に示すように多峰性を示した。また、この多峰性は、動径方向と偏角方向のシナジースコアから推定した手先平衡速度にも表れた。したがって、本手法で抽出した筋シナジーから求めた平衡点軌道は、手先剛性から求めた平衡点軌道と類似した特徴を持つことがわかった。これにより、動径方向のシナジーベクトルｕ_ρ（ｓ）および偏角方向のシナジーベクトルｕ_φ（ｓ）は、手先平衡点の調節に寄与していると考えることができる。

以上説明したように、本発明に係る解析方法及び装置は、上肢を机上で運動させる態様で説明したが、運動部位は、上肢の他、下肢でもよく、また上肢と下肢との双方を連動させた全身運動を含む広範な運動にも適用可能である。すなわち、本解析方法及び装置は、広範な運動に対し中枢神経系が選択する運動戦略（運動学習度や個人差に依存）を明確な基準の下、客観的かつ定量的に使用者へ提示し、運動学習支援のための具体的な手段を提供し得る。動作例としては、用途に応じて種々想定され得るが、各種スポーツの動きや技の学習、また、スポーツ以外の動きや立ち振る舞いの習得等にも適用可能である。すなわち、目標課題に対する効果的な運動訓練法や運動介入法の確立へ向けた運動学習支援方法、装置への適用では、例えば、予め模範者の運動時に特徴量を取得しておき、これを参考にして練習者は自己の特徴量を模範者に近づけるよう練習する練習支援である。あるいは、自己の特徴量の練習毎の変化を観察しながら練習する練習支援である。

また、前述の解析方法及び装置は、健常者や例えば脳疾患患者の運動生成機序や運動学習機序の解明へ向けた運動制御の評価方法、装置にも適用可能である。中枢神経系の運動指令に基づき推定される上記３つの特徴量から運動技能を診断評価し、附加的フィードバックを介して運動学習の促進を図る態様は、スポーツ分野及びリハビリテーション分野において有効的である。このように、本解析方法により得られた３つの特徴量は、運動制御評価、運動学習支援に効果的に利用され得る。

また、ヒトと物理的な相互作用を伴う機械、道具及び環境の設計と評価を目的としたインタフェース評価装置にも適用可能である。

１０信号処理部（演算手段）
１０３筋拮抗和算出部（第１の演算手段）
１０４筋拮抗比算出部（第１の演算手段）
１０５ヤコビ行列算出部（第２の演算手段）
１０６筋シナジー算出部（演算手段）
１０７平衡点算出部（演算手段）
１０８剛性楕円算出部（演算手段）
２０筋電位測定部（測定部）
３０運動測定部（測定部）
４０カメラ（撮像部）
１５表示部（出力部）
５０駆動装置（出力部）

Claims

運動者の筋活動及び身体運動を周期的に測定する測定部と、
前記測定部で得られた測定情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を筋骨格モデルに基づき算出する演算手段とを備えた運動解析装置。
前記測定部は、前記筋活動として筋電位を測定する筋電位検出部と、前記身体運動として関節の動きを測定する運動検出部とを備え、
前記演算手段は、筋拮抗比及び筋拮抗和を算出する第１の算出手段と、前記筋骨格モデルに基づいて身体運動の測定部位を前記動作点に写像するヤコビ行列を算出する第２の算出手段とを備え、算出された筋拮抗比、筋拮抗和及びヤコビ行列を用いて前記特徴量を算出することを特徴とする請求項１に記載の運動解析装置。
算出結果を出力する出力部を備えた請求項１又は２に記載の運動解析装置。
前記身体運動を観察する撮像部を備え、
前記出力部は、画像を表示する表示部を備え、前記表示部に、前記特徴量のうち、少なくとも前記平衡点及び前記剛性楕円の表示と、前記撮像部で撮像された身体運動の画像とを表示することを特徴とする請求項３に記載の運動解析装置。
前記出力部は、前記平衡点を前記身体運動の画像に対応させて表示することを特徴とする請求項４に記載の運動解析装置。
前記出力部は、前記剛性楕円を前記身体運動の画像の動作点に対応させて表示することを特徴とする請求項４又は５に記載の運動解析装置。
前記出力部は、前記特徴量に基づき支援対象に物理的支援を施す駆動部であることを特徴とする請求項３〜６のいずれかに記載の運動解析装置。
測定部によって運動者の筋活動及び身体運動を周期的に測定し、
前記測定部で得られた測定情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を筋骨格モデルに基づき算出する運動解析方法。
前記測定部で周期的に測定された運動者の筋活動及び身体運動の情報を取り込む測定処理手段、前記測定処理手段で取得された情報から、動作点における剛性楕円、平衡点、及び平衡点を記述する基底ベクトルである筋シナジーの各特徴量を筋骨格モデルに基づき算出する演算手段、としてコンピュータを機能させる運動解析プログラム。