ヒトの筋の粘弾性変化の推定 : モデルによるシミュレーション

愛 知 工 業 大 学 研 究 報 告 第2 6号B，平成3年

論 文

1

3

ヒトの筋の粘弾性変化の推定

ー モ デ ル に よ る シ ミ ュ レ ー シ ョ ン ー

A P U

l

p u

e u

H U HMm " H H 向 4 M W 側 聞 脚 " H H u ' π u n H U n v v -' 蜘 H U 悦 H H ﹄ TEL'AHU

・

l

a r T E A 4 1 L W A H V 白 川 ω

l

w M

，

V H H u t T L 欄 血

・l

・l

puFAU

-l

i

e o

白 u n d u

，

A H u t -A A H V n H U H M m

-n H V W 曹 J' A r u -h H H w q u . 曹 ' E 且 & 晶 伺 ‘ u " H M . ﹃ . 唱 lA 伺 K A u ' a411AU 曹 ' M

・

nHV'BEA n d u " H U 悦 H H A H v a n M U M -旬 a A -a + E E h ' n 凋 U . ， E E A T i n 4 H U W U ' 凸 U ' h H U 4 1 L W " H u n H H V 加 藤 厚 生 ¥ 伊 藤 晋 彦H Atsuo KATO+， Kunihiko lTO++

A

b

s

t

r

a

c

t

One approac主torea!Ize a co量p!I an t圃anIpu!ationrobot Is a I畢pedance

control. !n tlie IlUpedance contro! a d昌IUpIng and an e!astIcIt

v

.

coefIlcIents 01'

tlie robot jOInts are se!ected as a desIred va!ue. !ncIdentall.v， 品目lUanar量 see血S to acliIeve cOlUpliant量otionb.v contInuous regu!atIon 01'tlielIusc!e vIscos/ric -tiona! coel'f'jcient and elasticit

v

.

coel'l'icIent. To estiaate tlie vIsco-e!asticit_v cliange01'tlie IUUSC!e，曹ecarrIed0担tSO量豊 track ing哩，rperiaentson tlie liuaan曹rist jOInt and soae silUu!ationsI'or neuro-量uscu!ars

v

.

stesl aodel曹itlie!astic e自vi ron -lUen t. Tlien冒ecaae to a conc!世tion tliat tlie slusc!e vIsco-e!asticIt

v

.

increases 曹Itli tlie量usc!eactIvIt.y.事

e

tliInk tliat tlie d

v

.

nalUIcs 01'tlie coapl iant aanIpu!a -t ion roむot IUUS tむecontrolled b

v

.

tlie regu!ation01'tlieir joint visco-e!asticit

v

.

co員tinuous!

.

v

.

1 はじめに 工場でロボットは，ヒトに代わって危険な作業や 単調な繰り返し作業を行ない，その限りにおいて有 用である. しかしロボットは，塗装やアーク溶接の ように作業対象から離れて行う作業や，部材の搬送 のように対象が固定されない作業は巧みに行なうが クランクを回したり対象物の表面をなぞるような作 業対象との柔らかな接触を必要とする作業は苦手で ある.こうした作業は柔軟接触作業と呼ばれ，ロボ ット研究における今日の主要な課題の一つになって いる. 言うまでもなく， ヒトは接触作業を苦もなく実行 できる. ヒトと類似した形態をとるにもかかわらず ロボットが接触作業を苦手とする最大の理由は制御 + _{愛 知 工 業 大 学 電 子 工 学 科 ( 豊 田 市 〉} ++ 愛知県心身障害者コロニー(春日井市) 方式の違いにある.いま，工場で大量に使われてい るロボットの大部分は位置制御ロボットである. 位置制御されたロボットの剛性の高い手先を，硬い 環境(作業対象)との接触を保ちながら移動する場 合を考えてみよう.ロボットは指示された目標軌道 を正確にたどろうとし，環境も自分の位置を譲らな いから，環境の僅かな凹凸やロボットの軌道の片寄 りによって，接触が途切れたり環境を削り取ったり する不都合な現象が発生する. そこで，環境を傷っけない柔軟な接触作業を実現 するためには，位置だけでなくロボットの手先が環 境から受ける反作用力(相互作用力)も検出し，位 置と組み合わせて制御しなければならないけ. これは一般に力制御と言われる. 力制御法の一つにインピーダンス制御がある2) インピーダンス制御は，ロボットの手先の運動と手 先が環境から受ける相互作用力の関係を調節して安 定な接触作業を実現しようとする.運動と相互作用

14 愛 知 士 業 大 学 研 究 報 告 ， 第2 6号 B，平成 3年， Vol.26-B. Mar.1991 力の関係はインピ ダンスと言われ，具体的にはロ ボットの手先から見た↑貫性と粘性および弾性である. 手先から見たインピーダンスは，関節トルクと関節 粘弾性に変換できる3)ので，実際にはこれらを調節 することになる. ヒトは3 制御系として見たとき関節構造や筋群に 見られる複雑で冗長な構造や，神経系の信号伝播時 間に由来する不安定性をもっ上肢を自在に制御して， 自由空間でも環境の拘束を受けた空間でも，柔軟で 安定な運動を見事に実現している.その運動を遂行 するメカニスムは階層構造をなし，制御系としてみ れば上位の視覚による位置制御系と，下位の体制感 覚による運動制御系に大別される.運動の目標値は インパルス系列の形で中枢から出力され9 骨格筋へ 送られて筋張力を調節する一方，規範モデルとして 筋紡錘へフィードフォワードされる4)規範モデル から運動が逸脱すると脊髄レベルの伸張反射弓と縫 反射ループによるフィ ドパック制御によって修正 されるが，このとき神経系の信号伝播遅れ等によっ て発生する系の不安定性は筋の粘弾性によって抑制 される5) すなわち，筋は活動レベルを変化してそ の張力を変化するが，活動レベルの変化は同時に筋 粘弾性係数をも変化するため，神経筋系は運動の発 現に先立って系の安定性を予め高めておくことがで きる また，互いに措抗する筋が同時活動するとき 張力については差，粘様性係数については和の効果 をもたらすので引6) 関節トルクと関節粘弾性をほ ぼ独立に調節することもできる. ところで，任意の筋活動レベルにおける筋の粘弾 性係数については，摘出筋では知られている7)が， 生体筋では静的な等尺性収縮から微少な運動を発現 させた場合について知られているにすぎない. 筆者らはすでに，弾性負荷に拘束された手関節の 連続的な等速度屈曲運動に現れる振動周波数の推移 を実験によって確かめ，その原因を筋の粘弾性係数 変化にあるとしたB) しかしこの変化に粘性と弾性 の変化がどの様に関与するかを明らかにするには至 らなかった.本稿ではモデルによる計算機シミュレ ーションによって筋の粘弾性変化を考察する. 2 拘束運動する上肢のシミュレーション 2. 1 モテール 環 境 に 拘 束 さ れ て 運 動 す る 神 経 筋系のモデルを 図1のように定めた圃 モデルにおいて，主動筋の能動部は図の左半部に示 す直並列粘弾性モデルで表されるとした9) 筋は 一般的な2要素モデルにしたがって力発生部と直列 弾性要素 1

，

からなるとし.力発生部は力発生要素 f; と 並 列 弾 性 要 素 ん ， な ら び に 並 列 粘 性 要 素 ん で構成した.力発生要素と並列弾性要素は常に共同 して働き，次式の力を発生する。 f = f; ん ( 1 ) 各要素の粘弾性係数は筋活動によって変化する園 実験では弾性負荷を使ったためもっぱら主動筋のみ が活動し3 括抗筋の活動はほとんど認められなかっ た.従ってここでも詰抗筋の能動活動による粘弾性 要素を無視した園また，実験は主動筋が自然長から 収縮する方向に行なったので，主動筋の受動部の影 響は無いとして無視し，括抗筋の受動部はモデルに 図示しなかったが単純な並列粘弾性要素とし式に加 えて計算した固 (後述，式(7)， (11)) 図1 神 経 筋 環 境 系 の モ デ ル 筋紡錘は筋の伸縮長と伸縮速度を検出する能力を 持っとする.赤沢ら11)は筋紡錘について複雑な直並 列粘弾性モデルを提唱しているがここでは簡単なモ デルを採用した.入出力方程式は次式となる u， = g，(Xm+んdx./dt) 2 ( ) ただし lJS:筋紡錘出力 ん:筋の長さ g， ゲイン 九:速度検出ゲイン 筋紡錘の各パラメータは，

r

神経系による錘内筋の 制御を受けておおむね一定値に保たれるとした.

ヒ 卜 の 筋 の 粘 弾 性 変 化 の 推 定

1

建器官の特性は式 (3)で表わされるとした固 U， = g，[， ただし U

，

鍵 器 官 出 力 g

，

力検出ゲイン 神 経 系 の 伝 播 時 間 に よ る 遅 れ は 一 定 値 を と る た め 振 動 周 波 数 の 変 化 に は 影 響 を 与 え な い も の と し ， モ デルには加えなかった. モ デ ル の 運 動 方 程 式 は 式(4)となる. [;= C(1汁U，-U，) =sm dXm/ dt+ ImXm+ [， [，=/{，(Xm-x)=J!md2 X/ dt2+ [ [ =JI

，

d2 X/ dt2+ s

，

dx/ dt+/{，x これに式(l) ~(3) を加えラプラス変換すると，フロ ッ ク 線 図 図2と ， 伝 達 関 数 式(5)を得る. 図2 神 経 筋 環 境 系 モ デ ル の ブ ロ ッ ク 線 図 X cK， Cm (S)= - =一一一一一一 1 1 n3S3+ n2S2十nlstno ) 5 ( ただし， Jf=ん +J!， として， n3=Jlsm n2竺 JI(ん+1，(j+cg，))+sms，

n

1= s， (ん+l，(l十cg，))+ん (l，+ん) -cg， T， l， no= Iml，+(J+cg，) l，l，+l，lm-Cl，g， 次に，視覚を含む位置フィードパ、 yクループを加 え，図3の位置制御系を構成する.現実に生体系て、 は， ラ ン プ 入 力 に 対 す る 定 常 位 置 偏 差 が 有 限 値 に お さまることから，積分器を 1個加えた. 伝 達 関 数 は 式 (6)となる. X 1，/η (S)=一 一 = 一一一一一二 Xd H4S4+H3S3+H2S2+/flstHo ) n h υ (

1

5

) 3 ( 7こ7こし ，

H

4 = 113 H3= n2 H2 = nl 111 = nO 110= AcK，/[1 (4) 図3 位置制御系のプロ yク線図 2 _ 2 シミュレーション 計 算 機 シ ミ ュ レ ー シ ョ ン を 行 う に あ た り ， 筋 の 粘 弾 性 係 数 は 次 の よ う に 変 化 す る と し たa 先 ず 弾 性 特 性 に つ い て 図 4の筋長 収 縮 力 特 性 を 考 え る 司 図 は 筋 が 互 い に 括 抗 す る 場 合 を 模 式 的lこ線形 化 し て 描 い で あ る . ん は 筋 の 自 然 長 ， ん は 主 動 筋 の 能 動 要 素 の 最 大 収 縮 力 特 性 で あ る . 筋 の 最 大 収 縮 長 は 一 定 値 ん を と る た め ， 収 縮 力 特 性 曲 線 が 横 軸 と交わる位置は固定される圃 したがって，図から， 弾 性 負 荷 に 拘 束 さ れ た 筋 が 措 抗 す る 筋 の 協 同 活 動 を 得 な い で 筋 長 ま た は 筋 力 を 変 え る に は 自 ら の 弾 性 係 数 ん を 変 え な け れ ば な ら ず ， そ の た め に 自 身 の 活 動レベルを変化しなければならない. このことを め，A24a，と筋活動レベルを!順次低下さた場合につ いて示した， pは 受 動 要 素 の 特 性 で あ る 固 括 抗 筋 の 特 性 は ダ ッ シ ュ を 付 け て A' と p'で 示 し た . 図 中 右 下 が り の 直 線 ん は 弾 性 負 荷 の 負 荷 直 線 で あ るE 弾 性 負 荷 の 平 衡 位 置 と 筋 が 自 然 長 に な る 位 置 は 一 致 さ せた 筋 は 任 意 の 筋 長 で 負 荷 と つ り あ う 筋 力 を 発 生 で き る が ， そ の と き 筋 長 一 筋 力 関 係 は 常 に 負 荷 直 線 上 に な け れ ば な ら な い . 筋 の 弾 性 係 数 は 負 荷 弾 性 力 と筋弾性力が平衡する条件から次式によって決まる目 ん = ( f，+p )/( _ X，-X )+κ ( ) 7 た だ し ， ん : 筋 の 弾 性 係 数 f {， :負 荷 の 弾 性 係 数 X 弾性負荷の伸縮長さ(=筋収縮長) Xc:筋の最大収縮長さ κ:そ の 他 の 残 留 す る 弾 性 係 数

ただし， β :受動要素の粘性係数 α は定数で，シミュレーション結果に適度な振動が 見られるように選定した. ことができる.このときt AJ

，

んのように線形であ れば，式(7)で求めたんと最大収縮長 %cから筋活 動 レ ベ ル は ん%cで与えられる.そこで，筋の並列 粘性係数が筋活動レベルに依存すると仮定すれば次 式を得る. ) - E ム 1 A ( 2.3結果および考察 図6(a)， (b)， (C)は ら{s}の特性根の根軌跡， 図7は f!{s}の特性根の根軌跡である.いずれも， ん

，

l

，

お よ び / ま た は ん を 式(7)-(12)に従って 変化して描いた.選定した条件下では代表的な2根 para冊eters unit M.= 0.0009 Mm = 0.001 kg-m' B.= 0 Bm= 0.8-1.2 Nms/rad K

，

= 4 Km= 1-3.2 Nm/rad κ = 1 Kt= 5-16.2 Nm/rad g.= 1 gT= 1 T.=l.l T1 = 1/100 sec A = 2 C = 1 α = O.2/9 β = 0.8 シミュレーションに使ったパラメータ 表1に選定した定数および変数の範囲を示す. 100% 筋の粘性係数変化 愛知工業大学研究報告，第26号B，平成3年， Vol.26-B. Mar.1991 50

v

e

l

o

c

i

t

y

B.=αl.%c+β 図5

。

Q

2

50 Q3 100% 表1 む U ﹂ O ﹄ p

，

ま括抗筋の受動部の弾性特性であるが，簡単のた め，措抗筋の受動部の特性は並列弾性要素で表され， 張力Pが指数関数になり，更に自然長の60%伸展時 に自然長に於ける最大筋力 fjJD"xの約2倍なると仮 定して， とし

T

こ. ところで， トラッキング実験は 10deg/sの等速度屈 曲運動として行なったから，手関節の随意最大屈曲 角を 90deg とすればそこに歪る時間は 9sec を要 することから， ) o n u ( ) n u υ ( 筋の能動部分の弾性係数は直列弾性要素の方が大き く，並列弾性要素の5倍になる11)とした. ) n 川 ν 4 E A ( 次に筋の粘性特性について考える.図5は，摘出 筋から得た筋の収縮速度一収縮力特性刊を生体筋に 適用し模式的に再現したものである.右下がりの曲 線群は筋の粘性特性を表し，その傾きは粘性係数で ある.実曲線は筋の最大収縮特性を示し，Qoは静止 時最大収縮力である.ここでは静止時収縮力は筋活 動レベルの低下と共に Q/，Q2，Q3と低下し，同時に 筋の粘性係数も低下するとした. 図に示した垂直線分 t-t'は定速度で伸縮する弾性 負荷の負荷直線である.この負荷に拘束された筋が 負荷トルクの増減に対応して活動レベルを変えると その粘性係数も変化する. ところで，筋の活動レベルは，図4に示した筋力 特性曲線 Ao

，

A/

，

A2..が縦軸を横切る値で代表する 筋の活動レベルと弾性係数変化 p = erp( O. 3%)ー/ % = ( .rc/9) t

l

r

=

5l

旭 図4 16

1

7

I

m

2 与O ヒトの筋の粘弾性変化の推定 は何れも複素根になった 図6は位置フイ ドパ yク を か け な い 神 経 筋 系 の 正味の特性を表し， (a)が筋の弾性係数のみ変化し た場合， (b)は粘性係数のみ変化した場合， (c)は 七

- 一

ハ_u 回 部 V A l x x u λ V ぺ

2

。

ー社

real roots are aut of range ー持O ヌ ズ X

-

2

粘性係数と弾性係数を共 iこ変化した場合を示す. 図から弾性係数変化は周波数変化に大きな影響を与 えず，粘↑生係数変化は大きな影響を及ぼすことが判 る.図 7は視覚系を含む位置フィードパックのある 場合を表しているが，追加した積分器の効果を除け x x )()()(ぷ ば基本的な特性は図 6と同様になる. 位置制御系モデルの特性根軌跡 両図共に，

i

也の 1根または 2根は大きな負の実部を もち図示した範囲には入っていないため，振動現象 に与える影響は小さい. 図8は，実験に合わせて行ったランプ応答のシミ ュレーション結果の一例である. 計算区間を O.02secおきにとり，区間線形として逐 次高速逆ラプラス変換法を適用し算出した. 30T (deg) 図 7 Re

d

-40 向

/

〆

a.>20 CJl c (tJ E

-

x a.> ご10 ( a )筋の弾性係数を変化した場合 tlm巴 (sec) 位置制御系モテ'ルのランプ応答 2 Re -4 -2

real toots are

Out of range 一与0，

j

区X 叉 " )( " x >1<"><" x ~

。

図 8 図 9は，図 8から A Rモデル推定法によってパワ ースベクトルを計算した結果である. ( b )筋の粘性係数を変化した場合

バ 童 日 一

﹂ 由 z o a Re 2 -4 -2

real roots are out of range 一 ー井口， .. • ， . ヌ ー 鴨 . コ 巳 ， 支 援 1 同

。

<.J 0 15(Hz) ランプ応答のパワースベクトル 図9 ( c )筋の粘弾性係数を変化した場合 神 経 筋 環 境 系 の 特 性 根 軌 跡 図 6

1

8

愛知工業大学研究報告，第26号B，平成 3年，

V

o

l

.

2

6

-

B

，

M

a

r

.

1

9

9

1

図 10 ，こパワースペク卜ルがピークをとる周波数 を負荷トルクに対して描いた. kは ん

，

1

，

を， b は ん を ，

b

-

k

は ん，lt，んを変化した場合を示す. 根軌跡と同様，振動周波数の増加は筋の粘性係数の 増加に強く依存し，弾性係数にはほとんど依存しな いことが分かる.ここでは筋粘弾性係数を式 (7)~

(

1

1

)

により筋活動レベルに従って変化した.前述し たトラッキング実験の結果とよく一致していること が判る. (Hz) 7

M P

J

6 ﹁ ﹁ J L H H ， h U E U コ ロ 山 ﹂ U F 五回山口

。

2

l

o

a

d

e

d

t

o

r

q

u

e

(

N

m

)

図 10 ピーク周波数の推移 4. 結 論 先に，弾性負荷を付けた手関節で行った等速度ト ラッキング実験では， トラッキング偏差に現われる 振動のうち 5~10Hz の成分が負荷トルクに強い線形 相関を示して高域へシフトする現象を確認した. モデルによるシミュレーションからも，負荷トル クの増大と共に振動周波数を増加する結果を得た. 振動周波数の増加は筋の粘性係数の増加に強く依存 し，粘性係数と弾性係数を同時に増加したときも同 様の傾向を示した.弾性係数のみ増加したときは周 波数変化はほとんど認められなかった.筋の粘弾性 係数を筋活動レベルに合わせて増加したところ，実 験結果によく適合した. 従来，弾性負荷を付けた屈曲運動では，筋活動レ ベルの増加とともに筋弾性係数が増加すると考えら れている.われわれは，このとき同時に粘性係数も 増加していると結論付けた.筋粘性の増加は一般に 制御系の安定性向上に寄与する.このことは実験で も確認し，モテソレの式(6)(7)の根の変化からも確認 している. 本稿では，環境の特性を回定した単純な等速度運 動中にも，ヒトの上肢が筋粘弾性係数を連続的に変 化していることを明かにした.いま，接触作業ロボ ットの制御では，環境に応じて選定した目標関節粘 弾性を実現するようにアクチュエータトルクを調節 して行っているが，環境の特性が変わらない限り粘 弾性係数を連続的に変化する制御法は採用されてい ない.われわれが得た結論からは，環境の特性が変 化しない場合も関節粘弾性係数をアクチュエータト ルクに応じて連続的に調節するのが自然である.そ れがロボットの柔軟接触問題を解決へ向けて前進さ せる手がかりになると確信している. 参考文献

1

)

M

.

T

.

M

a

s

o

n

:

C

o

m

p

l

i

a

n

c

e

a

n

d

F

o

r

c

e

C

o

n

t

r

o

l

f

o

r

C

o

m

p

u

t

e

r

C

o

n

t

r

o

l

l

e

d

M

a

n

i

p

u

l

a

t

o

r

s

，

l

E

E

E

S

M

C

-

1

1

-

6

，

4

1

8

-

4

3

2

(

1

9

8

1)

2

)

N

.

H

o

g

a

n

I

m

p

e

d

a

n

c

e

c

o

n

t

r

o

l

:

A

n

a

p

p

r

o

a

c

h

t

u

m

a

n

i

p

u

l

a

t

i

o

n

:

1

I

I

m

，

A

S

M

E

D

S

M

C

，

1

0

7

，

1

-

2

3

(1

9

8

5

)

3

)

J

.

K

.

S

a

l

i

s

b

u

r

y

:

A

c

t

i

v

e

S

t

i

f

f

n

e

s

s

C

o

n

t

r

o

l

o

f

a

M

a

n

i

p

u

l

a

t

o

r

i

n

C

a

r

t

e

s

i

a

n

C

o

o

r

d

i

n

a

t

e

s

，

P

r

o

c

.

1

9

t

h

I

E

E

E

C

D

C

，

9

5

-

1

0

0

(1

9

8

0

)

4

)

伊藤宏司:筋運動制御機構，計測と制御，

2

5

，

2

，

2

3

-

2

7

(1

9

8

6

)

5)伊藤宏司，辻敏夫:筋骨格系の双線形特性と義肢 制御への応用，電気学会論文誌，

1

0

5

C

，

2

0

2

-

2

0

8

(

1

9

8

5

)

6

)

N

.

H

o

g

a

n

A

d

a

p

t

i

v

e

c

o

n

t

r

o

l

o

f

m

e

c

h

a

n

i

c

a

l

i

m

p

e

d

a

n

c

e

b

y

c

o

a

c

t

i

v

a

t

i

o

n

o

f

a

n

t

a

g

o

n

i

s

t

r

n

u

s

c

l

e

s

，

I

E

E

E

A

C

-

2

9

，

8

，

6

8

1

-

6

9

0

(1

9

8

4

)

7

)

H

.

M

a

s

h

i

r

n

a

，

K

.

A

k

a

z

a

w

a

，

a

n

d

H

.

K

u

s

h

i

r

n

a

T

h

e

f

o

r

c

e

-

l

o

a

d

-

v

e

l

o

c

i

t

y

r

e

l

a

t

i

o

n

a

n

d

t

h

e

v

i

s

-c

o

u

s

l

i

k

e

f

o

r

c

e

i

n

t

h

e

f

r

o

g

s

k

e

l

e

t

a

l

r

n

u

s

c

l

e

，

J

a

p

a

n

e

s

e

1.

o

f

P

h

y

s

i

o

l.，

2

2

，

1

0

3

-

1

2

0

(1

9

7

2

)

8)加藤厚生，伊藤晋彦:ヒトの筋の粘弾性変化の推 定，愛工大研究報告，

V

o

l.

2

5

8

，

1

9

-

2

6

(1

9

9

0

)

9) 赤沢堅造，藤井克彦，真島英信:骨格筋の力学特 性とそのモデルー収縮力への依存性 ，医用電子 と生体工学，

1

7

，

4

，

3

0

-

3

6

(

1

9

7

9

)

1

0

)

R

.

B

.

S

t

e

i

n

P

e

r

i

p

h

e

r

a

l

C

o

n

t

r

o

l

o

f

M

o

v

e

r

n

e

n

t

，

P

h

y

s

i

o

l

o

g

i

c

a

l

R

e

v

i

e

w

s

5

4

-

，1

2

1

5

-

2

4

3

(

1

9

7

4

)

1

1

)

赤沢堅造，楠本秀忠，藤井克彦ヒト骨格筋の収 縮力学定数の推定法¥医用電子と生体工学，

2

6

，

4

，

1

8

-

2

6

(

1

9

8

8

)

(受理平成3年 3月20日)