平成 13 年度
修士論文

上肢リハビリテーションロボットの制御方法の検討（制御法の提案と筋張力シミュレータの有効性の確認）

指導教員
副指導教員

井上
八田

喜雄
章光

教授
助教授

甲斐

高知工科大学 大学院工学研究科
基盤工学専攻 博士課程（前期）1045008 伊東 禅

義弘

助手

目次
１． 緒言・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・１
1.1 関節リハビリ機器の現状・・・・・・・・・・・・・・・・・・・・・・・・１
1.2 本研究の位置付けと研究目的・・・・・・・・・・・・・・・・・・・・・・１
1.3 本研究の内容・・・・・・・・・・・・・・・・・・・・・・・・・・・・・２

２．筋張力シミュレータ・・・・・・・・・・・・・・・・・・・・・・・・・・・３
2.1

2.2

上肢モデル・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・３
2.1.1

筋張力と腕先力の関係・・・・・・・・・・・・・・・・・・・・・･３

2.1.2

筋肉の最大張力・・・・・・・・・・・・・・・・・・・・・・・・･４

2.1.3

筋肉の歪み・・・・・・・・・・・・・・・・・・・・・・・・・・･４シミュレータ・・・・・・・・・・・・・・・・・・・・・・・・・・・・・4
2.2.1

双対法の概要・・・・・・・・・・・・・・・・・・・・・・・・・･５

2.2.2

目的関数及び，制約条件・・・・・・・・・・・・・・・・・・・・･５

2.2.3

固定パラメータ・・・・・・・・・・・・・・・・・・・・・・・・･６

2.2.4

シミュレータの処理手順・・・・・・・・・・・・・・・・・・・・･７

３．ロボットアームの制御法・・・・・・・・・・・・・・・・・・・・・・・・・８
3.1

3.2

3.3

仮想コンプライアンス制御・・・・・・・・・・・・・・・・・・・・・・・８
3.1.1

基礎制御式・・・・・・・・・・・・・・・・・・・・・・・・・・･８

3.1.2

仮想パラメータ・・・・・・・・・・・・・・・・・・・・・・・・･８

仮想コンプライアンス制御の改良・・・・・・・・・・・・・・・・・・・・９
3.2.1

制御法の考案とシミュレーション・・・・・・・・・・・・・・・・･９

3.2.2

シミュレーション結果とその考察・・・・・・・・・・・・・・・･１１関節角度推定プログラム・・・・・・・・・・・・・・・・・・・・・・・１２

４．実験・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・１４
4.1

4.2

筋張力シミュレータ実験・・・・・・・・・・・・・・・・・・・・・・・１４
4.1.1

実験装置・・・・・・・・・・・・・・・・・・・・・・・・・・･１４

4.1.2

実験方法・・・・・・・・・・・・・・・・・・・・・・・・・・･１４

4.1.3

実験結果・・・・・・・・・・・・・・・・・・・・・・・・・・･１５ロボットアーム制御実験・・・・・・・・・・・・・・・・・・・・・・・１５
4.2.1

実験装置・・・・・・・・・・・・・・・・・・・・・・・・・・･１５

4.2.2

実験方法・・・・・・・・・・・・・・・・・・・・・・・・・・･１６

4.2.3

実験結果・・・・・・・・・・・・・・・・・・・・・・・・・・･１６

５．考察・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・１７
5.1

筋張力シミュレータに関する考察・・・・・・・・・・・・・・・・・・・１７

5.2

ロボットアームの制御方法に関する考察・・・・・・・・・・・・・・・・１８

６．結言・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・２０
謝辞・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･２０
参考文献・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・２１
図表・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･２３
2 章図表・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･２３
3 章図表・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･２８
4 章図表・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･３５5 章図表・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･６５

付録・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・･６６
付録図表・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・・７３

1 緒言
1.1 関節リハビリ機器の現状
関節に障害を持ち，そのリハビリテーション（以下，リハビリ）を必要とする人は，交通事故や老人性の運動機能障害者を含めて，数十万人存在すると言われている．日本の高齢化の現状を加味すると，この数は今後更に増加するものと考えられる．一方，リハビリを実施する側である理学療法士（以下，療法士）を例にとると，100 床につき 1 人の割合(1)の配置でしかなく，これでは個々の患者にきめ細かい治療を長時間行うことは困難である．
この負担を軽減し，患者個々人に合わせたリハビリを実施するために，自動化された機器が必要とされている．現在，実用化されている機器に CPM（Continuous Passive Motion）
装置が有る．これは持続的に低速度の他動運動を繰り返す装置で，関節癒着防止などの効果が期待できる．しかし，CPM は等速で動く 1 自由度しか持たない単純構造で，運動の柔らかさを持たせる工夫も装置に遊びを持たせる程度に留まっている．
このような CPM の欠点を改善し，療法士の徒手による理学療法の柔らかさを，機械で再現させようとする研究が行われている．下肢に関しては膝用の多自由度型 CPM を安川電機が作成し (2) ，それに岡島ら (3) がインピーダンス制御を導入して，関節可動域（Range of
Motion 以下，ROM）改善のみならず，筋緊張抑制効果も持つ TEM（Therapeutic Exercise
Machine）を考案している．上肢のリハビリに関しては CPM では無いが，則次ら（4）のゴ
ム人工筋ロボットを用いて各種の運動訓練モードを実現させた研究や，坂口ら（5）の ER（電気粘性流体，Electrorheological Fluid）アクチュエータをリハビリ訓練システムに応用した研究がある．

1.2 本研究の位置付けと研究目的
先述した訓練機器制御に関する研究の多くは，インピーダンス制御を用いている．平林（6）

らが発表した仮想コンプライアンス制御も広義ではインピーダンス制御の一種であり，リハビリ等の人間と機械を協調させるための制御においては，インピーダンス制御が多用されている．また，人間の動きも機械インピーダンスで表現できる(7)ため，人間のインピーダンスを計測する研究（8）なども行われている．
インピーダンス制御は協調作業において有効な制御法である．しかし ROM 訓練においては問題が生じる場合もある．例えば，脳卒中などによる神経系の損傷に起因する痙縮（Appendix 用語説明参照）の上肢 ROM 訓練では，関節角速度を 3〜5[deg/sec]に抑えた
他動運動による関節の屈曲・伸展を行う．痙縮では筋が他動的に伸展すると，それに伴い抵抗が大きくなる．この抵抗は伸展速度が速ければ，比例して発生が早まる．また，ある程度伸展した時点で，筋はいわゆるジャックナイフ

反応で突然弛緩する．痙縮がこのような特性を持つため，現在の徒手によるリハビリでは，療法士が患者の負担にならない
ように速度を調整しているのである．
インピーダンス制御で他動運動を行う機器の移動量は，抵抗などの力と，目標位置と現1

在位置との偏差，この 2 つの要素から導かれる．大きな抵抗が加わると，力の釣り合いから，位置偏差が拡大する．拡大した位置偏差は，抵抗減少時に減少し，その結果移動量が増加する．痙縮等の症状がある場合は，一定以下の速度で動作することが望ましい．移動量（速度）の急激な増加は患者に負担を与える可能性が有り，適切な他動運動を実施することは困難である．よって，このような症状に対応した制御法が必要となる．
また，上記の研究は手先軌道の生成や，インピーダンスパラメータの決定などを主に行っており，患者の筋の活動状態については考慮していない．リハビリにおいて患者の筋の活動状態は，その進行具合をはかる指標の一つであり，リハビリの内容を査定するための判断基準の一つでも有る．これは機器を用いたリハビリにおいても軌道決定等で考慮すべきである．しかし，筋の活動状態を計測する，例えば筋電位計を用いるには，準備に時間が掛かり，運動療法を実施するたびに計測することは現実的でない．さらに，筋電位計を装着することで，患者に心理的，身体的ストレスを与えることも考えられる．
そこで，力センサを用いて患者の腕先に発生する力（以下，腕先力）を計測し，そこからシミュレーションにより筋張力を算出することで，筋の活動状態を推定できると考えた．
そのデータを用いて，運動療法実施前に，負荷・軌道の決定等を行い，患者個々人にとって最適な機器の制御を行うことも可能であろう．また，運動中の筋張力に関するシミュレ
ーション結果を解析することで，進行具合の確認や，使用するパラメータの設定等も効率良く行えるものと考えられる．さらに，患者に対するフィードバックとして，シミュレー
トしたデータを見てもらうことで，ストレスを与えること無くリハビリの効果を確認してもらい，リハビリに対する意欲の保持・増大も望める．
以上より，患者の筋活動を考慮し，痙縮等の症状にも対応した，上肢リハビリ機器の実現を最終的な目標とし，本研究では，筋活動を考慮した制御を行うために必要な筋張力シ
ミュレータと，痙縮に対応した制御法の考案及び，それぞれの有効性を確認することを研究目的とする．

1.3 本研究の内容
本研究では，2 章で上肢モデルの構築とそれを用いた筋張力シミュレータについて説明する．3 章ではロボットの制御手法について述べる．4 章では制御手法とシミュレータ，それぞれの有効性を確認するための実験方法について説明し，5 章でそれらの結果について考察を加え．最後に 6 章で 2 章から 5 章までを簡潔に要約し，結言を述べる．

2

2 筋張力シミュレータ
本章では，水平面運動における腕の筋張力シミュレータについて説明する．
作成したシミュレータは，腕の 2 自由度水平面運動における筋張力，関節角度及び，腕先力の関係をモデル化し，そこに腕先力と関節角度を入力して，2 次計画問題の解法の一つである双対法により，エネルギー最小となる各筋の筋張力を推定するものである．

2.1 上肢モデル
本研究では，岩本(9)が構築した 2 自由度水平面運動学モデルを参考として，右上肢のモデルを作成した．座標等を Fig.2.1.1 に示す．
腕の水平面運動とは，肩関節の水平屈伸運動と，肘関節の屈伸運動からなる運動で，手根関節部の運動は含まない．この運動は岩本や Piegon(10)らのモデルによれば 8〜9 個の筋が関与している．この中の上腕三頭筋（Triceps Branchii 以下，TB）は機能が肩関節と肘関節を動かす二関節筋と，肘関節のみを動かす単関節筋に分かれている．内側頭と外側頭が単関節筋，長頭が二関節筋として働く．Piegon らのモデルでは TB を一つの筋とみなしているが，本研究では肩関節と肘関節の運動を組み合わせた水平面運動を取り扱うので，TB を機能によって 2 つに分ける．一つは上腕三頭筋長頭：TBLo，もう一つを上腕三頭筋
外側頭（内側頭を含む）：TBLa とし，モデル全体で 9 個の筋を扱う．Fig.2.1.2 に扱う筋の位置， Table2.1.1 に扱う筋の一覧を示す．
以下に，モデル作成の過程を説明する．

2.1.1 筋張力と腕先力の関係
人間は筋の張力を制御することで，関節モーメントを発生させ，そこから腕先に力を発生させる．つまり，各筋の合力が，腕先の力となる．上肢モデル作成に当って，まずは重要な制約条件である筋張力を入力とし，腕先力を出力とするこの関係を定式化する必要がある．微小関節角を，微小腱変位に変換するモーメントアーム行列を T とすると，仮想仕事の原理より，筋の張力ベクトル f と，関節モーメントベクトル ô の関係はoˆ = T T ⋅ f（2.1.1）となる．同様に，微小関節角を，微小手先位置に変換するヤコビ行列を J とすると，仮想仕事の原理より，関節モーメントベクトル ô と，腕先力ベクトル F の関係は( ) ⋅ oˆ−1

F = JT（2.1.2）となる．最終的に入力：筋張力ベクトル f と，出力：腕先力ベクトル F の関係は，式（2.1.1）と式（2.1.2）より( ) ⋅T ⋅ f

F = JT−1T（2.1.3）と示すことができる．（ J ， T の算出は Appendix A を参照）

3

2.1.2 筋の最大張力筋には発揮できる限界の張力（最大張力）があり，シミュレータ作成に当っては目的関数にも制約条件にも用いるこの限界値を考慮しなければならない．最大張力は計測の仕方
によって幾つかの種類があるが，本研究では等尺性最大張力を用いる．等尺性最大張力とは，筋を任意の長さに固定した状態で，筋が発揮できる最大張力である．2.1.1 項で述べた
筋張力と手先力の関係では，各筋の等尺性最大張力は，全て等しく評価される．しかし，実際には各筋が出しうる等尺性最大張力は，それぞれ異なっており，筋の生理断面積（以下， PCSA ）に比例する(11)．よって，ある筋の等尺性最大張力 f max n は

f max n = k ⋅ PCSAn（2.1.4）と示せる． k は PCSA 1[cm２]あたりの比例定数で，一般人で 5〜6，スポーツなどで筋を鍛えた人間では 10 程度になる．

2.1.3 筋の歪み
2.1.2 項で述べたように，筋は発揮し得る最大の力，等尺性最大張力を持っている．この等尺性最大張力は，関節角度（姿勢）によって変化する筋の歪み h の影響を受ける．歪みとは，ある角度における筋の変位量と，自然長となる角度における変位量の差を自然長で割ったものである．自然長とは，筋が伸展も収縮もしていない状態を指す．これは等尺性最大張力が，筋の収縮要素による能動的最大張力と，並列弾性要素による静止張力との和によって構成されていることによる(12)．
（Appendix 用語説明参照）
筋が自然長状態となるときの等尺性最大張力を 1 として，無次元化した等尺性最大張力と筋の歪みの関係と，能動的最大張力と静止張力の特性曲線を Fig.2.1.3 に示す．図より能動的最大張力は，筋の歪みが 0 のとき最大値をとり，筋が自然長より長くても，短くても，値が小さくなることがわかる．また，静止張力は，筋を自然長以上に伸ばすと発生し，非線型なバネの性質を持っている．式（2.1.4）より，筋の歪みを考慮した，ある関節角度 θ nにおける筋の等尺性最大張力は

f max n (θ n ) = fˆ ( h) n ⋅ k ⋅ PCSAn（2.1.5）（ fˆ ( h) n のとなる． fˆ ( h) n は，筋の歪みより求まる無次元化した等尺性最大張力である(13)．
算出は Appendix B を参照）

2.2 シミュレータ
2.1 節で構築した上肢モデルに，2 次計画問題の解法の一種である双対法を用いて腕先力から各筋の筋張力を推定する．ここでは，双対法についてその概要を説明する．続いて，目的関数と制約条件を示す．次にシミュレータ内で用いた固定パラメータについて説明し，

4

最後にシミュレータの処理手順について述べる．

2.2.1 双対法の概要(14)2 次計画問題の一般的表記は，次のようなものである．
目的関数： F (x ) = c x + 1 x Gx 
→ 最小
T

T

2

制約条件： a i x = biTaiT x ≥ bi

i=1,2,…,me（2.2.1）i=me+1,…,m

ここで， a i と c は n 次元列ベクトル， bi は実数， G は n×n 対象行列である．ただし，等式制約条件をまったく含まない問題に対しては me =0，不等式制約条件を含まない問題に対しては m= me と考える．双対法は上記問題に対して，独立有効集合 J s の列と，それに対応する J s ‐最適解の列を生成する方法であり，概要は次のようになる．
［1］ 適当な独立有効集合 J s に対する J s ‐最適解を計算し，それを初期解とする．

x が式（2.2.1）の実行可能解ならば， x は最適解であるので停止．そうでなければ，
点 x が満たしていない制約条件の添字 s ∉ J s を一つ選ぶ．
［3］ 部分問題 P (J s ∪ {s}) が実行不可能ならば，式（2.2.1）は実行不可能なので停止．

［2］

そうでなければ［4］へ．
［4］

J s′ ⊆ J を満たし，かつ独立であるような制約条件の添字集合 J s′ ∪ {s} に対する
(J s′ ∪ {s}) ‐最適解 x ′ で F (x′) > F (x )を満足するものをみつける． J := J ′ ∪ {s}，x := x ′ とおいて［2］へ戻る．

この手順を有限回繰り返して，最適解を求める．なお，双対法の詳しいフローチャートをFig.2.2.1 に示す．

2.2.2 目的関数及び，制約条件
2.2.1 項で述べたように，双対法を用いるにあたり，目的関数と制約条件が必要となる．
人間の腕先力に対する各筋の出力比を表現する数式は幾つか有るが，次に示す各筋の出力比の総和が最小となる数式が有効であると言われているので，本研究ではこれを採用した．
2

 fn 

 
→ 最小
∑
n =1  f max n i(2.2.2)

ここに f n は各筋の張力， f max n は各筋の等尺性最大張力である．本研究で取り扱う筋は 9個なので，式(2.2.2)内の i =9 となる．関節角度による筋張力の変化を考慮して，筋の歪みより求まる無次元化した等尺性最大張力 fˆ ( h) n を f max n に加え，式(2.2.1)，式(2.2.2)より目的5関数を次のように定めた．
2



2



0
 fˆ (h) BS ⋅ k ⋅ PCSABS 
1 T ⋅ f 
O
⋅f ⋅
→ 最小

2
2

 

2


0

 fˆ (h) ⋅ k ⋅ PCSA  PCPC  


(2.2.3)上 式 に お い て ， 添 字 BS 〜 PC は 各 筋 の パ ラ メ ー タ で あ る ． ま た

f = [ f BS

f BR

f AN

f TBLaf TBLo

f BB

f DA

f DP

f PC ] である．
T

次に制約条件が，3 種類 20 個定まる．まず，筋張力と腕先力の関係から，次式に示すｘ方向の腕先力 Fx と，ｙ方向の腕先力 Fy の等式制約条件が 2 つ導かれる．

(J ) ⋅ T ⋅ f =  F 
T −1

Fx

T

(2.2.4) y

つぎに，筋は 1 自由度のマニピュレータであり，収縮時しか力を発生しない．作成したモデルでは収縮を正としている．これらから各筋において，次式に示す不等式制約条件が導かれる．

fi ≥ 0(2.2.5)（i＝BS，BR，AN，TBLa，TBLo，BB，DA，DP，PC）
最後に，各筋張力はその最大張力を超えることは無いので，不等式制約条件− f i ≥ − fˆ (h) i ⋅ k ⋅ PCSAi(2.2.6)が導かれる．

2.2.3 固定パラメータ
上肢モデルを用いて腕先力から各筋の張力を推定するシミュレータには計算中に値が変化しない以下の固定パラメータが必要となる．
［1］

基準関節角度 θ sn ,θ en

［2］

PCSA

［3］

基準となる骨（前腕骨と上腕骨）長 ln と筋の長さ mln［4］筋の構造指標 ia［5］

PCSA 比例定数 k筋の歪みが 0 となる基準関節角度は，
肘関節については An(15)らの報告から θ en ＝78［deg］

6

の時に各筋は自然長になるとした．肩関節に関しては，関節可動域が-30[deg]から 135[deg]であるので，中間の θ sn ＝52.5[deg]を基準関節角度とする．

PCSA については複数の文献 (16) 〜 (19) のデータを平均化した数値を使用した．一覧を
Table2.2.1 に示す．
筋は両端の腱を介して骨に付着しており，筋の長さと骨長には相互に関連性がある．おおむね骨長に筋の長さは比例する．そこで基準となる骨（リンク）長 ln と，筋の長さ ml i を
岩本と Wood(16)の報告より求めた．筋の構造指標 i a は岩本の報告より引用した．それらの一
覧を Table2.2.2 から Table2.2.3 に示す．

PCSA 比例定数は，一般人で 5〜6，筋肉を鍛えた人間では 10 となる．数値は実験の被験者の状態に合わせて調整する．

2.2.4 シミュレータの処理手順作成したシミュレータの処理手順は以下に示す通り．
[1] 2.2.3 項で説明した固定パラメータを読み込む．
[2] 骨長から，各筋の自然長 ml i を計算する．
[3] 腕先力 Fx ， Fy ，関節角度 θ s ， θ e を読み込む[4] 関節角度から，筋の歪みを考慮した等尺性最大張力 f max n (θ n ) と T ， J を計算する．
（各々の計算は Appendix A と B を参照）
[5] 導いた各パラメータを双対法アルゴリズムに入力し，推定筋張力 f n を出力する．
なお，[2]で行う計算は，次のようなものである．
筋の長さ mli を，次式に示すように骨長 l i と基準骨長 l ni の比率（筋長補正係数）から計算する．

lmli = mlni ⋅ ilni(2.2.7)

ここに， mlni は基準となる筋の長さである．
シミュレータのフローチャートを Fig2.2.2 に示す．

7

3 ロボットアームの制御法
本章では痙縮上肢の運動療法に対応したロボットアームの制御法について説明する．はじめに基本となる仮想コンプライアンス制御と仮想パラメータについて述べる．次に痙縮に対応した，仮想コンプライアンス制御を改良した各制御法について説明する．最後に筋張力シミュレータで用いる人間の関節角度を，ロボットアーム先端位置から推定するプロ
グラムについて説明する．

3.1 仮想コンプライアンス制御
運動療法を実施する際には，患者の抵抗に対して柔軟な対応をすることが，安全上必要になる．前述したようにロボットでこれを実現するには力制御，なかでも様々な運動に対応できるインピーダンス制御が有効であり，多くの研究で使用されている．痙縮の他動運
動では速度を一定以下に抑える必要があるので，本研究ではインピーダンス制御の一種で速度を目標値に持つ，仮想コンプライアンス制御を基本の制御法として用いる．

3.1.1 基礎制御式
目的とする運動が水平 2 自由度なので，Fig.3.1.1 に示すような仮想バネ・マス・ダンパ
ー系をロボットアーム先端ｘ方向，ｙ方向に設定する．この時のアームの運動方程式は，次に示すようになる．
F = M ⋅ &x& + C ⋅ x& + K ⋅ x(3.1.1)

ここに， F ：アームに加わる外力， &x& ：アームの加速度， x& ：アームの速度， x ：アームの位置（以上，2 次元ベクトル）， M ：仮想質量， C ：仮想粘性定数， K ：仮想バネ定数（以上，2 行 2 列の行列）である．
式(3.1.1)を離散値系で表すために

dx& v n +1 − v n=dtdt(3.1.2)とする． dt ：サンプリング時間， v n +1 ：１サイクル先の速度， v n ：現在の速度である．式(3.1.1)に式(3.1.2)を代入すると，

 v − vn 
Fn = M  n+1
 + C ⋅ v n + K ⋅ dx n dt 

(3.1.3)dx n ：アームの現在位置と目標値の差である．式(3.1.3)を変形すると，
v n +1 =dt(Fn − C ⋅ v n − K ⋅ dxn ) + vn
M

(3.1.4)

このような仮想コンプライアンス制御の制御式が導かれる．

3.1.2 仮想パラメータ
3.1.1 項で説明した基礎制御式で使われている仮想パラメータについて説明する．ロボッ

8

トアームに療法士の行う他動運動のような柔軟な動きを与え，上肢の運動療法を円滑に行うため，水平方向の負荷に対する療法士のインピーダンス(3)を仮想パラメータとして取り入れた． C ＝90[N/(m/s)]， K ＝1,000[N/m]である．また，人間の腕の質量を計測し，仮想質
量 M ＝2.5[kg]と設定した．

3.2 仮想コンプライアンス制御の改良
3.1 節で説明した仮想コンプライアンス制御は，痙縮等の運動療法には適していないと考える．この節ではその理由を述べ，問題解決のために考案した仮想コンプライアンス制御
改良型制御法の説明及び，そのシミュレーションと考察を行う．

3.2.1 制御法の考案とシミュレーション
仮想コンプライアンス制御を用いて他動運動を行う場合，ある１サイクルの移動量は軌
道上を移動する目標位置の移動量と，力の影響から定まる．よって，目標位置移動量よりも抵抗力の影響が大きい場合，移動量は 0 以下となり前進することはできない．この状態
は目標位置と現在位置との偏差が大きくなり，仮想バネ定数と位置偏差の積 K ⋅ dx が増大
して，抵抗力 F を上回ると解消される．つまり，仮想バネの伸縮力が抵抗力を上回ることで現在位置の前進が始まる．なお，この時の速度は微小であるため，速度項の影響は無視できるほど小さい．
このように，仮想コンプライアンス制御による他動運動の場合，抵抗力との釣り合いをとるために，位置偏差が発生する．この偏差は抵抗力が減少すると共に減少するが，この時，パラメータの値によって大小はあるが，目標位置に追い付こうと速度が上昇する．痙縮のように抵抗力が急減する場合は，速度の上昇も大きい．このような速度の上昇は他動運動においては望まれない現象であり，痙縮に対する他動運動においては伸張反射を誘発
し，患者に負担を与える可能性もある．以下で，この問題を解決するために幾つかの方法を考案し，シミュレーションを行った．
シミュレーションは 1 自由度で行い，現在位置を 0[m]から 0.6[m]まで移動させる．サン
プリングタイムは 0.01[sec]，Fig.3.2.1 に示すような痙縮をイメージした力を加えた．また
速度に関しては，痙縮の他動運動が 3〜5[deg/sec]で行われる事と，4 章で述べる実験においてリンク長 0.4[m]の DD ロボットを用いる事から，目標位置は速度 0.021[m/sec]で等速に移動させ，限界速度 0.035[m/sec]とする．
・通常の仮想コンプライアンス制御…改良した制御法との比較のために，通常の仮想コン
プライアンス制御でもシミュレーションを行った．結果を Fig.3.2.2 に示す．
・速度リミッタ追加…一定以下の速度の運動では，痙縮は伸張反射を起こさない．よって
速度限界値を設定し，これを越えた速度をカットする速度リミッタを追加した．シミュ

9レーションの結果を Fig.3.2.3 に示す．なお，速度リミッタの設定は限界値である
0.035[m/sec]とした．
・パラメータ変更…3.1.2 項で述べた仮想パラメータは，療法士を用いた実験から得られたものである．しかし，Fig.3.2.1 のような力を与えた場合，Fig.3.2.2 に示すように抵抗急
減時に速度が上昇してしまう．そこで，抵抗急減時に速度超過が起きないパラメータの組み合わせの一例を計算で求め，シミュレーションを行った．用いたパラメータは M ＝1000[kg]， C ＝1800[N/(m/s)]， K ＝600[N/m]である．結果を Fig.3.2.4 に示す．
・目標位置エネルギー制御(20)…目標位置等速で移動させるのではなく，抵抗力の大きさによって移動量を変動させることで，人間に近い柔らかな動作が行え，良い結果を生むと仮説を立案した．ある一定以上の力を検出した場合，目標位置の 1 サイクルの移動量を以下に示すエネルギーを導入した式により決定する．

dx n =E ⋅ dTFn（3.2.1）

ここに， dx n ：1 サイクルの目標位置移動量， E ：設定エネルギー値， dT ：サンプリン
グタイム， Fn ：外力である．これにより，力が大きいほど，目標位置の移動量は低下す
る．シミュレーションの結果を Fig.3.2.5 に示す．なお，エネルギーの設定値は-0.1[J]とした．
・目標位置切り替え制御…通常の目標位置に加えて，通常は目標位置と同じ動きをするが，一定以上の力が加わった場合，移動量が仮想コンプライアンス制御される仮想目標値を導入した．この仮想目標値の計算に使う仮想パラメータは，現在位置の計算に使われるものと同一である．速度が急増した場合，目標位置をこの仮想目標値に切り替え，位置偏差を小さくすることで，速度超過の時間を短時間で終わらせる．仮想目標位置 Xv と目
標位置 X 0 の関係を以下に示す．

X 0 = X 0 + D Xv = Xv + D
X0 = X0+ D Xv = Xv + VX n ⋅ dT

（ Fn > Fa & Vn < Vlimit ）(3.2.2)

 X 0 = Xv
（ Vn > Vlimit ） Xv = Xv + D
ここに， D ：1 サイクルの位置偏差量， Fa ：仮想コンプライアンス制御開始抵抗力，Vlimit ：目標位置切り替え速度である．また， VX n はあるサイクルの仮想目標位置の移動速度で式(3.1.4)より求める．この時 dx n には Xv と X 0 の差分を使う．

10

運動療法における他動運動は，目的の位置まで確実に移動することが重要であり，目標位置を切り替える（後退させる）ことで動作時間が長くなったとしても，極端に延びなければ大きな問題では無いと仮定し，この制御法を考案した．シミュレーション結果
を Fig.3.2.6 に示す．なお，仮想目標値が仮想コンプライアンス制御を開始する力を-10[N]，切り替え速度を 0.07[m/sec] とした．

3.2.2 シミュレーション結果とその考察
Fig3.2.2 から Fig.3.2.6 までのシミュレーション結果について制御法毎に考察し，4 章の実験において扱う制御法を決定する．
・通常の仮想コンプライアンス制御…抵抗力急減時に予想通り速度の急上昇が見られ，その最大値は限界速度の 10 倍以上を示した．抵抗力に対する応答性や目標値への追従性は高く，動作終了までに要した時間は 28.7[sec]と最短ではあるが，安全であるとは言い切れない．
・速度リミッタ追加…限界速度を超過しなかった点は評価できる．また，動作終了までに要した時間は通常の仮想コンプライアンス制御と同じ 28.7[sec]である．しかし，リミッタが働いている時間が約 3[sec]あり，この間は速度が限界値を振り切っていることを示している．
現実世界では様々なパターンの抵抗力が加わる事が考えられるが，速度リミッタが働いている間に抵抗力が加わった場合，実際には速度超過を起こしている状態なので，抵抗力に対する応答がリミッタ非動作時と異なるという問題が起こる．これは好ましい事
ではなく，リミッタの動作時間を短くする必要がある．
目標位置に追い付くには，当然その移動速度より速く動かなければならず，速度が増加するのは仕方ない．この事が痙縮に対する仮想コンプライアンス制御の根本的な問題であろう．
・パラメータ変更…あらかじめ限界速度を超過しないように計算して求めたパラメータだ
けに，最大速度は 0.0344[m/sec]と限界値を越えていない．動作終了までに要した時間は31.6[sec]で，パラメータ変更前と比べて若干長くなっている．抵抗力に対する応答は鈍く，また目標速度に到達する前に抵抗力が発生する位置まで進んでしまうなど，抵抗力に対する速度の応答性が低いことが問題である．
・目標位置エネルギー制御…動作終了までに要した時間は，66.5[sec]と最も長く，通常の仮想コンプライアンス制御の場合と比べて 2 倍以上の時間を要した．また，抵抗急減時の速 度上昇も 通常の仮 想コンプライ アンス制 御より 若干低 いと はい え，最大 値 が11

0.3609[m/sec]と限界速度の 10 倍以上になっている．抵抗力が大きい場合に移動量を低下させることは，安全性や患者の受ける心理的な側面を考えると良い考え方ではある．しかし，これ単体では目的を達成できないと推察された．
・目 標 位 置 切 り 替 え 制 御 … 動 作 終 了 ま で に 要 し た 時 間 は 30.7[sec] ， 最 大 速 度 は
0.021[m/sec]と，シミュレーションを行った各制御法の中で最もバランスの良い結果となった．目標位置切り替え時に速度が 0.013[m/sec]低下しているが，マイナスにはなっておらず，これが人を使った実験だったならば被験者はそれほどの抵抗を感じないであろう．
懸念された通常の仮想コンプライアンス制御との時間差も約 2[sec]で，大きな差は出なかった．目標位置を切り替えて位置偏差を縮めることで，目標位置に追い付こうと増速する仮想コンプライアンス制御の問題点を上手く回避できた．
以上，各制御法のシミュレーション結果を説明した．各制御を単体で用いた場合，速度超過が起こる，動作時間が長時間になる等の問題があり，複数の方法を組み合わせる事で目的達成の可能性が高まると考えられる．特に速度リミッタは安全のために全ての制御法に付加するべきであろう．この結果を受けて 4 章の実験では以下に示す制御法を用いる．
1) 他 の 制 御 方 法 と の 比 較 の た め に ， 通 常 の 仮 想 コ ン プ ラ イ ア ン ス 制 御 （ VirtualCompliance Control）
2) 通常の仮想コンプライアンス制御＋速度リミッタ（Virtual Compliance control withvelocity Limiter 以下，VCL と略す）
3) 通常の仮想コンプライアンス制御＋パラメータ変更＋速度リミッタ（Parameterchanges Virtual Compliance control with velocity Limiter 以下，PVCL と略す）
4) 目標位置切り替え制御＋速度リミッタ（Switching target position control with velocity
Limiter 以下，SL と略す）
5) 目標位置切り替え制御＋目標位置エネルギー制御＋速度リミッタ（Switching target
position and Energy control with velocity Limiter 以下，SEL と略す）

3.3 関節角度推定プログラム
2 章で説明した筋張力シミュレータは，人間の関節角度をパラメータとして必要とする．
人間の関節角度を計測するにはゴニオメータなどの機器を使うのが一般的ではあるが，装着に掛かる時間や，装着後の煩わしさを軽減するために，ロボットアームの先端位置から逆運動学によって，人間の肩と肘の関節角度を推定するプログラムを作成した．
腕の動きが水平 2 自由度なので，ロボットアームと腕の位置関係は Fig.3.3.1 に示すようになる．まず，腕先位置 x H ， y H を求める．
x H = dX − (x A + d x )(3.3.1)y H = dY − y A + d y(3.3.2)()
12ここで( dX ， dY )はロボットアーム原点と人間の肩関節の距離，( x A ， y A )はロボットアー
ム先端位置，( d x ， d y )はロボットアーム先端位置と腕先位置との距離（オフセット量）である．
求めた腕先位置と関節角度 θ s ， θ e の関係は，
x H = l1 cos θ s + l 2 cos(θ s + θ e )(3.3.3)y H = l1 sin θ s + l 2 sin (θ s + θ e )(3.3.4)と示される． l1 と l 2 は上腕の長さと前腕の長さである．式(3.3.3)，式(3.3.4)から，

cosθ e =

(x

2
H

) (

+ y H − l12 + l 22
2l1 ⋅ l 2
2)

(3.3.5)−1

式(3.3.5)から cos を用いて θ e を求めることができるが， θ e の可動域は‐5 から 145[deg]−1であり，cos や sin−1−1

は 0，90 あるいは 180[deg]付近で誤差が大きくなり，また cos は 0

から 180[deg]まで， sin−1

は‐90 から 90[deg]までしか求められない．

そこで， sin θ 
(3.3.6)
θ = tan −1  cosθ 
を利用し， cos と sin を符号も含めて独立に与え，±180[deg]の範囲で θ e を求める．

()

いま， sin 2 θ e = 1 − cos 2 θ e より，

sin θ e = ±(2l1 ⋅ l 2 )2 − {x H2 + y H2 − (l12 + l 22 )}2(3.3.7)2l1 ⋅ l 2式(3.3.5)，(3.3.6)，(3.3.7)より，肘関節角度 θ e は，±
θ e = tan −1(2l1 ⋅ l 2 )2 − {x H2 + y H2 − (l12 + l 22 )}2 (x + y ) − (l + l )
2
H

2
H

2
1

2
2





(3.3.8)となる．なお， θ e は計算上負も存在するが，角度が負ということは肘が反り返った状態であり．人間の場合危険なので，本研究では．正の値のみを用いる．
求めた θ e より，肩関節角度 θ s は，
 (l + l cos θ e ) ⋅ y H − l 2 ⋅ x H sin θ e θ s = tan −1  1 2
 (l1 + l 2 cos θ e ) ⋅ x H + l 2 ⋅ y H sin θ e (3.3.9)となる．
式(3.3.8)と式(3.3.9)に，3.1 節や 3.2 節で説明した制御法で動作するロボットアームの先端位置を逐次入力として取り込み，計算することで人間の関節角度をリアルタイムに推定
できる．

13

4

実験
作成した筋張力シミュレータと，痙縮にも対応できるように考案したロボットアームの制御方法の有効性を確認するために実験を行った．

4.1 筋張力シミュレータ実験
実験の概要は，ロボットアーム先端に取りつけた力センサで，右腕の屈曲時と伸展時の腕先力を計測し，シミュレータを用いて推定筋張力を求める．同時に，個々の筋の筋張力を計測する事は困難であるため，表面筋電位計（Appendix 用語説明参照）を用いて，体表面の筋の表面筋電位（以下，IEMG）を計測し，両者を比較することで，シミュレータの有効性を確認する．

4.1.1 実験装置
実験に使用した機材は，力検出装置として 6 軸 Force/Torque センサ IFS-67M25A 50-I40（NITTA 製），筋電位計測装置としてサイナアクト MT11（NEC メディカルシステムズ製），
力センサの位置調整と腕先へ負荷を与えるために 7 軸ロボットアーム PA-10A-ARM（三菱重工業製）を用いた．各装置の仕様を Table4.1.1 から Table4.1.3 に示す．
この実験の概要図を Fig.4.1.1 に示す．ロボットアーム先端に力センサと，被験者が力を加え易いようにグリップが自由に回転するハンドルを取りつけた．また，被験者の肩・肘関節角度の検出は 3.3 節で説明した関節角度推定プログラムを利用する．そのためにアーム
原点から，被験者の肩関節までの距離を計測した．

4.1.2 実験方法
被験者は，年齢 23 歳の健常男性とした．被験者の上腕及び前腕の長さと，筋長補正係数を Table4.1.4 に示す．シミュレータで用いる筋の生理断面積（ PCSA ）の比例定数 k は，被験者が健常男性であるため， k ＝5.5 と設定した．
計測する筋は，シミュレータで扱っており，かつ体表面に存在し表面筋電位計で計測できる上腕二頭筋：Biceps Brachii（以下 BB，屈曲時使用）と，TB（伸展時使用）とした．
計測までの準備として，アルコールで皮脂を除去した上で，Delagi(21)らの測定位置で筋肉の Motor Point を特定し，被験者の BB と TB にディスポーザブル電極を貼り付けた．
Fig.4.1.2 にその貼り付け位置を示す．次に，被験者にイスに着座してもらい，右肩関節を90[deg]外転，右肘関節を 90[deg]屈曲した姿勢をとらせ，その姿勢でハンドルを把持できるようにロボットアームを移動させる．Fig.4.1.3 に実験の様子を示す．
計測は，ロボットアームを停止させて，腕先位置を変動させない状態で，水平面での腕の屈曲と伸展をそれぞれ行う場合と，ロボットアームをスティフネス・ダンピング制御し，腕先位置を変動させながら，水平面での腕の屈曲と伸展をそれぞれ行う場合．以上の 4 パ
ターンについてそれぞれ 3 回ずつ行った．1 回の計測は 20[sec]とし，5[sec]間無理のない程

14

度の屈曲・伸展を行い，次の 5[sec]間は力を抜く．これを繰り返し，ステップ状の力を発生させた．筋電位計と力センサのサンプリングタイムは 0.001[sec]．ロボットアームの制御周期は，アーム性能の限界のため 0.01[sec]とした．よって，スティフネス・ダンピング制御
時には力センサの情報は 0.01[sec]毎にロボットアームに反映される．
なお，実験で用いたスティフネス・ダンピング制御の制御式は以下に示すものである，

C
dt


(4.1.1)⋅  Fn + K ⋅ x0 + ⋅ x n K ⋅ dt + C dt
ここに， x n +1 ：目標位置， x n ：現在位置， x0 ：基準位置， dt ：サンプリングタイム， Fn ：x n +1 =力センサで検出した外力， C ：仮想粘性定数， K ：仮想バネ定数である．実験ではパラメ
ータ C ＝90， K ＝1000 と設定した．これは 3.1.2 項で説明したものと同じものである．

4.1.3 実験結果
実験により得られた腕先位置から推定関節角度を求め，力センサで計測した腕先力と共に，シミュレータに入力して各筋の推定筋張力を得た．また，計測した IEMG はローパスフィルタに通し，全波整流（絶対値化）した後に積分して，積分筋電波形に加工した．積分筋電波形では無負荷時でも 15〜25[μV]前後の電圧が計測された．これは姿勢や腕先位置の維持のために発生している IEMG か，あるいは電極の貼り付け状態等によるノイズであり，腕先力の発生には関与し無いと考えられる．各筋の無負荷時の平均電圧は BB＝16.2[μV]，TB＝25.9[μV]であり，これらを腕先力が発生する有効な筋活動の有無の目安として，
積分筋電波形内に表示する．
筋電整流波形，積分筋電波形，腕先力，推定筋張力，腕先位置と角度の変動値，以上のデータを実験パターン毎にまとめて，Fig.4.1.4 から Fig.4.1.7 に示す．なお，図中の TB は，
上腕三頭筋長頭：TBLo と上腕三頭筋外側頭（内側頭を含む）：TBLa の合計値である．

4.2 ロボットアーム制御実験
実験の概要は，痙縮のモデル(22)を与えた 2 軸の Direct Drive（DD）ロボットを痙縮右上肢に見立て，これを先端に力センサを取りつけたロボットアームで他動的に伸展させる．
そして，抵抗が急減した時のロボットアームの挙動を 3.2 節で述べた制御方法毎に調べ，関節角速度が許容値以内であるか等を確認し，考案した制御方法の痙縮に対する有効性を確かめる．また，DD ロボットの関節角度：Jrad は 3.3 節で説明した関節角度推定プログラ
ムを用いて求めるが，それとは別にプログラムの正確さを確認するために，DD ロボットのエンコーダ出力も記録する．

4.2.1 実験装置
実験に使用した機材は，療法士に見立てた 7 軸ロボットアーム PA-10A-ARM（三菱重工
業製），力検出装置として 6 軸 Force/Torque センサ IFS-67M25A 50-I40（NITTA 製），痙

15

縮右上肢に見立てた 2 軸 DD ロボット SR-402DD(S)（東京エレクトロニックシステムズ製）
以上である．DD ロボットの仕様を Table4.2.1 に，実験の光景を Fig.4.2.1 に，概要図を
Fig.4.2.2 に示す．
ロボットアーム先端に力センサを取り付け，DD ロボットはクランプで台座に固定した．
また，2 台のロボットの接触点には衝撃緩和のためにスポンジゴムを貼り付けた．

4.2.2 実験方法
まず，DD ロボットについて説明する．DD ロボットは第 2 軸を固定して，リンク長：

link =0.4[m]の 1 軸ロボットとして用い，初期角度を 1.5708[rad]（90[deg]）とした．また，DD ロボットに導入した痙縮モデルは本質的にインピーダンス制御法の制御式と同じ物であり，その制御式は式（3.1.1）になる．
仮想パラメータは本来ならば実際に痙性麻痺の患者から計測したデータを使うことが望ましい．本研究では研究の初期段階であること．また，患者の同意を得ることが困難である た め 以 下 の パ ラ メ ー タ を 仮 定 し た ． 実 験 開 始 時 か ら の DD ロ ボ ッ ト の 角 偏 差 が0.5236[rad]（30[deg]）となった時点で，DD ロボットの質量や摩擦も含めた抵抗力が 40[N]となるようにパラメータを，第 1 関節が M =0，C =15，K ＝32，第 2 関節が M =0，C =30，

K ＝64 として関節トルクを制御する．また，抵抗がある時点で急減するという痙縮の特徴を再現するために，DD ロボットの関節角度が 1.0472[rad]（60[deg]）以下となった時点で関節トルク 0 となるように設定した．
以上のように動作する DD ロボットの関節角度を 1.5708[rad]（90[deg]）から設定目標角度：rad=0.7854[rad]（45[deg]）まで動かすようにロボットアームに軌道を与え，3 章で説明した制御方法で動作させデータを計測する．この時，DD ロボットのリンク中心線からロボットアーム先端位置までの間に，力センサ，
スポンジゴム，DD ロボットの幅等で 0.095[m]の距離が存在する．これをオフセット量と考えて，ロボットアームの軌道を決定した．
仮想コンプライアンス制御のプログラムフローチャートを Fig.4.2.3 に示す．VCL 及びPVCL のプログラムフローチャートを Fig.4.2.4 に，SL のプログラムフローチャートを
Fig.4.2.5 に，SEL のプログラムフローチャートを Fig.4.2.6 に示す．
痙縮に対する他動運動は関節角速度 3〜5[deg/sec]で行われるので，各制御法とも抵抗力が発生しない場合は，0.021[m/sec]（3[deg/sec]）の速度で DD ロボットを動かし，許容限界速度は 0.035[m/sec]（5[deg/sec]）とした．
サンプリングタイムは DD ロボットが 0.001[sec]，ロボットアームが 0.01[sec]である．

4.2.3 実験結果
実験により得られた推定角度と DD ロボットのエンコーダから出力された角度，抵抗力，
ｘ・ｙ方向の目標位置と現在位置，ｘ・ｙ方向の現在位置と目標位置の差，ｘ・ｙ方向の速度．以上のデータを制御法毎に Fig.4.2.7 から Fig.4.2.11 に示す．

16

5 考察
5.1 筋張力シミュレータに関する考察
本節では Fig.4.1.4 から Fig.4.1.7 に示す，筋張力シミュレータ実験の結果について考察
する．
Fig.4.1.4 は腕先位置を固定した状態で BB の作用によって肘関節を屈曲させようとした運動である．BB の推定筋張力図と積分筋電波形の全般的な傾向を見ると，2 度のステップ
状の反応がほぼ同時刻に見られ，これらは腕先力の発生とも一致し，両者に関連が有ることを示している．
TB に関しては，推定筋張力図の値は 0 で，活動していない事を示しているが，積分筋電波形では BB に比べれば低いもののステップ状の反応が見られる．この反応のピークは BBのそれと比較すると少し遅れて現れている．これは腕先位置を固定した状態で実験を行ったので，関節を屈曲させる力に対する，力センサを取りつけたロボットアームからの反力が大きく影響していると推察される．また，力を加えていない時に計測されている 15〜
20[μV]前後の積分筋電波形は腕先力の発生には関与していないので，これ以下の電圧は無視して良いと考えられる．
Fig.4.1.5 は腕先位置を固定した状態で TB の作用によって肘関節を伸展させようとした運動である．Fig.4.1.4 と同じく，TB の推定筋張力図と積分筋電波形には一定の傾向が有ることが推察できる．また，腕先力が Fig.4.1.4 に比べて減少しているにも関わらず，TB の
積分筋電波形の最大値が Fig4.1.4 の BB のそれより大きい．皮膚抵抗等の条件が異なるので数値的な比較はできないが，肘関節の伸展と屈曲では使用する筋の数に違いがあり，伸展時の方が少ない．したがって，同じ腕先力を発生させるには，個々の伸展筋は屈曲筋よりも大きな張力を要求されると考えられ，この実験結果は傾向としては正しいと推察される．TB の推定筋張力図をみるとこの現象を良く表しており，Fig4.1.4 の BB のそれより大きな数値を算出している．
Fig4.1.6 と Fig.4.1.7 は腕先位置を変動させて，腕の屈曲と伸展を行ったもので，運動に伴い肩関節も若干変化している．定性的な比較のために，屈曲時の BB と伸展時の TB の，
積分筋電波形と推定筋張力，それぞれの最大値を 1 として無次元化したものを Fig.5.1.1 とFig.5.1.2 に示す．推定筋張力図と積分筋電波形の関係は，腕先位置を固定した場合と比較
してピーク等のズレが小さく，より高い類似性を示した．また，その運動では使わない筋の積分筋電波形は，腕先位置固定時のようなステップ状の反応を示さなかった．これは加えた力が腕先位置の変化に使われ，過度の反力が発生しなかったためと考えられる．
以上の結果より，本研究で作成した筋張力シミュレータは，全般的傾向において，積分
筋電波形との間に類似性が見られ，腕先位置の固定の有無によらず，筋の働きを良く表していると言える．もともと，IEMG は数値一つ一つを見るので無く，筋活動の有無やその強さ，活動の開始時刻や持続時間，あるいは動作に伴う複数の筋相互間の協調性や相反性を知るためにある一定の時間内での全体的な変動を見るものである．したがって，作成し17

たシミュレータは，表面筋電計のように準備・装着に時間を取られることなく，回復課程の評価に用いる指標を提供できる可能性が高いと言える．また運動の分析や作業における筋負担の人間工学的評価の推定に有用であり，患者に適した軌道の決定を行えるものと推察される．
このシミュレータの今後の課題としては，腕先位置を固定した場合に本来動いていない
はずの IEMG が発生し，推定筋張力との間に相違が見られた点を改善すべきである．これは反力の影響だと考えられるが，作成したシミュレータは力センサの出力である腕先力のみを用いて計算を行うので，この現象に反応しなかった．また，今回の実験では推定筋張力：[N]と，積分筋電波形：[μV]との相関関係を統一的に分析しなかったが，今後多くの被験者で実験を行い，多変量解析を行って筋張力の推定を行えば，シミュレータの有意性が一層増大するであろう．

5.2 ロボットアームの制御方法に関する考察
本節では Fig.4.2.7 から Fig.4.2.11 に示す，ロボットアーム制御実験の結果について考察する．
まず関節角度については，3.4 節で説明した関節角度推定プログラムの出力と，DD ロボットのエンコーダ出力の比較を行った．結果は最も差の大きかった通常の仮想コンプライ
アンス制御において，瞬間最大誤差 0.1[rad]（5.73[deg]），平均誤差 0.035[rad](2.0[deg])となった．逆に誤差の最も小さかった SL では平均誤差 0.008[rad]（0.46[deg]）である．
平均誤差を距離に換算すると 0.003〜0.014[m]であり，エンコーダ出力の変動に対してもよく追従している．したがって，このようなケースにおいて関節角度推定プログラムはゴニ
オメータ等の代役を十分果たしていると考えられる．
以下，各制御法の結果について考察する．Fig.4.2.7 は他の制御方法との比較のために行った通常の仮想コンプライアンス制御の結果である．シミュレーションの段階で判明していたが，実験においても抵抗急減時に速度が急上昇しており，ｘ・ｙ両方向の速度を合成した数値は最大で 0.18[m/sec]となった．これは許容限界速度の 5 倍以上であり，速度超過の状態は 0.5[sec]間続いた．この間アームは進行方向に約 0.04[m]進んでおり，実験におけるトータルの移動量が約 0.32[m]であることを考えると，全体の 10%以上の行程を速度超過の状態で進んでいる．このような状態で適切な運動療法が実施できるとは考え難く，通常の仮想コンプライアンス制御は痙縮患者の他動運動のようなケースには向いていないことが実証できた．なお，動作終了までに要した時間は 15.7[sec]である．
Fig.4.2.8 は VCL の実験結果である．動作終了までに要した時間は 15.8[sec]，このうち
速度リミッタが動作していた時間は約 3.0[sec]とシミュレーションと同様の結果となった．
よって，通常の仮想コンプライアンス制御と速度リミッタの組み合わせでは，リミッタの動作時間が長くなり，シミュレーション結果で述べた問題が現実に起こることがわかった．
Fig.4.2.9 は PVCL の実験結果である．動作終了までに要した時間は 19.9[sec]と VCL に

18

比べて長い．仮想パラメータ M ， C を大きくした結果速度変動が小さく，抵抗力急減時の速度上昇も許容限界を越えていない．しかし，あまりにパラメータが大きいために抵抗力に対する応答が極めて低く，抵抗力が加わっても速度が上昇を続ける，抵抗力のピークと速度のピークが 2.1[sec]ずれるなどの現象が現れた．限界速度を越えないという目的は達成できたが，力に対する応答性がこれほど低いと，力制御の意味が無いと言っても過言ではない．しかし，パラメータの数値を下げると Fig.4.2.7 に示すように限界速度を越えてしまう．仮想コンプライアンス制御のパラメータを調整するだけでは想定した状況に対応できないことが確認された．VCL，PVCL 共に許容限界速度を超過していないが，力に対する応答に問題がある．これは人間が使用した場合に違和感を生じさせるものと考えられ，心理的側面に悪影響を及ぼしかねない．
Fig.4.2.10 は SL の実験結果である．動作終了までに要した時間は 17.7[sec]で，通常の仮想コンプライアンス制御よりも若干長い．抵抗力急減による速度上昇時に，速度リミッタが動作した時間は 0.03[sec]と短時間であった．これは VCL における速度リミッタ動作時間の 100 分の 1 に過ぎない．この程度の時間であれば，速度リミッタ動作時の抵抗力に対する応答の変化も殆ど問題にならないものと推察される．目標位置切り替え時に，x 方向の速度が 1.1[sec]に渡って低下しているが，これはシミュレーションでも同様の結果が出たように，切り替えによって位置偏差が急激に低下したためである．しかし，他動運動において
は目標位置まで確実に移動する事が重要であり，SL はこれを達成している．よって，一時的な速度低下により動作時間が若干長くなっても，深刻な問題にはならないと考える．以上より，SL は VCL，PVCL に比べて抵抗力に対する応答が良く，速度超過も起こさないという好結果が得られた．
Fig.4.2.11 は SEL の実験結果である．動作終了までに要した時間は 68.9[sec]と他の制御法に比べて約 4 倍の時間が掛かっている．これは設定エネルギー値 E を低く設定したために目標位置の進行が遅くなった事が原因であり， E を調整することで解消されると考えられる．基本的に SEL は SL と同様の動作特性を持つので，動作終了までの時間が長いことを除けば，他に問題となる点は無い．
以上の結果を要約すると，痙縮にも対応したロボットアームによる上肢運動療法に要求される

動作中常に抵抗力に対して敏感に反応しながら目標位置まで対象を移動させ，抵抗急減時に限界速度をこえないという動作を実現させるには，通常の仮想コンプライア

ンス制御に速度リミッタの追加，パラメータの変更を施しても困難であり，目標位置の切り替えを導入することで要求を満足できた．
共に要求を満たす SL と SEL の優劣については，動作時間の短さから言えば SL の方が優れていると言えるが，抵抗力が大きければ移動量を低下させるという SEL の特性は，安全性や患者の受ける心理的な面を考えると利点が有るように思える．今後，エネルギー設
定値を調整して人間を用いた実験を行い，感性面での評価データを揃えた上でどちらの制
御方法が優れているかを検討する必要があろう．

19

6 結言
本研究では，痙性麻痺にも対応したロボットアームによる上肢運動療法の実現を目標に，
筋張力シミュレータと，仮想コンプライアンス制御を改良した制御法を考案し，基礎的実験を行うことでその有効性を確認した．その結果，以下のような結論を得た．
(1)

表面筋電位計を用いて，人の BB 及び，TB の IEMG を計測し，同条件のシミュレー
タの結果と比較した．この結果，両者には関連性が見られ，作成したシミュレータの出力である推定筋張力は運動療法の効果等の評価指標として使用できる可能性が高いと考えられる．また，ロボットアームを用いた運動療法を実施する前に，シミュレー
タに軌道と腕先力を入力することで，各筋の予測活動状態を推定することが可能であり，アームの軌道決定など機能回復訓練の効率を上昇させる効果も望める．

(2)

擬似的な痙縮の特性を持たせた DD ロボットを患者の腕に見立て，ロボットアームを用いた上肢の伸展運動を様々な制御方法で行った．この結果，通常の仮想コンプライ
アンス制御では抵抗力急減時に，速度超過となる等の問題が発生し，速度リミッタの追加やパラメータの変更を行っても満足できる結果が得られないことも確認された．

(3)(2)を受けて，問題解消のために目標位置の切り替えを導入した改良型仮想コンプライ
アンス制御法 SL，SEL を考案した．これらを用いた実験を行い，抵抗急変時においても，速度が許容範囲内に収まり，抵抗力に対する応答も良好であることを確認した．
したがって，両制御法とも，痙縮に対応したリハビリを実施できる可能性が高い．

以上より，考案したシミュレータと制御法は，上肢リハビリにおいて有効であると考えられる．今回行った実験は健常者やロボットを用いて行ったものであり，また制御法とシミュレータはそれぞれ独立した状態でその有効性を確認している．今後の課題としては，シミュレータの結果からアームの軌道を生成し，考案した制御法で動作させるシステムを構築した上で，人間を用いた実験や，実際にリハビリを必要とする患者のデータ解析を行い，より実用的なレベルの検討を行う必要がある．

謝辞
本研究を進めるにあたり，終始丁寧なご指示，ご指導を賜りました指導教員の井上喜雄
教授，甲斐義弘助手，副指導教員の八田章光助教授に，深く感謝の意を表し，お礼の言葉とさせていただきます．また，医学的見地から貴重な助言をして下さった川崎医療福祉大学の谷岡哲也先生，菅原憲一先生に感謝の意を表します．
実験に参加して頂いた，大学院生の川澤延弘氏，松田拓也氏，森井康友氏，学部生の戸次亨弐氏をはじめ，研究のサポートをして頂いた知能機械力学研究室の研究生の皆さん，
ご協力ありがとうございました。

20

参考文献
(1)

厚生省大臣官房統計情報部： 平成１１年 医療施設（静態・動態）調査・病院報告の概況 ，（2000）(2)

永田他： 多自由度・力制御可能なパラレル型ＣＰＭ装置 ，第 10 回リハ工学カンファレンス抄録集，（1995），pp.323-326(3)岡島他： 間接可動域訓練装置：柔らかさを与える機構の検討 ，総合リハ，vol.26，（1998），pp.363-369．
no.4，

(4)

則次，安藤，山中： ゴム人工筋を用いたリハビリテーション支援ロボット（第 1 報

インピーダンス制御による訓練運動モードの実現），日本ロボット学会誌，vol.13，no.1，（1995），pp.141-148．
(5)

坂口，古荘，元田： ER アクチュエータを用いたリハビリテーション訓練システムの開発に関する基礎研究 ，日本ロボット学会誌，vol.19，no.5，（2001），pp.612-619．

(6)

平林，杉本他： 多自由度ロボットの仮想コンプライアンス制御 ，計測自動制御学会
論文集，vol.22，no.3，（1986），pp.343-350．

(7)

林原，園田他： 人とロボットによる長尺物の協調運搬（人の協調挙動に基づく鉛直平
面内の制御方法の検討） ，日本機械学会論文集（C 編），vol.67，no.653，(2001)，
pp.162-169．

(8)

五味，川人： 人腕運動中の手先インピーダンスの計測 ，ロボティクス・メカトロニ
クス講演会 講演論文集，no.940‐21，（1994），pp.835-840．

(9)

岩本： 筋肉の特性を考慮した操作力多角形による上肢の力学的評価 ，同志社大学大学院修士論文，（1997）(10) P.Piegon，L.H.Yahia and A.G.Feldman： Moment Arms and Lengths of Human
Upper Limb

Muscles as Functions of Joint Angles ，J.Biomechanics，vol.29，

no.10，(1996)，pp.1365-1370(11) Rutherford O.M and Jones D.A： Contractile Properties and Fatiguability of the
Human Adductor Pollicis and the Effects of Two Chronic Stimulation Patterns,

J.Neuirology and Science，vol.85，(1988)，pp.319-331(12) 日本機械学会： 生体工学 ，オーム社，（1991）(13) K.R.Kaufman，K.N.An and E.Y.S.Chao： Incorporation of Muscle Architecture into
the Muscle Length-Tension Relationship ，J.Biomechanics，vol.22，no.8/9，(1989)，
pp.943-948(14) 茨木，福島： Fortran77 最適化プログラミング ，岩波書店，(1991)(15) K.N.An，K.R.Kaufman and E.Y.S.Chao： Physiological Considerations of Muscle
Force Through the Elbow Joint ，J.Biomechanics，vol.22，no.11/12，(1989)，
pp.1249-1256(16) J.E.Wood ， S.G.Meek and S.C.Jacobsen ： Quantitation of Human Shoulder

21

Anatomy for Prosthetic Arm Control – 1.Surface Modeling ，J.Biomechanics，vol.22，
no.3，(1989)，pp.273-292(17) K.N.An，F.C.Hui，B.F.Morrey，R.L.Linscheid and E.Y.Chao： Muscle Across the
Elbow Joint a Biomechanical Analysis ，J.Biomechanics，vol.14，no.10，(1981)，
pp.659-669(18) H.E.J.Veeger，Bing Yu，K.N.An and R.H.Rozendal： Paramaters for Modeling theUpper Extremity ，J.Biomechanics，vol.30，no.6，(1997)，pp.647-652(19) H.E.J.Veeger ， F.C.T.Van Der Helm ， L.H.V.Van Der Woude ， G.M.Pronk and
R.H.Rozendal： Inertia and Muscle Contraction Parameters for Muscleskeletal
Modeling of the Shoulder Mechanism ，J.Biomechanics，vol.24，no.10，(1991)，
pp.615-629(20) 山崎，横川，原： 柔らかな物体に接触する関節形マニピュレータに対する位置と力の制御モードの一切換え法 ，同志社大学理工学研究報告，vol.39，no.1，（1998），
pp.21-32．
(21) Edward F.Delagi，Aldo Perotto： 筋電図のための解剖ガイド‐四肢‐ ，西村書店，(1985)（1973）(22) 加藤他： 続人間の手足の制御 ，学献社，

22Fig.2.1.1 Coordination and Term Definition

Fig.2.1.2 Muscle Model of Human Arm

23

Table2.1.1 List of Human Arm Muscle
Name of Muscle

Abbreviated Title

Japanese

Brachialis
Brachioradialis
Anconeus
Triceps Brachii Caput Laterale
Triceps Brachii Caput Longum
Biceps Brachii
Anterior Portion of the deltoid
Posterior Portion of the deltoid
Clavicular Portion of Pectoralis

BS
BR
AN
TBLa
TBLo
BB
DA
DP
PC

上腕筋
腕橈骨筋
肘筋
上腕三頭筋外側頭
上腕三頭筋長頭
上腕二頭筋
三角筋前部
三角筋後部
大胸筋鎖骨部

Nomalized Muscle Tension

1.4

Active Muscle Tension
Passive Muscle Tension
Maximum Isometric Tension

1.2
1
0.8
0.6
0.4
0.2
0
-0.5

0

0.5

1

Muscle Strain

Fig.2.1.3 Relationship between Muscle Strain and Normalized Muscle Tension

24

StartSelecting Initialize active set J

Calculate

J ‐Optimal Solutions ∉ J exists thatNO

meets the conditionc s (x ) ＜0

YES

End：
OptimalSolutionz = H J as
r = AJ+ a s

z=0

YES

NO

t1 =


 v (x )v k (x )ri > 0 , i ∈ J ∩ K = min  k
rk

 rk

t2 = −i ∈ J ∩ K existsthat meets thec s (x )NO

a sT zcondition ri ＞0

End：
YES

Calculatet1 ≥ t 2NO

Impossible

YES

x := x + t 2 z

x := x + t1 z

x not change

J := J − {k }J := J − {k }(full step)(partial step)(partial step)J := J ∪ {s}Fig.2.2.1 Flowchart of Dual Method Algorithm
25MuscleWood(14)

BS
BR
AN
TBLa
TBLo
BB
DA
DP
PC

11.01
2.56

Table2.2.1 Data of PCSAPCSA[ｃｍ２］
An(15)
Veeger(16)Veeger(17)
6.97
1.47
2.48
12.13
6.70
4.62

9.66
7.10
10.30
6.46
4.54
8.38

5.34
2.78
1.08
7.42
3.71
5.16

Average
7.77
2.27
1.78
9.74
6.09
6.59
6.93
6.72
6.76

6.84
6.29
7.39
8.89
5.13

Table2.2.2 Standard Phalanx Length of Arm
Upper Arm
Forearmln1
ln2
Length
310
287
[mm]

Table2.2.3 Parameters of Arm Muscle
Name of Muscle

Abbreviated title

mln [mm]ia

Brachialis
Brachioradialis
Anconeus
Triceps Brachii Caput Laterale
Triceps Brachii Caput Longum
Biceps Brachii
Anterior portion of Deltoid
Posterior portion of Deltoid
Clavicular portion of Pectoralis

BS
BR
AN
TBLa
TBLo
BB
DA
DP
PC

98.0
170.8
50.0
152.4
207.2
169.3
137.9
174.1
137.7

0.95
0.72
0.67
0.51
0.51
0.67
0.67
0.67
0.67

26

Startθ sn , θ en , PCSA , ln , mln , ia , kl
mli = mlni ⋅ ilni

Fx ， Fy ， θ s ， θ eCalculate

f max n (θ n ) ， T ， J

Dual Method：
Algorithm

Calculatefn

EndFig.2.2.2 Flowchart of Estimated Muscle Tension Program

27

Fig.3.1.1 Virtual Mechanical Impedance

10
0
0

0.1

0.2

0.3

0.4

Force[N]

-10
-20
-30
-40
-50
Position[m]Fig.3.2.1 Relationship between Force and Position

28

0.5

0.6

10
Force[N]

0
-10 0

5

10

15

20

25

30

20

25

30

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

Current Position
Target Position

0

5

10

15
Time[sec]
Position

Distance[m]

0
-0.01

0

5

10

15

20

25

30

25

30

-0.02
-0.03
-0.04
-0.05
Time[sec]

Velocity[m/sec]

Distance between Current Position and Target Position

0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
0

5

10

15
Time[sec]

20

Velocity

Fig.3.2.2 Simulation Results of Virtual Compliance Control
29

10Force[N]

0
-10 0

5

10

15

20

25

30

-20
-30
-40
-50

Position[m]Time[sec]
Force

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

Current Position
Target Position

0

5

10

15
Time[sec]
Position

20

25

30

10

15

20

25

30

25

30

0
Distance[m]

-0.01

0

5

-0.02
-0.03
-0.04
-0.05
Time[sec]

Velocity[m/sec]

Distance between Current Position and Target Position

0.04
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0
0

5

10

15
Time[sec]

20

Velocity

Fig.3.2.3 Simulation Results of Virtual Compliance Control with Velocity Limiter30

10

Force[N]

0
-10

0

5

10

15

20

25

30

35

20

25

30

35

20

25

30

35

30

35

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1 0Current Position
Target Position

5

10

15
Time[sec]

Distance[m]Position

0-0.02 0
-0.04
-0.06
-0.08
-0.1
-0.12
-0.14

5

10

15

Time[sec]

Velocity[m/sec]

Distance between Current Position and Target Position

0.04
0.035
0.03
0.025
0.02
0.015
0.01
0.005
0
0

5

10

15

20

25

Time[sec]
Velocity

Fig.3.2.4 Simulation Results of Virtual Compliance Control (Parameter Change)3110Force[N]

0
-10 0

10

20

30

40

50

60

70

40

50

60

70

40

50

60

70

60

70

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

Current Position
Target Position

0

10

20

30
Time[sec]
Position0
Distance[m]

-0.01

0

10

20

30

-0.02
-0.03
-0.04
-0.05
Time[sec]

Velocity[m/sec]

Distance between Current Position and Target Position

0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-0.05 0

10

20

30

40

50

Time[sec]
Velocity

Fig.3.2.5 Simulation Results of Energy Control
32

10Force[N]

0
-10 0

5

10

15

20

25

30

35

20

25

30

35

20

25

30

35

30

35

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

Current Position
Target Position
Virtual Target Position

0

5

10

15
Time[sec]
Position

0.01
Distance[m]0-0.01 0

5

10

15

-0.02
-0.03
-0.04
-0.05
Time[sec]
Distance between Current Position and Target Position

Velocity[m/sec]

0.025
0.02
0.015
0.01
0.005
0
0

5

10

15

20

25

Time[sec]
Velocity

Fig.3.2.6 Simulation Results of Switching Target Position Control
33

Fig.3.3.1 Coordination and Term Definitionin Estimated Joint Angle Program34

Table4.1.1 Specification of 6 axes F/T Sensor
6 軸力覚センサ IFS-67M25A 50-I40（NITTA 製）
測定限界ｘ，ｙ方向力
測定限界ｚ方向力
測定限界各軸モーメント
外径
高さ
重量
データ出力周期200[N]
400[N]
12.5[N･m]
67[mm]
25[mm]
180[g]
8000[Hz]

Table4.1.2 Specification of Electromyograph
多用途テレメータ サイナアクト MT11（NEC メディカルシステムズ製）
チャンネル数
14 チャンネル，マーク 1 チャンネル
電極入力数
最大 28
測定現象
脳波，心電図（心拍数），表面筋電図，眼球運動，呼吸曲線（1 チャンネル）伝送情報
バッテリエンプティ，ナースコール，モンタージュ番号，測定条件
電波形式
特定小電力医用テレメータ E 型
電波の到達距離約 30[m]以上（見通し可能な場所）
周囲温度…-10～60℃
保存温度
相対湿度…30～95%
周囲温度…10～40℃
使用環境
相対湿度…30～85%
入力インピーダンス
10MΩ以上最大入力電圧±3m[V]/±300μ[V]以上
入力重畳直流電圧±0.4[V]
同相成分除去比70[dB]以上（入力平衡）雑音
7μVP‐P（f=60[Hz]）周波数特性0.1～200[Hz]利得60/80[dB]
時定数
0.3/1.5[s]
出力電圧
0.5[V]/500μ[V]；50μ[V]（レコーダ用出力）
校正信号
500/50μ[V]（インターバル 4 秒の方形波）
送信機…138×82×37.5[mm]寸法
受信機…210×285×226[mm]送信機…300[g]（電池含まず）
重量
受信機…6.8[kg]35Table4.1.3 Specification of 7 axes Robot Arm
7 軸ロボットアーム PA-10A-ARM（三菱重工業製）
アーム本体重量
35[kg]
先端取扱重量10[kg]
関節数7
取りつけ面より R-P-R-P-R-P-R
関節構成（R は回転関節，P は旋回関節を示す）
ＡＣサーボモータ＋電磁ブレーキ＋駆動方式
アブソリュート角度検出器（レゾルバ）
手先合成最大速度
1550[mm/s]
位置繰り返し精度±0.1[mm]
肩リーチ 315[mm]+上腕 450[mm]
全長+下腕 500[mm]+手首リーチ
80[mm]=1345[mm]アーム長950[mm]（関節間距離）
アーム本体保護構造
防塵、防滴型（ＩＰ６５）
コントローラ質量25[kg]
各軸値入力，位置姿勢入力：６方向キ
教示方式
ーイン，ダイレクトティーチング
ＰＴＰ（点間）制御（直線補間，円弧補経路制御方式間，各軸），ＣＰ（連続軌道）制御
オープンアーキテクチャ（ＩＳＡバス内蔵
システム構成ＤＯＳ／Ｖ仕様）
ＡＣ100±10%，50/60[Hz]電源
1.5ｋVA 以下
0℃～50℃，
周囲温度
耐放射線 1Ｘ106[rad]（アーム本体）
ISA バス内蔵 DOS/V、ISA バス内蔵
コントローラバリエーション
Windows95/98/NT/英語版 NT36

PC (Muscle Tension Simulator)

Simulated Output:
Estimated Muscle Tension

Moving Direction

Force &Arm Position

Robot Arm
PA-10
6 axesForce/TorqueySensorRight Armx
Direction of the Position

Myogenic Potential of& Force Coordinate

Electromyograph

Biceps Brachii

Myogenic Potential of
Triceps Brachii
Data of Myogenic Potential
PC

Electromyogram of Muscles

Fig.4.1.1 System Structure of Simulator Experiment37

Table4.1.4 Phalanx Length of Subject Arm
Upper Arm

Forearm

Correction Factor of Muscle Length

ln1

ln2

ｌi/lni

280

230

0.852

Length[mm]

Triceps
BicepsBrachii

Brachii
Earth

Fig.4.1.2 Affix Point of ElectrodeFig.4.1.3 Experiment Scene

38Volt[μV]

150
100
50
0
0

2

4

6

8

10

12

14

16

18

20

14

16

18

20

Time[sec]
Biceps Brachii Myogenic Potential

Volt[μV]

100

50

0
0

2

4

6

8

10
Time[sec]

12

Triceps Brachii Myogenic Potential

Volt[μV]

150
Borderline of Effective
Muscle Activity

100
50
0
0

2

4

6

8

10
12
14
16
Time[sec]
Biceps Brachii Myogenic Potential (Post-Processing)

18

20

Volt[μV]

100

Borderline of Effective
Muscle Activity
50

0
0

2

4

6

8

10
12
14
16
Time[sec]
Triceps Brachii Myogenic Potential (Post-Processing)

18

20

Fig.4.1.4 (1) Experiment Results of Curvature Movement in Case of Fixed Arm Position
39

10

Force[N]

0
-10 0

2

4

6

8

10

12

14

16

18

20

10
12
14
Time[sec]
Estimated Biceps Brachii Tension

16

18

20

16

18

20

-20
-30
-40
-50
Time[sec]

Muscle Tension[N]Force of Arm Tip

150
100
50
0

Muscle Tension[N]

0

2

4

6

8

100

50

0
0

2

4

6

8

10
Time[sec]

12

14

Estimated Triceps Brachii Tension

Fig.4.1.4 (2) Experiment Results of Curvature Movement in Case of Fixed Arm Position40

Volt[μV]

100

50

0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

14

16

18

20

Biceps Brachii Myogenic Potential

300
Volt[μV]

250
200
150
100
50
0
0

2

4

6

8

10
Time[sec]

12

Triceps Brachii Myogenic Potential

Volt[μV]

100
Borderline of Effective
Muscle Activity
50

0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

Volt[μV]

Biceps Brachii Myogenic Potential (Post-Processing)

300
250
200
150
100
50
0

Borderline of Effective
Muscle Activity

0

2

4

6

8

10
Time[sec]

12

14

16

18

20

Triceps Brachii Myogenic Potential (Post-Processing)Fig.4.1.5 (1) Experiment Results of Extension Movement in Case of Fixed Arm Position41

Force[N]

30
20
10
0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

14

16

18

20

14

16

18

20

Muscle Tension[N]Force of Arm Tip

100

50

0
0

2

4

6

8

10
Time[sec]

12

Muscle Tension[N]Estimated Biceps Brachii Tension

350
300
250
200
150
100
50
0
0

2

4

6

8

10
Time[sec]

12

Estimated Triceps Brachii Tension

Fig.4.1.5 (2) Experiment Results of Extension Movement in Case of Fixed Arm Position42

Volt[μV]

100

50

0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

14

16

18

20

Biceps Brachii Myogenic Potential

Volt[μV]

100

50

0
0

2

4

6

8

10
Time[sec]

12

Triceps Brachii Myogenic Potential

100
Volt[μV]

Borderline of Effective
Muscle Activity
50

0
0

2

4

6

8

10
12
14
16
Time[sec]
Biceps Brachii Myogenic Potential (Post-Processing)

18

20

Volt[μV]

100
Borderline of Effective
Muscle Activity
50

0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

Triceps Brachii Myogenic Potential (Post-Processing)Fig.4.1.6 (1) Experiment Results of Curvature Movement in Case of Moved Arm Position43Force[N]

0
-5 0
-10
-15
-20
-25
-30
-35

2

4

6

8

10

12

14

16

18

20

12

14

16

18

20

14

16

18

20

Time[sec]

Muscle Tension[N]Force of Arm Tip

150
100
50
0
0

2

4

6

8

10
Time[sec]

Muscle Tension[N]Estimated Biceps Brachii Tension

100
50
0
0

2

4

6

8

10
Time[sec]

12

Estimated Triceps Brachii Tension

Fig.4.1.6 (2) Experiment Results of Curvature Movement in Case of Moved Arm Position

44

Position[m]0-0.005 0-0.01-0.015-0.02-0.025-0.03-0.035

2

4

6

8

10

12

14

16

18

20

14

16

18

20

14

16

18

20

Time[sec]
Position of Arm Tip

0.008
Angle[rad]

0.006
0.004
0.002
0
-0.002

0

2

4

6

8

10

12

Time[sec]
Shoulder Joint Angle

Angle[rad]

1.72
1.68
1.64
1.6
1.56
0

2

4

6

8

10
Time[sec]

12

Elbow Joint Angle

Fig.4.1.6 (3) Experiment Results of Curvature Movement in Case of Moved Arm Position45

Volt[μV]

100

50

0
0

2

4

6

0

2

4

6

8

10
12
14
Time[sec]
Biceps Brachii Myogenic Potential

16

18

20

8

16

18

20

Volt[μV]

250
200
150
100
50
0
10
12
14
Time[sec]
Triceps Brachii Myogenic Potential

Volt[μV]

100
Borderline of Effective
Muscle Activity
50

0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

Biceps Brachii Myogenic Potential (Post-Processing)250
Borderline of Effective
Muscle Activity

Volt[μV]

200
150
100
50
0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

Triceps Brachii Myogenic Potential (Post-Processing)Fig.4.1.7 (1) Experiment Results of Extension Movement in Case of Moved Arm Position46

Force[N]

30
25
20
15
10
5
0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

14

16

18

20

10
12
14
Time[sec]
Estimated Triceps Brachii Tension

16

Muscle Tension[N]Force of Arm Tip

100
50
0
0

2

4

6

8

10
Time[sec]

12

Muscle Tension[N]Estimated Biceps Brachii Tension

300
250
200
150
100
50
0
0

2

4

6

8

18

20

Fig.4.1.7 (2) Experiment Results of Extension Movement in Case of Moved Arm Position

47

Position[m]

0.03
0.025
0.02
0.015
0.01
0.005
0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

14

16

18

20

14

16

18

20

Position of Arm Tip

0.005
Angle[rad]

0.004
0.003
0.002
0.001
0
-0.001 0

2

4

6

8

10

12

Time[sec]

Angle[rad]Shoulder Joint Angle

1.58
1.56
1.54
1.52
1.5
1.48
1.46
1.44
0

2

4

6

8

10
Time[sec]

12

Elbow Joint Angle

Fig.4.1.7 (3) Experiment Results of Extension Movement in Case of Moved Arm Position48

Table4.2.1 Specification of Direct Drive（DD） Robot
2 軸 DD ロボット SR-402DD(S)（東京エレクトロニックシステムズ製）＜アーム機構部仕様＞
自由度
2
駆動方式
DD モータアーム長…200+200[mm]寸法
全高…約 250[mm]
本体重量
約 40[kg]
動作範囲
1 軸…±90[deg]，2 軸…±120[deg]最大先端力（伸張時）
71[N]最大動作速度
1 軸…720[deg/s]，2 軸…864[deg/s]
最大合成速度（無負荷時）
8.0[m/s]＜モータ部仕様＞
型名
最大出力トルク
最大回転数
定格回転数
エンコーダ分解能
ロータイナーシャ
重量

1軸

2軸

DR1070E
70[Nm]
2.0[rps]
1.5[rps]614,400[p/rev]
85×10-3[kgm]
22[kg]

DR1015B
15[Nm]
2.4[rps]
2.0[rps]507,904[p/rev]
21×10-3[kgm]8[kg]Fig.4.2.1 Experiment Scene

49

PC

Force

6 axesForce/TorqueController

Control Signal

SensorCushion

Arm Position

Robot Arm PA-10xMoving Direction

θy
2 axes DD

Direction of the Position &Robot

Force Coordinate
PC

Control Signal.
Controller
Data

Fig.4.2.2 System Structure of Control Method Experiment50

StartInput Vn in Robot Arm.
Robot Arm moved.

dT,M,K,C,
link,X0,Xn
Read Xn from
Robot Arm

Loop Start
X0=X0+D
Vo=Vn

Calculate Jrad by using
Xn and link.

Yes
Jrad>radNo
Read F from
Force Sensor.

Yes
Push any Key?

No

dX=Xn-X0

Loop End
Vn=dt/M*(F-C*Vo-K*dX)+Vo

End

Fig.4.2.3 Flowchart of Virtual Compliance Control Program51

StartInput Vn in Robot Arm.
Robot Arm moved.

dT,M,K,C,
link,X0,Xn,
VlimitRead Xn from
Robot Arm

Loop Start

X0=X0+D
Vo=Vn
Calculate Jrad by using
Xn and link.

Yes
Jrad>rad
Read F from
Force Sensor.

No
Yes
Push any Key?

dX=Xn-X0

No
Vn=dt/M*(F-C*Vo-K*dX)+Vo

Vn<Vlimit

Loop End

Yes
End

No
Vn=Vlimit

Fig.4.2.4 Flowchart of VCL and PVCL Program52

Start1

dT,M,K,C,
link,X0,Xn,Xv,
Vlimit,Fa

X0=X0+D
Xv=Xv+D

Loop Start
Vn<Vlimit

2

Yes

No

Calculate Jrad by using
Xn and link.

Vn=Vlimit

Read F from
Force Sensor.

Input Vn in Robot Arm.
Robot Arm moved.

dX=Xn-X0
DX=Xv-X0
Read Xn from
Robot Arm
Vn=dT/M*(F-C*Vo-K*dX)+Vo
VXn=dT/M*(F-C*VXo-K*DX)+VXo

F>Fa &
Vn<Vlimit

Vo=Vn
VXo=VXn

YesJrad>radNo

Xv=Xv+VXn*dT
X0=X0+D

Vn>Vlimit

Yes

No

Push any Key?

2

Yes

Yes

No
No

Xv=Xv+D
X0=Xv
dXn=Xn-X0
Vn=dT/M*(F-C*Vo-K*dX)+Vo

1

Loop EndEnd

2

Fig.4.2.5 Flowchart of SL Program53

Start1

dT,M,K,C,
link,X0,Xn,Xv,
E,Vlimit,Fa

X0=X0+D
Xv=Xv+D

Loop Start
Vn<Vlimit

2

Yes

No

Calculate Jrad by using
Xn and link.

Vn=Vlimit

Read F from
Force Sensor.

Input Vn in Robot Arm.
Robot Arm moved.

dX=Xn-X0
DX=Xv-X0
Read Xn from
Robot Arm
Vn=dT/M*(F-C*Vo-K*dX)+Vo
VXn=dT/M*(F-C*VXo-K*DX)+VXo

F>Fa &
Vn<Vlimit

Vo=Vn
VXo=VXn

YesJrad>radNo

Xv=Xv+VXn*dT
X0=X0+E*dT/F

Vn>Vlimit

Yes

Yes

No

Push any Key?

2

Yes

No
No

Xv=Xv+D
X0=Xv
dXn=Xn-X0
Vn=dT/M*(F-C*Vo-K*dX)+Vo

1

Loop EndEnd

2

Fig.4.2.6 Flowchart of SEL Program54

Angle[rad]2
Encoder Output
Estimated Angle

1.5
1
0.5
0
0

2

4

6

8
10
Time[sec]

12

14

16

18

12

14

16

18

12

14

16

18

12

14

16

18

Joint Angle

10

Force[N]

0
-10 0

2

4

6

8

10

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

Current Position
Target Position

0

2

4

6

8
10
Time[sec]
X Direction Position

0.04

Position[m]0.03

Current Position
Target Position

0.02
0.01
0
-0.01 0

2

4

6

8

10

-0.02

Time[sec]
Y Direction Position

Fig.4.2.7 (1) Experiment Results of Virtual Compliance Control55

0Distance[m]

-0.01

0

2

4

6

8

10

12

14

16

18

16

18

-0.02
-0.03
-0.04
-0.05
Time[sec]
Distance between Current Position and Target Position on X Direction

Distance[m]

0.005
0
0

2

4

6

8

10

12

14

-0.005-0.01-0.015
Time[sec]
Distance between Current Position and Target Position on Y Direction

Velocity[m/sec]

0.2
0.15
0.1
0.05
0
-0.05

0

2

4

6

8

10

12

14

16

18

12

14

16

18

Time[sec]
X Direction Velocity

Velocity[m/sec]

0.08
0.06
0.04
0.02
0
-0.02 0

2

4

6

8

10

-0.04
Time[sec]
Y Direction Velocity

Fig.4.2.7 (2) Experiment Results of Virtual Compliance Control

56

Angle[rad]2
Encoder Output
Estimated Angle

1.5
1
0.5
0
0

2

4

6

8
10
Time[sec]

12

14

16

18

12

14

16

18

12

14

16

18

12

14

16

18

Joint Angle

10

Force[N]

0
-10 0

2

4

6

8

10

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

Current Position
Target Position

0

2

4

6

8
10
Time[sec]
X Direction Position

0.04
Current Position
Target Position

Position[m]

0.03
0.02
0.01
0
-0.01

0

2

4

6

8

10

-0.02
Time[sec]
Y Direction Position

Fig.4.2.8 (1) Experiment Results of VCL
57

0Distance[m]

-0.01

0

2

4

6

8

10

12

14

16

18

-0.02
-0.03
-0.04
-0.05
Time[sec]
Distance between Current Position and Target Position on X Direction

Distance[m]

0.005
0
0

2

4

6

8

10

12

14

16

18

-0.005-0.01-0.015
Time[sec]
Distance between Current Position and Target Position on Y Direction

Velocity[m/sec]

0.04
0.03
0.02
0.01
0
-0.01

0

2

4

6

8

10

12

14

16

18

12

14

16

18

Velocity[m/sec]

Time[sec]
X Direction Velocity

0.03
0.02
0.01
0
-0.01 0
-0.02
-0.03
-0.04

2

4

6

8

10

Time[sec]
Y Direction Velocity

Fig.4.2.8 (2) Experiment Results of VCL
58

Angle[rad]2
Encoder Output
Estimated Angle

1.5
1
0.5
0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

12

14

16

18

20

12

14

16

18

20

14

16

18

20

Joint Angle

10
Force[N]

0
-10

0

2

4

6

8

10

-20
-30
-40
Time[sec]
Force

0.5
Position[m]

0.4
0.3

Current Position
Target Position

0.2
0.1
0
-0.1 0

2

4

6

8

10

Position[m]Time[sec]
X Direction Position

0.06
0.05
0.04
0.03
0.02
0.01
0
-0.01 0Current Position
Target Position

2

4

6

8

10
Time[sec]

12

Y Direction Position

Fig.4.2.9 (1) Experiment Results of PVCL59

0Distance[m]-0.02 0

2

4

6

8

10

12

14

16

18

20

18

20

-0.04
-0.06
-0.08
-0.1
-0.12
Time[sec]
Distance between Current Position and Target Position on X Direction

Distance[m]

0.01
0
0

-0.01

2

4

6

8

10

12

14

16

-0.02
-0.03
-0.04
Time[sec]
Distance between Current Position and Target Position on Y Direction

Velocity[m/sec]

0.04
0.03
0.02
0.01
0
-0.01

0

2

4

6

8

10

12

14

16

18

20

14

16

18

20

Time[sec]
X Direction Velocity

Velocity[m/sec]

0.03
0.02
0.01
0
-0.01

0

2

4

6

8

10

12

-0.02
Time[sec]
Y Direction Velocity

Fig.4.2.9 (2) Experiment Results of PVCL
60

Angle[rad]2

Encoder Output
Estimated Angle

1.5
1
0.5
0
0

2

4

6

8

10
Time[sec]

12

14

16

18

20

12

14

16

18

20

10
12
Time[sec]
X Direction Position

14

16

18

20

8

14

16

18

20

Joint Angle

10

Force[N]

0
-10 0

2

4

6

8

10

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

Current Position
Target Position
Virtual Target Position

0

2

4

6

8

Position[m]0.06
Current Position
Target Position
Virtual Target Position

0.04
0.02
0
0

2

4

6

10

12

-0.02
Time[sec]
Y Direction Position

Fig.4.2.10 (1) Experiment Results of SL

61

0.01
Distance[m]0-0.01 0

2

4

6

8

10

12

14

16

18

20

18

20

-0.02
-0.03
-0.04
-0.05
Time[sec]
Distance between Current Position and Target Position on X Direction

Distance[m]

0.005
0
0

2

4

6

8

10

12

14

16

-0.005-0.01-0.015Velocity[m/sec]

Time[sec]
Distance between Current Position and Target Position on Y Direction

0.03
0.025
0.02
0.015
0.01
0.005
0
-0.005 0

2

4

6

8

10

12

14

16

18

20

14

16

18

20

Time[sec]

Velocity[m/sec]

X Direction Velocity

0.03
0.02
0.01
0
-0.01 0
-0.02
-0.03
-0.04

2

4

6

8

10

12

Time[sec]
Y Direction Velocity

Fig.4.2.10 (2) Experiment Results of SL62

Angle[rad]2

Encoder Output
Estimated Angle

1.5
1
0.5
0
0

10

20

30

40

50

60

70

Time[sec]
Joint Angle

10

Force[N]

0
-10 0

10

20

30

40

50

60

70

40
Time[sec]
X Direction Position

50

60

70

50

60

70

-20
-30
-40
-50
Time[sec]

Position[m]

Force

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

Current Position
Target Position
Virtual Target Position

0

10

20

30

Position[m]0.06
Current Position
Target Position
Virtual Target Position

0.04
0.02
0
0

10

20

30

40

-0.02
Time[sec]
Y Direction Position

Fig.4.2.11 (1) Experiment Results of SEL

63

0.01
Distance[m]0-0.01 0

10

20

30

40

50

60

70

-0.02
-0.03
-0.04
-0.05
Time[sec]
Distance between Current Position and Target Position on X Direction

Distance[m]

0.005
0
0

10

20

30

40

50

60

70

-0.005-0.01-0.015
Time[sec]
Distance between Current Position and Target Position on Y Direction

Velocity[m/sec]

0.025
0.02
0.015
0.01
0.005
0
-0.005 0

10

20

30

40

50

60

70

50

60

70

Time[sec]
X Direction Velocity

0.015
Velocity[m/sec]

0.01
0.005
0
-0.005 0

10

20

30

40

-0.01
-0.015-0.02
Time[sec]
Y Direction Velocity

Fig.4.2.11 (2) Experiment Results of SEL

64

Myogenic Potential

Estimated Muscle Tension

Normalized Outputs

1
0.8
0.6
0.4
0.2
0
0

5

10
Time[s]15

20

Fig.5.1.1 Experiment Result of Biceps Brachii in Case of Moved Arm Position

Myogenic Potential

Estimated Muscle Tension

Normalized Outputs

1
0.8
0.6
0.4
0.2
0
0

5

10
Time[s]15

20

Fig.5.1.2 Experiment Result of Triceps Brachii in Case of Moved Arm Position

65

・付録
パラメータT ， J の導出
A パラメータ
本研究では，岩本(9)が構築した 2 自由度水平面運動学モデルを参考として，上肢のモデルを作成した．ここでは，本文で述べた腕先力と各筋張力の関係で用いられるパラメータ T ， J を導出する．

A-1 筋の変位と関節角度の関係
各筋の変位量 ei （i＝BS，BR，AN，TBLa，TBLo，BB，DA，DP，PC）は，関節角度の関数
であり，肩関節角度，肘関節角度をそれぞれ θ s ， θ e とすると，次のように表される．
上腕筋（BS）

e BS = r11θ e + r12θ e2 + r13θ e3 + r14θ e4（A.1）

腕橈骨筋（BR）

e BR = r21θ e + r22θ e2 + r23θ e3 + r24θ e4（A.2）

肘筋（AN）e AN = r31θ e + r32θ e2 + r33θ e3 + r34θ e4 + r35θ e5 + r36θ e6（A.3）

上腕三頭筋外側頭（TBLa）

eTBLa = r41θ e + r42θ e2 + r43θ e3 + r44θ e4 + r45θ e5 + r46θ e6（A.4）

上腕三頭筋外側頭（TBLo）

eTBLo = r51θ e + r52θ e2 + r53θ e3 + r54θ e4 + r55θ e5 + r56θ e6 + r5θ s（A.5）

上腕二頭筋（BB）

e BB = r61θ e + r62θ e2 + r63θ e3 + r64θ e4 + r6θ s（A.6）

三角筋前部（DA）

e DA = r7θ s（A.7）

三角筋後部（DP）

e DP = r8θ s（A.8）

大胸筋鎖骨部（PC）

e PC = r9θ s（A.9）

ここでは，筋の変位は筋が縮小するときに正の値をとり，筋が伸張するときは負の値を取るものとする．また，各式における定数 rij については，岩本の作成したモデルより引用した．定数の数66

値を TableA.1 に示す．

A-2 パラメータT
パラメータ の導出
A-1 節で示した式を行列化すると

 e BS   0
e  
 BR   0
 e AN   0 eTBLa   0
eTBLo  =  r5  e BB  r6
 e  r
 DA   7
 e DP   r8
 e  r
 PC   9


r11 + r12θ e + r13θ e2 + r14θ e3
23
r21 + r22θ e + r23θ e + r24θ e

2
3
4
5
r31 + r32θ e + r33θ e + r34θ e + r35θ e + r36θ er41 + r42θ e + r43θ e2 + r44θ e3 + r45θ e4 + r46θ e5 θ s 
r51 + r52θ e + r53θ e2 + r54θ e3 + r55θ e4 + r56θ e5    θ e 
r61 + r62θ e + r63θ e2 + r64θ e3


0


0

0


と表せる．筋の変位ベクトルを e = [e BS
節角度ベクトルを Θ =[θ s

e BR

e AN

eTBLa

eTBLo

e BB

（A.10）e DA

e DP

e PC ] ，関T&&e ，関節角の微小変位ベクト
θ e ]T とする．筋の微小変位ベクトルを A&&Θ とすると， A&&e と A&&Θ の関係は
ルを A&&e = T ⋅ A&&ΘA（A.11）&&Θ から A&&e への変換行列であり，次式のようになる．
となる．ここで， T は（9×2）の行列で， A
0

0
0

0
T =  r5

 r6
r
7
 r8

 r9


r11 + 2r12θ e + 3r13θ e2 + 4r14θ e3
23
r21 + 2r22θ e + 3r23θ e + 4r24θ er31 + 2r32θ e + 3r33θ e2 + 4r34θ e3 + 5r35θ e4 + 6r36θ e5 
r41 + 2r42θ e + 3r43θ e2 + 4r44θ e3 + 5r45θ e4 + 6r46θ e5 
r51 + 2r52θ e + 3r53θ e2 + 4r54θ e3 + 5r55θ e4 + 6r56θ e5 
r61 + 2r62θ e + 3r63θ e2 + 4r64θ e3


0


0

0


（A.12）

A-3 パラメータ J の導出
手先の位置ベクトルを r = [xy ] ，関節角度ベクトルを Θ =[θ s θ e ]T とする．それぞれの微T&&r ， A&&Θ とすると，両者の関係は次式となる．
小変位ベクトルを A

67

&&r = J ⋅ A&&ΘA（A.13）

よって，ヤコビ行列 J は，式（A.13）と Fig.2.1.1 より，

− l sin θ s − l 2 sin (θ s + θ e ) − l 2 sin (θ s + θ e )J = 1 l1 cosθ s + l 2 cos(θ s + θ e ) l 2 cos(θ s + θ e ) となる． l1 ， l 2 は，上腕骨と前腕骨の長さである．

68（A.14）

B 無次元化された等尺性最大張力の数式モデル
等尺性最大張力の計算において，筋の歪みの影響を表すパラメータとして無次元化した等尺性
最大張力が用いられる．この無次元化した等尺性最大張力は，無次元化した能動的最大張力と，無次元化した静止張力の和として求められる．初めに，無次元化した等尺性最大張力のパラメータである筋の構造指標 ia と筋の歪み h について説明する．次に無次元化した能動的最大張力と，無次元化した静止張力の数式モデルを示し，最後に無次元化された等尺性最大張力の数式モデルを示す．

B-1 筋の構造指標
筋の構造指標 ia とは，筋繊維長 fl と筋節長 ml の比率であり，
ia = fl ml（B.1）と表される．この指標を用いる事で，筋の長さと能動的最大張力の関係を一般的に表現できる(13)．
また，この指標の変化により，筋の長さと静止張力の関係を表す特性曲線がある傾向を示すことが報告されている．

B-2 筋の歪み
筋の収縮によって腕の運動が発生する事から，筋の歪みは関節角度の関数と言える．Piegon(10)らの筋の変位量 ei (θ s )（i＝DA，DP，PC），ei (θ e )（i=BS，BR，AN，TBLa），ei (θ s ,θ e )（i=TBLo，BB）は筋が収縮するときに，正の値をとる．よって，筋の歪みは，式（A.1）から式（A.9）と各筋の自然長 ml i を用いるとe (θ ) − ei (θ sn )hi (θ s ) = − i s
mli

（B.2）

e (θ ) − ei (θ en )hi (θ e ) = − i e
ml i

（B.3）

e (θ ,θ ) − ei (θ sn ,θ en )hi (θ s ,θ e ) = − i s e
mli

（B.4）

と表される． θ sn ， θ en は各筋が自然長となるときの関節角度を表す．

B-3 無次元化した能動的最大張力
Kaufman(13)らは，自然長時の能動的最大張力を１として，筋の歪みと，無次元化した能動的最大張力 fˆa (h ) の関係を次のように表した
2{0 . 96343 [1−1 i a ]}
− 1 .0  
ˆf (h ) = exp  −  (h + 1)


a
 
0 . 35327 (1 − i a )
 
 

69

ia < 1

（B.5）

[

2
fˆa (h ) = exp − {2.727 ⋅ ln (h + 1)}

ただし]

ia = 1F (h )fˆa (h ) = a
Fmax

（B.6）

（B.7）

である．上式において， Fa (h) は任意の歪みにおいて筋が発揮し得る能動的最大張力， Fmax は最大張力， h は筋の歪み， ia は筋の構造指標を表す．

B-4 無次元化した静止張力
岩本は，式（B.5）を参考に，筋の歪み h と， Fmax で無次元化した静止張力 fˆ p (h ) の関係を次のように表した2 ˆf (h ) = 0.563exp2.239 h   − 1.0
p
i 
 

 a  



h≥0

（B.8）

fˆ p (h) = 0

h<0

（B.9）

Fp (h )fˆ p (h ) =Fmax（B.10）

上式において， F p (h) は任意の歪みにおいて筋が発揮し得る静止張力， Fmax は最大張力， ia は筋の構造指標を表す．

B-5 無次元化した等尺性最大張力
等尺性最大張力は，B-3 節で示した能動的最大張力と，B-4 節で示した静止張力の和で求める事ができる．よって，筋の歪み h と， Fmax で無次元化した等尺性最大張力の関係を数式化すると，h ≥ 0 のときは
2
2{0.96343[1−1 ia ]} h  − 1 .0  
ˆf (h ) = exp  −  (h + 1)
  + 0.563exp 2.239   − 1.0
i 
 
0.35327 (1 − ia )


 

 a   ia < 1（B.11）2ˆf (a ) = exp − {2.727 ⋅ ln (a + 1)}2 + 0.563exp 2.239 a   − 1.0
i 



 a  
[]ia = 1（B.12）

h < 0 のときは70

2{0 . 96343 [1−1 ia ]}− 1 .0  
ˆf (h ) = exp  −  (h + 1)
  
0 .35327 (1 − i a )
 
 

[

2
fˆ (h ) = exp − {2.727 ⋅ ln (h + 1)}]

ia < 1（B.13）

ia = 1（B.14）となる．式（B.11）から式（B.14）までを用いて，筋の歪み h と無次元化した等尺性最大張力fˆ (h ) の関係を， ia をパラメータとして 0.25，0.50，0.75，1.00 と変化させて Fig.B.1 に示す．

71

用語解説
spasticity）
・痙縮（
）
・痙縮
（spasticity
筋肉が収縮された状態で保たれ，伸展に対して過度に反応する状態．
（感度が増加している）
痙縮は人間の中枢神経系が損傷した後に見られる最もありふれた病気の１つである．臨床的には痙縮は受動動作に対する抵抗として特徴づけられ，この抵抗は筋が高い張力に達するまで伸展す
ればする程強度が増す．このとき，筋はいわゆる 折りたたみナイフ（ジャックナイフ） 反応で突然弛緩する．受動動作に対するこの抵抗は位置ばかりでなく速度変化にも敏感で，手足が早く伸展すれば，早く生成される．

・収縮要素
収縮要素とは，最も一般的な筋肉モデルである 3 要素モデル(11)内の 1 つである．収縮要素が発生する張力は能動的張力と呼ばれ，その限界値は筋肉の長さによって変化する．筋肉が自然長となるときの筋肉の歪みを 0 とすると，筋肉の歪みと能動的最大応力の関係は Fig.2.1.2 に示すような曲線となる．図より能動的最大張力は，筋肉の歪みが 0 のとき最大値をとり，筋肉が自然長よ
り長くても，短くても，値が小さくなる事がわかる．

・並列弾性要素
並列弾性要素とは，最も一般的な筋肉モデルである 3 要素モデル(11)内の 1 つである．筋肉を自然長以上に伸ばすと，Fig.2.1.2 に示すような静止張力（受動的張力）が発生する．これは並列弾性要素が発揮する張力として説明される．並列弾性要素はフックの法則に従わない非線型なバネの性質を持ち，発生する静止張力は筋肉の長さの関数となる．

・表面筋電位計
体表にある筋肉の筋電位を計測する機器で，その結果から表面筋電図が作られる．表面筋電図
からは筋活動の有無やその強さ，活動の開始時刻や持続時間，あるいは動作に伴う複数の筋肉相互間の協調性や相反性を知る事ができる．こうした表面筋電図の分析は，運動の分析や作業における筋負担の人間工学的評価に，あるいはリハビリテーション訓練において回復課程の評価などに用いられている．

72

TableA.1 Parameters of Human Arm Model
Elbow

Muscle [ j ]
BS[1]
BR[2]
AN[3]
TBLa[4]
TBLo[5]
BB[6]
DA[7]
DP[8]
PC[9]

rj1

rj2

rj3

rj4

5.549
19.490-5.345-23.287-23.287
14.660

6.612
4.779-6.544-8.676-8.67612.9842.563
11.0359.224
14.10114.101
1.975-0.965-3.065-6.738-8.978-8.978-1.405

Shoulder
rj5

rj6

2.252
2.862
2.862

-0.281-0.362-0.362sj-25.400
29.21033.020-78.740
50.800Normalized Muscle Tension

1.4
1.2
1
0.8
0.6
i=0.25i=0.5
i=0.75
i=1

0.4
0.2
0
-1

-0.5

0
Muscle Strain

0.5

1

Fig.B.1 Relationship of Normalized Muscle Tension and Differing Architecture73