計測自動制御学会東北支部 第 244 回研究集会 (2008.7.2)
資料番号 244-3

筋電データを考慮した下肢筋張力推定法
Muscle tension estimation of the lower extremities
using EMG data
○塙秀祐 ∗ ，山藤和也 ∗ ，平元和彦 ∗∗ ，巖見武裕 ∗∗ ，島田洋一 †○ Syusuke Hanawa∗ , Kazuya Santo∗ , Kazuhiko Hiramoto∗∗ ,Takehiro Iwami∗∗ , Yoichi Shimada†*秋田大学大学院工学資源学研究科, **秋田大学工学資源学部，† 秋田大学医学部* Graduate School of Engineering and Resource Science, Akita University,** Faculty of Engineering and Resource Science, Akita University,† Faculty of Medicine, Akita University
キーワード : 筋張力推定 (Muscle tension estimation)，筋電データ (EMG data)，線形行列不等式 (LMIs: Linear
matrix inequalities)，リハビリテーション (Rehabilitation)，ローイング運動 (Rowing)

連絡先 : 〒 010-8502 秋田市手形学園町 1-1 秋田大学工学資源学部機械工学科
平元和彦，Tel.: (018)889-2348，Fax.: (018)837-0405，E-mail: hira@ipc.akita-u.ac.jp

1.

はじめにし，それらのリンクを回転関節で接続したモデルが多く用いられている．この剛体リンクモデルを医療における効果的かつ安全なリハビリテーショ用いて，筋張力は以下の手順で推定されてきた．

ン計画や，スポーツにおける効率的なトレーニン
グ計画の策定のため，運動時の身体各部にかかる

Step 1: 身体各部の運動や力を計測し，解剖学的なデータから得られる身体部位パラメータ値負担を解析することは非常に重要である．

を用いて各関節の関節モーメントを求める．

身体各部への負担を解析するためのバイオメカニクスからのアプローチとして，身体筋骨格モデ

Step 2: 張力を推定する筋を，剛体リンクモデル

ルに基づく関節モーメントの導出や，筋張力の推の関節やリンクに付着させ，step 1 で求めた定を挙げることができる．特に，筋張力は，関節関節モーメントを発生するような筋張力をへの負担の定量的指標である関節間力を求めたり，推定する．

運動時に筋の消費するエネルギーを計算したりするために必須な量であり，その正確な推定が求め

Step 2 において，関節まわりのモーメントの釣合い式の数に比べて，推定する筋張力の数は一般にられている．
筋張力を推定するための身体筋骨格モデルとして，運動に関連する身体各部を剛体リンクと仮定多くなるため，推定のために（効率的な筋活動を表現する何らかの）評価関数を設け，それを最適化することによって筋張力を推定してきた 1) ．評

–1–



価関数としては，筋活性度や筋応力の二乗和または三乗和が一般的で，従来主に傾斜法 たとえば 2) をA = diag(a1 , . . . , an ): 対角行列A = 

a1..


．
.
an

用いて張力推定が行われていたが，得られる結果には，以下のような問題点があった．

1) 前述の最適化手法では，筋張力は正の値をとるという制約しか考えられていなかった．

2.

ローイング運動のモデル化

2.1

実験環境

その結果，考慮している運動において，筋電本報告で取り扱うローイング運動によるリハビ

データから拮抗筋作用が認められたとしてリテーション機器の構成を図1に示す．装置の全長，も，その結果を最適化の際考慮できなかった．

全幅はそれぞれ 2.492[m] および 0.600[m] であり，ト2) 発生可能な筋張力は，筋の長さや伸縮速度にレーニング時の座面の高さは 0.303[m] である．

依存していることが知られていたが，これらマシン本体に，被験者の手と足の力を測定するを最適化計算の際考慮した例は少なかった．

為に 6 軸力覚センサを取り付ける．6 軸力覚センサ

本報告では，上記の問題点を解決する新しい筋張力推定法を提案する．はじめに，筋張力を推定するための最適化問題が，筋張力に関するLMI（LinearMatrix Inequality: 線形行列不等式）問題に帰着され，大域的な最適解を効率的に求めることが可能であることを示す．この結果を用いることによって，筋張力を求める際，計測された筋電データから従来の推定手法では張力の発生しない筋に（筋は，互いに直交する 3 つの軸方向の力と各軸回りのトルクを計測することができる．データはセンサ部から直接 DSP ボードへ送られ，データ処理用PC によって記録される．
3 次元運動解析システムによって，被験者の運動を記録する．被験者の身体表面には標点として反射マーカを取り付け，6 台のカメラで運動を記録し，そのデータから各標点の 3 次元座標データを計算できる．

電データを反映した）張力を発生させながら筋張力推定を行うことができる．さらに，筋の長さの一部の筋に関しては，無線式生体計測装置を使用して筋電位を計測し，そのデータを記録する．

変化やその伸縮速度を考慮した張力推定も可能である．本報告では，安全かつ負担の小さいリハビリテーションメニューとして最近注目されている

2.2

剛体リンクモデル

運動の解析に用いる剛体リンクモデルを図 2 にローイング運動3) を取り上げ，提案する手法によっ示す．足関節，膝関節および股関節モーメントをて筋張力推定を行い，結果を考察する．
本報告の構成を以下に示す．2 節では，本報告で取り扱うローイング運動によるリハビリテーション機器の実験環境を紹介し，解析に用いる剛体リンクモデルを示す．3 節で，本報告で扱う筋張力推定のための筋骨格モデルを示し，問題設定と推定手法を述べる．4 節で，実験データに基づく筋張力推定の結果を示し，5 節でまとめを行う．
本報告で用いる記号は，以下の通りである．g: 重力加速度．AT : 行列Aの転置．In : n次元単位行列．

–2–Fig. 1

ローイングマシン

Fig. 2

剛体リンクモデル

Fig. 3筋骨格モデルそれぞれ M a , M k , M h とすると，各関節モーメン

9 個の筋群（1. 大腿直筋，2. 大殿筋，3. ハムス

トは次式で求められる．

トリングス，4. 広筋，5. 大腿二頭筋長頭，6. 腓腹筋，7. ヒラメ筋，8. 前脛骨筋，9. 腸腰筋）をM a = I1 θ̈1 + fp (yp − y1 ) + fa (y1 − ya )− np (xp − x1 ) − na (x1 − xa )考慮する．それぞれの筋に発生する筋張力，その(1)M k = I2 θ̈2 − fa (y2 − ya ) − fk (yk − y2 )− na (xa − x2 ) − nk (x2 − xk ) + Ma

(i = 1, . . . , 9) とおく．各々の筋が付着する関節の(2)レバーアーム長を Lji (i = 1, . . . , 9, j = a, k, h) とする．なお，筋が付着していない関節に関するレM h = I3 θ̈3 + fk (yk − y3 ) + fh (y3 − yh )− nk (xk − x3 ) − nh (y3 − yh ) + Mk

最大および最小値をそれぞれ Fi , Fimax および Fimin(3)

バーアーム長は 0である．このとき，筋張力と関節モーメントの釣合い関係式は，次式で与えられる．

fa = fp − m1 ẍ1 , na = np − m1 (ÿ1 + g),

Mj =fk = fa − m2 ẍ2 , nk = na − m2 (ÿ2 + g),

9∑

Fi sji Lji , j = a, k, h(4)i=1

fh = fk − m3 ẍ3 , nh = nk − m3 (ÿ3 + g),

ここで，sji (i = 1, . . . , 9, j = a, k, h) は，関節に作ここで，mj , Ij , xj , yj (j = 1, 2, 3) は，それぞれ足用するモーメントの符号によって 1 または −1 の値

部，下腿部，大腿部の質量，慣性モーメント，重心を取る定数である．本報告で筋張力を推定する筋位置の x および y 座標であり，文献 4) のデータを用の数は 9 個であるが，利用可能な条件式は式 (4) でいて体重から推定されたものを用いている．xi , yi表される 3 つしかない．このため，このままでは，fi , ni (i = p, a, k, h) は，それぞれ計測された COP

筋張力推定問題は数学的に不定問題となり，式 (4)（圧力中心），足関節，膝関節および股関節の x, yを満足する筋張力 Fi (i = 1, . . . , 9) の組み合わせを座標および各部に作用する並進力の x, y 方向成分唯一に求めることは不可能である．従来研究では，である．

唯一の筋張力の組み合わせを求めるために，運動に対する効率的な筋活動を表すなんらの評価関数

3.

筋張力推定法3.1

筋張力推定問題を設定し，その評価関数を最適化して唯一の筋張力の組み合わせを求めたものが多く，本報告でもそのアプローチを採用する．本報告では，評価関前節で求められる関節モーメントを用いて，本数を次式のようなある時刻における筋活性度の二報告では，図 3 の筋骨格モデルに示されるような–3–乗和 J と定義する．

J :=さらに，等式拘束条件である関節モーメントに関)29 (∑FiFimax
i=1

する釣合い条件（式 (4)）は，微小正数 0 < ϵ ≪ 1(5)を用いて，G に関する以下の線形不等式拘束条件に近似できる．

以上の準備の下で，本報告では，下肢筋張力推定問題を以下のように定式化する．
¶下肢筋張力推定問題−ϵ ≤ M j − GT HN j ≤ ϵ, j = a, k, h³(9)ここで，H := diag(F1max , . . . , F9max )であり，Nj :=[]T(j = a, k, h) である．式 (7),sj1 Lj1 . . . sj9 Lj9ローイング運動開始時刻から終了時刻までの区間でサンプルされた時刻それぞれに対して，式(4)で与えられるモーメントのつりあいに関(8), (9) は，未知筋張力 Fi (i = 1, . . . , 9) を含む未知する等式拘束条件および筋張力に関する不等式拘束条件 Fimin ≤ Fi ≤ Fimax (i = 1, . . . , 9)

変数ベクトル G に関する連立 LMI である．よって，の下で，式 (5)J を最小化する筋張力 Fi (i =拘束条件式 (8) および (9) を満たしながら，式 (7)µ1, . . . , 9) を求めよ．
µ´を大域的に最適化するような G を，効率的解法の実装されたソフトウエアパッケージ 6) を用いて求めることが可能である．言い換えると，前節で設

3.2

筋張力推定法

3.2.1

LMI 問題への帰着定された下肢筋張力推定問題に対して，式 (5) を大域的に最適化する筋張力 Fi (i = 1, . . . , 9) の唯一の式(5)は，以下のように書き換えることができる．

J = GT G(6)組み合わせを，必ず求めることが可能である．

3.2.2

ここで，G は，次式で定義される 9 次元の列ベクトルである．

[G :=F1F1max

F2F2max...

F8F8max

F9F9max

筋電データの考慮

従来の筋張力推定は，式(5)のような評価指標を，傾斜法などによって最適化する際，筋 i に発生する]T筋張力の最小値 Fimin を 0 としていた．本報告では，筋電データが取得されている筋 i に対して，以下がベクトル G は，筋張力 Fi (i = 1, . . . , 9) の線形関数成立していると仮定する 7) ．

であることに注意する．式(6)およびSchur comple-ment の補題 5) より，J がある正数 µ > 0 よりも小さくなる条件 µ − GT G > 0 は，以下のような G（Fi ,i = 1, . . . , 9）および µ > 0 に関する LMI（線形行列不等式）に帰着される．
][µ GTGI9qi =FiEMGi=Fimax
EMGmaxi(10)ここで，qi , EMGi および EMGmax
は，それぞれ筋ii の活性度（0 ≤ qi ≤ 1），筋電データおよび最大筋電を表す．式 (10) は，筋活性度を筋張力と筋電>0(7)さらに，Fi に関する不等式拘束条件は，G に関する以下の線形不等式拘束条件に書き換えられる．

Gmin ≤ G ≤ Gmax ,[ min]Tmin
Gmin := FF1max . . . FF9max,
19[]TGmax := 1 . . . 1(8)で評価した場合，両者が等しくなっているという意味がある．式 (10) を変形すると次式となる．
)(EMGiFimax(11)
Fi =EMGmaxi
本報告では，筋電データ及び最大筋張力から得られる式 (11)のFi を，最適化時に筋 iの最小筋張力 Fimin
とおく．これにより，筋 i の推定される筋張力は，()EMGi
少なくとも EMG
Fimax 以上となり，（式 (10) のmaxi–4–意味で）筋電データを考慮した筋張力推定が可能

k(ξi)=0.32+0.71exp{−1.112(ξi−1)} sin{3.722(ξi−0.656)}
1.2となる．
10.8

筋長・筋伸縮速度の考慮
k(ξi)

3.2.3

筋張力は，その自然長からの変動や，伸縮速度0.4によって変化することが知られている．長谷 8) お0.2

よびその引用文献 9) では，筋 i の筋張力 Fi は，以下00.6のように与えられている．

Fi = k(ξi )h(ηi , ξi )Fimax qi

0.6

0.8(12)

11.2
ξi

1.4

1.61.8正規化筋長 ξi と筋力の関係 k(ξi )Fig. 4

ここで，Fimax は，解剖学的見地から決定される筋h(ηi,ξi)={1+tanh(a1ηi)}/a2−a3exp{−2.6(ξi−1)}iの長さや速度に依存しない定数の最大筋張力である．ξi および ηi は，それぞれ次式で定義される筋 i

2.5の正規化された筋長および筋伸縮速度である．

1.5ξ˙iliξi = n , ηi = max , vimax = 3.0livih(η,ξ)

21(13)

0.5

式 (13) において，li および lin は，それぞれ筋 i の長−0.52011.5さおよび自然長であり，vimax は最大筋伸縮速度で0.501ξある．関数k(ξi )およびh(ηi , ξi )は，次式で定義され

i0.5−0.5−1ηiている 9) （図 4, 5）．

Fig. 5k(ξi ) = 0.32 + 0.71e{−1.112(ξi −1)}
× sin {3.722(ξi − 0.656)}h(ηi , ξi ) =正規化筋長 ξi , 正規化筋伸縮速度 ηi と筋力の関係 h(ηi , ξi )(14)

1 + tanh(a1 ηi )− a3 e{−2.6(ξi −1)} (15)a2a1 , a2 , a3 : 筋繊維の種類によって定まる定数

4.

結果

4.1

運動データと関節モーメント

各筋の ξi は，直立状態で筋は自然長になり，筋は被験者は健常男性1名であり，3秒間でローイング

図 3 のように滑車に巻きついて自由に伸縮可能での 1 ストロークを完了するように指示されている．

あると仮定することにより，運動データより得ら運動を計測するための標点となる反射マーカは，れる各関節角度から求めることができる．さらに，頭頂点，第 7 頚椎棘突起（C7），胸骨剣状突起と右ηi は ξi を差分近似して求めることができる．本報

半身の鎖骨近位端，肩峰点，上腕骨外側上顆，橈

告では，式 (12) を考慮して，ある姿勢および運動骨茎状突起，尺骨茎状突起，第 5 中手骨頭，大転子速度（ξi , ηi ）での最大筋張力 Fimax (ξi , ηi ) を次式と点，膝外側点，足外踝点，弟 5 中足骨頭点の計 13 箇する．

所に貼付した．
計測データより，足関節，膝関節および股関節

Fimax (ξi , ηi ) := k(ξi )h(ηi , ξi )Fimax(16)角度（θa , θk , θh ）を計算したものを図6に示す．2.2

最大筋張力を式 (16) のようにおくことで，筋長や節の剛体リンクモデルから，各関節モーメントを筋伸縮速度を考慮した筋張力推定が可能になる．

式 (1)-(3) から計算したものを図 7 に示す．

–5–q1

θ [deg]150a100200150100500212334time [s]4time [s]5678q515678q6kθ [deg]

50
00.1
0.05
001234time [s]56780.1
0.05
001234time [s]56780.1
0.05
001234time [s]56780.1
0.05
001234time [s]5678θ [deg]150hq810050012Ma [Nm]Fig. 620100−100

34time [s]5678Fig. 8 筋 電 デ ー タ か ら 計 算 さ れ た 筋 活 性 度
qi = EMGi /EMGmax(i = 1, 5, 6, 8)i

関節角度 θi (i = a, h, k)

4.3
1234time [s]567筋長，筋伸縮速度と最大筋張力

8

図 3 の筋骨格モデルを仮定し，ローイングを行

M [Nm]

100

わせた場合の各筋の筋長変化の様子を図 9 に示す．

k503.2.3節でも述べたように，各筋の自然長は直立状態

001234time [s]5678における筋長とし，三次元画像解析ソフトMimicsr

Mh [Nm]50を用いて，MRI によって断層撮影された画像から，0下肢の筋の領域を抽出・積層して作成した筋モデ−5001234time [s]5678ルのデータを参考にして決定された．筋長は，図 3のように各筋は滑車に巻きついて自由に伸縮可能

Fig. 74.2

関節モーメント M i (i = a, h, k)であると仮定することにより求めた．筋伸縮速度は，得られる筋長をデータ取得サンプル周期 T =筋電データと最小筋張力1/60[s] で差分を取って求めた．式 (16) で計算されローイング運動時の 1. 大腿直筋，5. 大腿二頭筋長頭，6. 腓腹筋および 8. 前脛骨筋で取得された表面筋電データと，別途行われた最大力試験で得られたそれらの筋の最大筋電の値を用い，式 (10)る最大筋張力 Fimax (ξi , ηi ), i = 1, . . . , 9 を図 10 に示す．解剖学的見地から決定される一定の最大筋張力値 Fimax , i = 1, . . . , 9 も図 10 に合わせて示す．それらの値は，それぞれ以下のようである．

を用いて求めた筋活性度 qi , i = 1, 5, 6, 8 を図 8 に示す．この場合の筋活性度は，概ね 10%未満であり，F1max = 975, F2max = 4500, F3max = 3000,比較的低い．この値を用いて，最適化の際に用いF4max = 2250, F5max = 750, F6max = 1875,られる最小筋張力 Fimin , i = 1, 5, 6, 8 が式 (11) からF7max = 2100, F8max = 1125, F9max = 2025 [N]計算される．

–6–0.4
12345678l2 [m]0.2Fmax[N]
200.4

10000123456780.420.412345678

l4 [m]

0.3

Fmax
[N]
400.35123456780.3
123456780.3
0.2012345678

l7 [m]

0.35
0.3
0.25

12345678

l8 [m]

0.35
0.3
0.25

1234

56

78

l9 [m]

0.4
0.3
0.2

1234
time [s]5677801234

56

7

8

01234

56

7

801234567801234

567

801234567801234567801234567801234
time [s]56783000
2000

3000
2000

1500
1000

3000
2000
1000

06800

500

0510001000

042000

1000
Fmax
[N]
7l6 [m]

0.4

Fmax
[N]
6032500

600

Fmax
[N]
80.2

Fmax
[N]
9l5 [m]

0.4

Fmax
[N]
502

3000

1500

0.25

13500

2500

0.38

0

5000

00l3 [m]

1200

800

0.3

Fmax
[N]
3l1 [m]Fmax[N]
1Muscle length
Narural length

0.58Fig. 10
Fig. 94.4

筋長 li , i = 1, . . . , 9筋長，筋伸縮速度を考慮した最大筋張力

Fimax , i = 1, . . . , 9（式 (16)）

筋張力推定ていることが読み取れる．この結果は，筋電計測前節までの結果を用いて，筋張力を推定した．の結果を筋張力推定に取り込むことが可能となっ結果を図 11 に示す．比較のために，従来手法と同ていることを示しており，提案する手法による筋様の設定である，Fimin = 0, ∀ i = 1, . . . , 9（筋電張力推定が，従来推定法より実験で得られる筋電データを考慮しない）かつ Fimax = Fimax （筋長やデータを反映しながら張力推定を可能にしている

筋伸縮速度によって最大筋張力を変化させない）ことを示している．なお，他に表面筋電を計測しの下で式(5)を最適化して推定された筋張力を図12た筋（1. 大腿直筋および 8. 前脛骨筋）でも，筋電に示す．

データに基づく不等式拘束が満足された筋張力計算が行われている．

4.5

さらに，最大筋張力を筋長や筋伸縮速度を考慮考察して変化させた提案手法は，1. 大腿直筋で発生
図 11 および 12 より，筋電データを考慮する事によって，従来の方法では張力の発生していない 5.
大腿二頭筋長頭および 6. 腓腹筋に張力が発生しする筋張力がやや減少する傾向が見られ，とくに，ストロークとストロークの間の股関節と膝関節が屈曲から伸展に切り替わる区間（2.5-3.0, 5.5-6.0[s]

–7–

300
F [N]

200
100001234567F [N]0123456708F [N]012345674F4 [N]F [N]0123456708F [N]505F5 [N]0123456701234567801234567801234567801234567801234567801234567801234567801234time [s]567850
08150
F [N]

100

6F6 [N]

8100150
5001234567100
50
08150
F [N]

150
100

7F7 [N]

71000

100

50
01234567100
50
08300
F [N]

300
8200
10001234567200
100
08500

500

9F [N]

F8 [N]

620001000

F9 [N]

55008200004100

3F3 [N]

500315010002500150011000

2F2 [N]

500001500

100001000815000200

1F1 [N]

300001Fig. 11234time [s]56708

筋電データ，筋長・筋伸縮速度を考慮した筋張力推定値

Fig. 12

従来の条件での筋張力推定値（筋電デー

タ，筋長・筋伸縮速度を考慮せず）• 関節間力や筋消費エネルギーの評価付近）で，その減少の度合が顕著であることが分かる．他の筋に関しては提案手法と従来法には大• ローイング運動以外での筋張力推定

きな相違は見られなかった．

• より正確な筋のモデルの導入

5.

おわりに参考文献

ローイング運動を例として，実験で得られた筋電データを反映可能な筋張力推定法を提案した．

1) Pedotti, A., Krishnan, V. V. and Stark, L.:

筋長と筋伸縮速度も同時に考慮することが可能でOptimization of Muscle Force Sequencing inあり，従来法と比べてより妥当な推定ができるとHuman Locomotion, Mathematical Biome-期待できる．

chanics, 38, 57/76 (1978).

今後の課題として，以下を挙げる．

2) 宮脇和人，巖見武裕，大日方五郎，近藤康夫，• 被験者数を増やして，提案手法の妥当性を確沓澤圭一，小笠原雄二，西村信三：電動歩行機を用いた高齢者歩行の評価（第 1 報，平地認すること–8–歩行），日本機械学会論文集（C 編），65-640,4759/4766 (1999).
3) 山藤和也，巖見武裕，宮脇和人，平元和彦，島田洋一，大日方五郎：FES ローイング・マシンの開発と身体負荷の検討，福祉工学シンポジウム 2007 予稿集 (2007).

4) 阿江通良，湯海鵬，横井孝志：日本人アス

リートの身体部分慣性特性の推定，バイオメ
カニズム，11, 23/33 (1992).

5) Skelton, R.E., Iwasaki, T. and Grigoriadis,K.M.: A Uniﬁed Algebraic Approach to Linear
Control Design, 255, Taylor & Francis, London(1998).
6) Gahinet, P., Nemirovski, A., Laub, A.J. andChilali, M.: LMI Control Toolbox for Use
with MATLABr , The Mathworks, Inc., Natick, MA (1995).
7) 山本澄子，赤沢堅三，六馬信之，江原義弘，Pandy, M.G.: 合同セミナー 第四回 筋の力学モデルと筋張力推定テキスト，日本バイオメカニクス学会，日本臨床バイオメカニクス学会，バイオメカニズム学会臨床歩行分析懇談会，33/57 (1994).

8) 長谷和徳：身体運動モデル解析システム仕様
書 Ver. 2.0, 41/42 (1995).

9) Hatze, H.: A myocybernetic control modelof skeletal muscle, Biological Cybernetics, 25,
103/119 (1977).

–9–