筋疲労の評価方法、筋疲労度評価装置、および、使用者の生理学的状況をリアルタイムで反映する運動支援システム
【課題】運動時の関節角度や床反力などの人体外部の力学的諸量を用いて人体内部の状況を推定し、筋疲労度を評価する方法、そのための装置及びそれらを用いた安全且つ効果的な運動支援システムを提供すること。
【解決手段】使用者に装着したセンサーを介して得られた関節角度θ及び床反力Ｆを用いて筋張力を計算し、次いで得られた筋張力を前記センサーを用いた測定開始から測定終了までの時間範囲で積分する筋疲労度の評価方法、そのための装置及びそれらを用いたシステム。
前記計算は下記一般式(1)に基づいてなされる。
一般式(1)
Fce=k(Vce,Lce)・(q-αffatigue(∫Fcedt))
但し、式中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉ごとに決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式である。
【解決手段】使用者に装着したセンサーを介して得られた関節角度θ及び床反力Ｆを用いて筋張力を計算し、次いで得られた筋張力を前記センサーを用いた測定開始から測定終了までの時間範囲で積分する筋疲労度の評価方法、そのための装置及びそれらを用いたシステム。
前記計算は下記一般式(1)に基づいてなされる。
一般式(1)
Fce=k(Vce,Lce)・(q-αffatigue(∫Fcedt))
但し、式中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉ごとに決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式である。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はコンピュータを用いた筋肉疲労度の評価、及びその結果を用いた使用者の生理学的状況をリアルタイムで反映する運動支援システムに関し、特に、使用者の運動部位と対応する関節角度及び床反力データから、運動による使用者の筋肉疲労度をコンピュータを用いて算出する方法、そのための装置、それを利用した使用者の生理学的状況をリアルタイムで反映する運動支援システムに関する。
【背景技術】
【０００２】
従来、リハビリテーションは理学療法士によって、助言・指導がなされているが、理学療法士の負担が大きいために、高齢化社会、高福祉社会を迎えた今日、需要に充分答えることができないというのが実状である。そこで、従来の健康増進用機器のみならず、更に高度な技術が要求されるリハビリテーション用等の、医療機器の開発が望まれている。即ち、健康増進用機器の場合には、機器に要求される性能の任意性を広く許容することができるが、健常者とは異なる弱者に対する機器の場合には、使用者の状態に機器の方が対応できることが望まれるとともに、安全性に対する要求が一段と厳しいものとなる。
【０００３】
従来、安全性を確保する観点から、臥位で動作を行わせることによって転倒を防止したり、ハーネスなどを用いて体の自由度をある程度拘束する方法が採られており(特許文献1-4)、機器自体が大型化せざるを得ないという欠点があった。また、これまでは、人体内部の状況に対応して機器の方を調整するという発想が提案されることはなかった。しかしながら使用者の人体内部の情報を容易に推定することができ、推定した情報を機器にリアルタイムで反映することができれば、リハビリテーションなどの効果を最大にすることができる筈である。
【特許文献１】特開2000-102576号公報
【特許文献２】特開2005-192695号公報
【特許文献３】特開2005-211086号公報
【特許文献４】特開2005-211328号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
従って本発明の第一の目的は、運動時の関節角度や床反力などの人体外部の力学的諸量を用いて人体内部の状況を推定し、筋疲労度を評価する方法を提供することにある。
本発明の第二の目的は、運動時の関節角度や床反力などの人体外部の力学的諸量を用いて、人体内部の情報をデータ化すると共に、筋疲労度を評価するための装置を提供することにある。
更に本発明の第三の目的は、使用者の状況をリアルタイムで反映し、使用者の生理学的状況をリアルタイムで反映させることにより、無理なく安全であると共に効果的に運動を行うことのできる運動支援システムを提供することにある。
【課題を解決するための手段】
【０００５】
本発明者等は上記の諸目的を達成すべく鋭意研究した結果、特定の計算式を用いてコンピュータによる計算を行わせることにより、人体外部の力学的量である関節角度、及び床反力から使用者の筋疲労度を推定することができることを見出し、本発明に到着した。
即ち本発明は、使用者に装着したセンサーを介して得られた関節角度θ及び床反力Fを用いて筋張力を計算し、次いで得られた筋張力を前記センサーを用いた測定開始から測定終了までの時間範囲で積分する筋疲労度の評価方法であって、前記計算が下記一般式(1)に基づいてなされることを特徴とする筋疲労度の評価方法、その評価方法に使用する装置、及び、その評価方法を利用した使用者の生理学的状況をリアルタイムで反映する運動支援システムである。
一般式(1)
本発明においては、前記関節角度として、足関節、膝関節、股関節の各角度を採用することが好ましい。上半身が該当箇所である場合はその部位に至るまでの主要な各関節角度を入力しなければならない。該当箇所が下半身の末端に近ければ膝関節、股関節の角度は必ずしも必要ではない。さらに、人体モデルを詳細にし、入力対象となる関節数を増やすことによってより精緻な結果を得ることができる。また、安全を確保するために、設定した筋疲労度に達したときに運動装置が停止するように、安全装置を設けることが好ましい。
【発明の効果】
【０００６】
本発明によれば、人体外部の力学的量を用いて筋疲労度を推定するので、容易であるにもかかわらず使用者に負担をかけない上、得られたデータをリアルタイムで運動装置に反映させることができるので、安全に且つ最大限の効果を得ることが容易である。
【発明を実施するための最良の形態】
【０００７】
以下、図面を参照しながら本発明を詳述する。図1は本発明の実施例の概念図である。図において、符号1は足、膝、股関節角度及び床反力検出センサーを有し使用者に強制歩行させるための運動具、符合2は、得られた各関節角度および床反力に基づいて使用者の筋張力、筋疲労度を算出するためのコンピュータ、符合3は、算出された筋疲労度を運道具に反映させ、使用者の運動を制御するための制御部である。
【０００８】
本発明においては、使用者の骨格間の角度を測定し、これらの角度を、コンピュータ内にとりこまれた人間の骨格を有する人体モデルに入力して多体系問題を解き、次いで得られた関節モーメントを下記一般式(1)を含む筋モデルに適用し、最適化法を用いて使用者の筋張力決定し、最後にこの筋張力から筋疲労度を算出する（図2参照）。
一般式(1)
【０００９】
例えば、歩行時における下肢の筋肉を対象の部位とした場合には、足、膝および股関節の各角度を指定することによって、人間の骨格を有するコンピュータ内の人体モデルの姿勢を一義的に決定することができる。従って、使用者が装着する運動具に取り付けたセンサーによって使用者の各関節角度θfoot，θknee，θhipを検出し、これをコンピュータに入力すれば、コンピュータ内の人体モデルは使用者と同じ姿勢をとることになる。尚、上記のセンサーとしては、公知のセンサーの中から適宜選択して使用することができる。
【００１０】
コンピュータ内における人体モデル各部位の質量は使用者に合わせて設定され、人体モデルにかかる重力は自動的に計算されるため、人体モデルは使用者と同様の負荷を受ける。これに加え、使用者の体が外部と接触する点から受ける外力を入力すれば、人体モデルの力学的状態は使用者と同じになる。歩行の場合、運道具を含む使用者が外部と接触する点は足の裏と床の接触点のみであるから、左右の床反力FrおよびFlを入力することによって、使用者の力学的状態をコンピュータ2中に再現することができる。
【００１１】
コンピュータ2中の人体モデルは使用者の運動を反映するので、各タイムステップにおける各筋張力および筋疲労度を出力することができる。図1においては結果を出力する出力画面を有しているが、更にプリントアウトするためのプリンターを有していても良いことは当然である。各関節角度θ、床反力Fから筋張力、筋疲労度を計算するシステムフローを図2に示す。
【００１２】
(人体モデルについて)
人体各部位を適当な質量と慣性モーメントをもった剛体と仮定し、各関節を、同等の自由度と機能を持った機械的関節とみなすと、人体を一連の剛体の集まりとして数学的に表すことができ、運動時の各関節に発生する力及びトルクを得る問題は、多体系に対する逆動力学問題に帰着する。従って、各部位の質量及び慣性モーメントを非侵襲的に実測することは不可能であるが、過去の研究と統計学的アプローチによれば、身長及び体重から推定値を得ることが可能である（非特許文献１）。
【非特許文献１】阿江通良,横井孝志:日本人アスリートの身体部分慣性特性の推定 バイオメカニズム11,23-33東京大学出版会
【００１３】
また、機構学の分野においては、各関節における任意の機能を数学的に表現することがなされてきた。例えば、人体の肘や膝関節は、歩行や投擲などの大変位の運動においては、単一回転自由度のみを持つ回転関節と仮定しても充分にその運動時の挙動を推定することができることが知られている（非特許文献２）。
【非特許文献２】山崎信寿，長谷和徳:自由歩行における歩調・歩幅の生体力学的決定基準 バイオメカニズム11,179-189東京大学出版会
【００１４】
(拘束式の例 回転拘束,加速度式)
二次元平面で軸周りの回転だけを許容する回転関節は、数学的に
と表すことができる。但し、式中のｒは各剛体の並進座標、sは剛体の中心から関節の拘束中心までのベクトル、添え字iおよびjは関節によって結合される二つの剛体をそれぞれ表す。回転関節によって結合される物体iおよびjは、常にこの式を満たす範囲で運動せねばならない。
【００１５】
また、この式を時間で偏微分することにより、速度に関する拘束式（速度式）を得ることができる。速度式は回転関節によって結合される剛体が、運動中に満たさなければならない速度の方程式である。また速度式をさらに偏微分することによって、同様に加速度式を得ることができる。
【００１６】
三次元の回転拘束式は
と表される。式中のΦｓは二次元の回転拘束式と同様のものであり、剛体iおよびj上の関節中心PiおよびPjが同一点であることを意味している。またΦplはi、j双方で定義される回転の軸hiおよびhjが常に並行であることを意味している。概念図を図3に示す。具体的には、ΦｓおよびΦplのそれぞれは下記のように表される。
上式中、Aは各剛体からグローバル座標系へ変換する座標変換マトリクス、ｈは図に示される関節の回転方向を定義するベクトル、fおよびgは、P点を原点としベクトルhと直行する、関節座標系をなすベクトルである。
【００１７】
関節の拘束式に加え、時間ごとの関節角度を指定することによって、各時間における関節の状態を一義的に決定することができる。入力する関節角度の情報は、ドライバ拘束といわれるものであり、
で表される。三次元であっても、拘束式を時間微分することによって、速度式及び加速度式を得ることができる。
実際のプログラム内では、拘束式、速度式および加速度式は代数学的にライブラリ化されており、関連する物体の状態量を入力することによって計算され、時間微分を行わず直接計算することにより、高速に処理される。
【００１８】
(運動方程式)
人体の各関節を拘束式Φi(i=0,1,Λnj)(njは関節の数)で表すと、人体全体の
Newton-Euler運動方程式は
と表すことができる。式中のMは、人体の各部位の質量を対角成分にもつ対角行列、J'は各部位のローカル座標系周りの慣性モーメントを対角成分にもつ対角行列であり、「’」はローカル座標系を対象としていることを表している。ΦｒおよびΦπ'は、それぞれ並進自由度に関する拘束式およびEuler角で表したローカル座標系の回転自由度に関する拘束式である。
【００１９】
左辺の二つめの[]は変数行列であり、ｒは各部位の並進方向の加速度、ω'はローカル系での角加速度、λは運動方程式と加速度式を結びつけるLagrange乗数である。逆動力学解析においてはこのλを求めることにより、各関節に発生する力及びトルクを決定することができる。
【００２０】
右辺の行列におけるFΑは一般化外力を表す。一般化外力とは、系に作用する外力を、作用する各部位のローカル座標系で表したものである。二項目は一般化トルクを表しており、N'Αは外部から作用する外トルクである。またγは、各拘束式から得られる前述の加速度式である。
【００２１】
(解法について)
上述の運動方程式は、左辺の係数マトリクスが対称行列であり、左辺のベクトルはモデルに入力された諸条件から決定されるため、ガウス消去法、LU分解法などによって解くことできる。そして、ラグランジェ乗数ベクトルλが決定されると、各関節に生じる拘束反力および反トルクは、それぞれ
と表すことができる。式中の「”」は、これらの式が関節座標系での表現であることを表している。また、Cは関節座標系から関連する物体への座標変換マトリクス、上付きの「〜」は歪行列を表している。
【００２２】
以上の、人体モデルの構築と運動の指定および解法を用いることにより、運動中の関節反力と関節トルクの計算を行うことができる。この方法はマルチボディ・ダイナミクスに基づいており（非特許文献3）、得られた関節トルクは、複数の筋肉の筋張力によって関節周りに発生するトルクの和である。関節トルクを、適当な筋肉ごとのパラメータを持った筋モデルを用い、最適化法により解くことによって、各筋肉の筋張力を推定できることが知られている（非特許文献４）。
【非特許文献３】E.J.ハウグ:コンピュータを利用した機構解析の基本,1996,大河出版
【非特許文献４】長谷和徳,山崎信寿(1995):汎用三次元筋骨格モデルの開発 日本機械学会論文集(C編) 61巻
【００２３】
（筋モデルについて）
人体モデルから得られた関節トルクを入力したときにその関節周りの各筋が発生する筋張力を、適当な筋モデルを用い、最適化法によって推定することができる。下記一般式(1)は本発明者等が用いている疲労特性を持つ筋モデルである。
式中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋ごとのパラメータである。kは筋肉の速さ−力関係式であり、収縮要素の収縮時には
FACTOR=Min(1,3.33・q)
と表される。ここでArelおよびBrelはそれぞれ0.4および5.0前後の値をとり、Vceは収縮要素速度を収縮要素長さで割った値である。
【００２４】
収縮要素の伸張時における力と速さの関係は、その低伸張速度部においては、
と表現することができる。ここで、Fasympt（発揮張力の漸近レベル）並びにSlopefactor（速度ゼロにおける力-速さ曲線の伸張部と収縮部の傾きの比）はパラメータであり、それぞれ1.5及び2.0前後の値とすることが妥当である。
【００２５】
一方、高伸張速度部においては
と表現することができる。ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータであり、200前後の値を取る。低速度部と高速度伸張部の境界は
となる点である。以上の筋肉の力学的特性については多くの研究がなされている（非特許文献５〜７）。
【非特許文献５】Hill AV (1938): The heat of shortening andthe dynamic constants of muscle. Proceedings of the Royal Society of London,Series B 126:136-195
【非特許文献６】Felix E. Zajac(1989): Muscle and Tendon:Properties, Models, Scaling And Application to Biomechanics and Motor Control,Critical Reviews in Biomedical Engineering Volume17, Issue 4:359-411
【非特許文献７】Nagano and Gerritsen (2001): Effect ofneuromuscular strength training on vertical jumping performance-a computersimulation study. Journal of Applied Biomechanics 17:113-128
【００２６】
(疲労特性について)
前記一般式（１）中のffatigueは、本発明者等が実験によって導出した表面筋電位-疲労特性であり、
ffatigue=ｍ[(-5.0y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]＋E
と表される。式中のｍは筋肉ごとに決まるパラメータ、ｘは運動開始から現在時点までに筋が発した力の積分値であり(∫Fcedt)、yは運動の周期[秒]を表し、Eは周期が０付近での誤差を少なくするための補正項を表す。
【００２７】
上記ffatigueは、軽負荷における反復運動時に運動条件を変え、疲労に伴う積分筋電位（iEMG）の変化を観察し、得られた筋肉疲労特性のデータを回帰分析することによって決定した、表面筋電位上昇特性式である。この式は、筋張力を計算する際に、運動開始時からの積分値を引数とすると共に、前記一般式(1)式で表されるように、疲労特性を筋肉の収縮要素モデルに加えることによって疲労に伴う表面筋電位の上昇を推定するものであり、本発明者等によって始めて導出された、文献未載の式である。
【００２８】
上記した本発明の筋疲労の評価方法を実施する装置は、使用者に装着させることにより、装着した箇所の関節角度θ及び床反力を測定することのできる手段、前記θ、F及び前記一般式(1)を用いて使用者の筋疲労度を算出する手段、並びに算出された筋疲労度を出力する出力手段とからなる（図1参照）。
【００２９】
上記関節角度θ及び床反力を測定することのできる手段としては、公知のセンサーの中から適宜選択されたセンサーを使用することができ、装置を移動可能とすることもできるが、単なる足踏みとしたり、トレッドミルのような強制歩行装置を用いて非移動式とすることもできる（図6参照）。測定されたθ、Ｆ及び前記一般式(1)を用いて使用者の筋疲労度を算出する手段は、キーボードのようなθやＦを入力するための入力手段、および前記一般式(1)の演算プログラムがインストールされたコンピュータからなり、算出された筋疲労度を出力する出力手段とは、画像として表示することのできる各種ディスプレイ及び/またはプリンターである。
【００３０】
本発明の運動支援システムは、(1)使用者を強制的に運動させる、コンピュータ制御可能なトレッドミルのような運動装置、(2)使用者に装着させ、所望箇所の関節角度θ及び床反力Fを測定する測定手段、(3)得られた前記θ及びFのデータを用いて筋疲労度を算出する手段とを有する運動支援システムである。特に、前記筋疲労度のデータを前記運動装置(1)のコンピュータ制御部に入力することにより、該運動装置の強制運動に使用者の筋疲労度が反映されて制御される。このように、本発明の運動支援システムは使用者の生理学的状況をリアルタイムで反映する運動支援システムである点が最大の特徴である。これによって、本発明の運動支援システムを使用することにより、使用者に無理のない範囲で、安全且つ効率よくリハビリ等の運動を支援することができる。また、安全性を更に高める上で、使用者に許容される疲労度を予め設定し、その値に達したときに自動的に運動システムが停止するように安全装置を備えることが好ましい。
以下本発明を実施例によって更に詳述するが、本発明はこれによって限定されるものではない。
【００３１】
実施例
一般式(1)で表される筋モデル、それを含む人体モデルを用い、下記の様にして歩行時の筋疲労の推定を行った。
トレッドミル上で4km/hの速度で10分間歩行した際の表面筋電位と、その際の歩行に周期をあわせた（T≒1秒）歩行データ（各関節角度、床反力）を、一般式(1)の演算プログラムをインストールしたコンピュータに入力し、シミュレーションによって疲労度を推定した。上記表面筋電位の計測を内側広筋について行った。計測のための電極の位置を図5に、計測時の電極や筋電計等の配置を図6に示す。推定結果と実験値の比較は図4に示した通りである。
【産業上の利用可能性】
【００３２】
本発明の筋疲労の評価方法は、使用者の運動部位と対応する関節角度及び床反力データから、運動による使用者の筋肉疲労度をコンピュータを用いて算出するので、筋電位を実測する必要がない上装置が小型になるので産業上極めて有意義である。また、本発明のシステムは使用者の状況がリアルタイムで反映されるので、安全である上支援効率が向上するという大きな利点を有する。
【図面の簡単な説明】
【００３３】
【図１】本発明の構成を説明するための説明図である。
【図２】筋疲労度を計算するためのルーチンフローである。
【図３】回転関節のモデル図である。
【図４】本発明を用いたシミュレーションの結果と、筋電位の実測値を比較したグラフである。
【図５】実施例において筋電位を測定したときの、電極位置を示す図である。
【図６】実施例における計測機器の配置を示す図である。
【符号の説明】
【００３４】
１ 運道具
２ 疲労度を計算するコンピュータ
３ 制御部
【技術分野】
【０００１】
本発明はコンピュータを用いた筋肉疲労度の評価、及びその結果を用いた使用者の生理学的状況をリアルタイムで反映する運動支援システムに関し、特に、使用者の運動部位と対応する関節角度及び床反力データから、運動による使用者の筋肉疲労度をコンピュータを用いて算出する方法、そのための装置、それを利用した使用者の生理学的状況をリアルタイムで反映する運動支援システムに関する。
【背景技術】
【０００２】
従来、リハビリテーションは理学療法士によって、助言・指導がなされているが、理学療法士の負担が大きいために、高齢化社会、高福祉社会を迎えた今日、需要に充分答えることができないというのが実状である。そこで、従来の健康増進用機器のみならず、更に高度な技術が要求されるリハビリテーション用等の、医療機器の開発が望まれている。即ち、健康増進用機器の場合には、機器に要求される性能の任意性を広く許容することができるが、健常者とは異なる弱者に対する機器の場合には、使用者の状態に機器の方が対応できることが望まれるとともに、安全性に対する要求が一段と厳しいものとなる。
【０００３】
従来、安全性を確保する観点から、臥位で動作を行わせることによって転倒を防止したり、ハーネスなどを用いて体の自由度をある程度拘束する方法が採られており(特許文献1-4)、機器自体が大型化せざるを得ないという欠点があった。また、これまでは、人体内部の状況に対応して機器の方を調整するという発想が提案されることはなかった。しかしながら使用者の人体内部の情報を容易に推定することができ、推定した情報を機器にリアルタイムで反映することができれば、リハビリテーションなどの効果を最大にすることができる筈である。
【特許文献１】特開2000-102576号公報
【特許文献２】特開2005-192695号公報
【特許文献３】特開2005-211086号公報
【特許文献４】特開2005-211328号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
従って本発明の第一の目的は、運動時の関節角度や床反力などの人体外部の力学的諸量を用いて人体内部の状況を推定し、筋疲労度を評価する方法を提供することにある。
本発明の第二の目的は、運動時の関節角度や床反力などの人体外部の力学的諸量を用いて、人体内部の情報をデータ化すると共に、筋疲労度を評価するための装置を提供することにある。
更に本発明の第三の目的は、使用者の状況をリアルタイムで反映し、使用者の生理学的状況をリアルタイムで反映させることにより、無理なく安全であると共に効果的に運動を行うことのできる運動支援システムを提供することにある。
【課題を解決するための手段】
【０００５】
本発明者等は上記の諸目的を達成すべく鋭意研究した結果、特定の計算式を用いてコンピュータによる計算を行わせることにより、人体外部の力学的量である関節角度、及び床反力から使用者の筋疲労度を推定することができることを見出し、本発明に到着した。
即ち本発明は、使用者に装着したセンサーを介して得られた関節角度θ及び床反力Fを用いて筋張力を計算し、次いで得られた筋張力を前記センサーを用いた測定開始から測定終了までの時間範囲で積分する筋疲労度の評価方法であって、前記計算が下記一般式(1)に基づいてなされることを特徴とする筋疲労度の評価方法、その評価方法に使用する装置、及び、その評価方法を利用した使用者の生理学的状況をリアルタイムで反映する運動支援システムである。
一般式(1)
本発明においては、前記関節角度として、足関節、膝関節、股関節の各角度を採用することが好ましい。上半身が該当箇所である場合はその部位に至るまでの主要な各関節角度を入力しなければならない。該当箇所が下半身の末端に近ければ膝関節、股関節の角度は必ずしも必要ではない。さらに、人体モデルを詳細にし、入力対象となる関節数を増やすことによってより精緻な結果を得ることができる。また、安全を確保するために、設定した筋疲労度に達したときに運動装置が停止するように、安全装置を設けることが好ましい。
【発明の効果】
【０００６】
本発明によれば、人体外部の力学的量を用いて筋疲労度を推定するので、容易であるにもかかわらず使用者に負担をかけない上、得られたデータをリアルタイムで運動装置に反映させることができるので、安全に且つ最大限の効果を得ることが容易である。
【発明を実施するための最良の形態】
【０００７】
以下、図面を参照しながら本発明を詳述する。図1は本発明の実施例の概念図である。図において、符号1は足、膝、股関節角度及び床反力検出センサーを有し使用者に強制歩行させるための運動具、符合2は、得られた各関節角度および床反力に基づいて使用者の筋張力、筋疲労度を算出するためのコンピュータ、符合3は、算出された筋疲労度を運道具に反映させ、使用者の運動を制御するための制御部である。
【０００８】
本発明においては、使用者の骨格間の角度を測定し、これらの角度を、コンピュータ内にとりこまれた人間の骨格を有する人体モデルに入力して多体系問題を解き、次いで得られた関節モーメントを下記一般式(1)を含む筋モデルに適用し、最適化法を用いて使用者の筋張力決定し、最後にこの筋張力から筋疲労度を算出する（図2参照）。
一般式(1)
【０００９】
例えば、歩行時における下肢の筋肉を対象の部位とした場合には、足、膝および股関節の各角度を指定することによって、人間の骨格を有するコンピュータ内の人体モデルの姿勢を一義的に決定することができる。従って、使用者が装着する運動具に取り付けたセンサーによって使用者の各関節角度θfoot，θknee，θhipを検出し、これをコンピュータに入力すれば、コンピュータ内の人体モデルは使用者と同じ姿勢をとることになる。尚、上記のセンサーとしては、公知のセンサーの中から適宜選択して使用することができる。
【００１０】
コンピュータ内における人体モデル各部位の質量は使用者に合わせて設定され、人体モデルにかかる重力は自動的に計算されるため、人体モデルは使用者と同様の負荷を受ける。これに加え、使用者の体が外部と接触する点から受ける外力を入力すれば、人体モデルの力学的状態は使用者と同じになる。歩行の場合、運道具を含む使用者が外部と接触する点は足の裏と床の接触点のみであるから、左右の床反力FrおよびFlを入力することによって、使用者の力学的状態をコンピュータ2中に再現することができる。
【００１１】
コンピュータ2中の人体モデルは使用者の運動を反映するので、各タイムステップにおける各筋張力および筋疲労度を出力することができる。図1においては結果を出力する出力画面を有しているが、更にプリントアウトするためのプリンターを有していても良いことは当然である。各関節角度θ、床反力Fから筋張力、筋疲労度を計算するシステムフローを図2に示す。
【００１２】
(人体モデルについて)
人体各部位を適当な質量と慣性モーメントをもった剛体と仮定し、各関節を、同等の自由度と機能を持った機械的関節とみなすと、人体を一連の剛体の集まりとして数学的に表すことができ、運動時の各関節に発生する力及びトルクを得る問題は、多体系に対する逆動力学問題に帰着する。従って、各部位の質量及び慣性モーメントを非侵襲的に実測することは不可能であるが、過去の研究と統計学的アプローチによれば、身長及び体重から推定値を得ることが可能である（非特許文献１）。
【非特許文献１】阿江通良,横井孝志:日本人アスリートの身体部分慣性特性の推定 バイオメカニズム11,23-33東京大学出版会
【００１３】
また、機構学の分野においては、各関節における任意の機能を数学的に表現することがなされてきた。例えば、人体の肘や膝関節は、歩行や投擲などの大変位の運動においては、単一回転自由度のみを持つ回転関節と仮定しても充分にその運動時の挙動を推定することができることが知られている（非特許文献２）。
【非特許文献２】山崎信寿，長谷和徳:自由歩行における歩調・歩幅の生体力学的決定基準 バイオメカニズム11,179-189東京大学出版会
【００１４】
(拘束式の例 回転拘束,加速度式)
二次元平面で軸周りの回転だけを許容する回転関節は、数学的に
と表すことができる。但し、式中のｒは各剛体の並進座標、sは剛体の中心から関節の拘束中心までのベクトル、添え字iおよびjは関節によって結合される二つの剛体をそれぞれ表す。回転関節によって結合される物体iおよびjは、常にこの式を満たす範囲で運動せねばならない。
【００１５】
また、この式を時間で偏微分することにより、速度に関する拘束式（速度式）を得ることができる。速度式は回転関節によって結合される剛体が、運動中に満たさなければならない速度の方程式である。また速度式をさらに偏微分することによって、同様に加速度式を得ることができる。
【００１６】
三次元の回転拘束式は
と表される。式中のΦｓは二次元の回転拘束式と同様のものであり、剛体iおよびj上の関節中心PiおよびPjが同一点であることを意味している。またΦplはi、j双方で定義される回転の軸hiおよびhjが常に並行であることを意味している。概念図を図3に示す。具体的には、ΦｓおよびΦplのそれぞれは下記のように表される。
上式中、Aは各剛体からグローバル座標系へ変換する座標変換マトリクス、ｈは図に示される関節の回転方向を定義するベクトル、fおよびgは、P点を原点としベクトルhと直行する、関節座標系をなすベクトルである。
【００１７】
関節の拘束式に加え、時間ごとの関節角度を指定することによって、各時間における関節の状態を一義的に決定することができる。入力する関節角度の情報は、ドライバ拘束といわれるものであり、
で表される。三次元であっても、拘束式を時間微分することによって、速度式及び加速度式を得ることができる。
実際のプログラム内では、拘束式、速度式および加速度式は代数学的にライブラリ化されており、関連する物体の状態量を入力することによって計算され、時間微分を行わず直接計算することにより、高速に処理される。
【００１８】
(運動方程式)
人体の各関節を拘束式Φi(i=0,1,Λnj)(njは関節の数)で表すと、人体全体の
Newton-Euler運動方程式は
と表すことができる。式中のMは、人体の各部位の質量を対角成分にもつ対角行列、J'は各部位のローカル座標系周りの慣性モーメントを対角成分にもつ対角行列であり、「’」はローカル座標系を対象としていることを表している。ΦｒおよびΦπ'は、それぞれ並進自由度に関する拘束式およびEuler角で表したローカル座標系の回転自由度に関する拘束式である。
【００１９】
左辺の二つめの[]は変数行列であり、ｒは各部位の並進方向の加速度、ω'はローカル系での角加速度、λは運動方程式と加速度式を結びつけるLagrange乗数である。逆動力学解析においてはこのλを求めることにより、各関節に発生する力及びトルクを決定することができる。
【００２０】
右辺の行列におけるFΑは一般化外力を表す。一般化外力とは、系に作用する外力を、作用する各部位のローカル座標系で表したものである。二項目は一般化トルクを表しており、N'Αは外部から作用する外トルクである。またγは、各拘束式から得られる前述の加速度式である。
【００２１】
(解法について)
上述の運動方程式は、左辺の係数マトリクスが対称行列であり、左辺のベクトルはモデルに入力された諸条件から決定されるため、ガウス消去法、LU分解法などによって解くことできる。そして、ラグランジェ乗数ベクトルλが決定されると、各関節に生じる拘束反力および反トルクは、それぞれ
と表すことができる。式中の「”」は、これらの式が関節座標系での表現であることを表している。また、Cは関節座標系から関連する物体への座標変換マトリクス、上付きの「〜」は歪行列を表している。
【００２２】
以上の、人体モデルの構築と運動の指定および解法を用いることにより、運動中の関節反力と関節トルクの計算を行うことができる。この方法はマルチボディ・ダイナミクスに基づいており（非特許文献3）、得られた関節トルクは、複数の筋肉の筋張力によって関節周りに発生するトルクの和である。関節トルクを、適当な筋肉ごとのパラメータを持った筋モデルを用い、最適化法により解くことによって、各筋肉の筋張力を推定できることが知られている（非特許文献４）。
【非特許文献３】E.J.ハウグ:コンピュータを利用した機構解析の基本,1996,大河出版
【非特許文献４】長谷和徳,山崎信寿(1995):汎用三次元筋骨格モデルの開発 日本機械学会論文集(C編) 61巻
【００２３】
（筋モデルについて）
人体モデルから得られた関節トルクを入力したときにその関節周りの各筋が発生する筋張力を、適当な筋モデルを用い、最適化法によって推定することができる。下記一般式(1)は本発明者等が用いている疲労特性を持つ筋モデルである。
式中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋ごとのパラメータである。kは筋肉の速さ−力関係式であり、収縮要素の収縮時には
FACTOR=Min(1,3.33・q)
と表される。ここでArelおよびBrelはそれぞれ0.4および5.0前後の値をとり、Vceは収縮要素速度を収縮要素長さで割った値である。
【００２４】
収縮要素の伸張時における力と速さの関係は、その低伸張速度部においては、
と表現することができる。ここで、Fasympt（発揮張力の漸近レベル）並びにSlopefactor（速度ゼロにおける力-速さ曲線の伸張部と収縮部の傾きの比）はパラメータであり、それぞれ1.5及び2.0前後の値とすることが妥当である。
【００２５】
一方、高伸張速度部においては
と表現することができる。ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータであり、200前後の値を取る。低速度部と高速度伸張部の境界は
となる点である。以上の筋肉の力学的特性については多くの研究がなされている（非特許文献５〜７）。
【非特許文献５】Hill AV (1938): The heat of shortening andthe dynamic constants of muscle. Proceedings of the Royal Society of London,Series B 126:136-195
【非特許文献６】Felix E. Zajac(1989): Muscle and Tendon:Properties, Models, Scaling And Application to Biomechanics and Motor Control,Critical Reviews in Biomedical Engineering Volume17, Issue 4:359-411
【非特許文献７】Nagano and Gerritsen (2001): Effect ofneuromuscular strength training on vertical jumping performance-a computersimulation study. Journal of Applied Biomechanics 17:113-128
【００２６】
(疲労特性について)
前記一般式（１）中のffatigueは、本発明者等が実験によって導出した表面筋電位-疲労特性であり、
ffatigue=ｍ[(-5.0y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]＋E
と表される。式中のｍは筋肉ごとに決まるパラメータ、ｘは運動開始から現在時点までに筋が発した力の積分値であり(∫Fcedt)、yは運動の周期[秒]を表し、Eは周期が０付近での誤差を少なくするための補正項を表す。
【００２７】
上記ffatigueは、軽負荷における反復運動時に運動条件を変え、疲労に伴う積分筋電位（iEMG）の変化を観察し、得られた筋肉疲労特性のデータを回帰分析することによって決定した、表面筋電位上昇特性式である。この式は、筋張力を計算する際に、運動開始時からの積分値を引数とすると共に、前記一般式(1)式で表されるように、疲労特性を筋肉の収縮要素モデルに加えることによって疲労に伴う表面筋電位の上昇を推定するものであり、本発明者等によって始めて導出された、文献未載の式である。
【００２８】
上記した本発明の筋疲労の評価方法を実施する装置は、使用者に装着させることにより、装着した箇所の関節角度θ及び床反力を測定することのできる手段、前記θ、F及び前記一般式(1)を用いて使用者の筋疲労度を算出する手段、並びに算出された筋疲労度を出力する出力手段とからなる（図1参照）。
【００２９】
上記関節角度θ及び床反力を測定することのできる手段としては、公知のセンサーの中から適宜選択されたセンサーを使用することができ、装置を移動可能とすることもできるが、単なる足踏みとしたり、トレッドミルのような強制歩行装置を用いて非移動式とすることもできる（図6参照）。測定されたθ、Ｆ及び前記一般式(1)を用いて使用者の筋疲労度を算出する手段は、キーボードのようなθやＦを入力するための入力手段、および前記一般式(1)の演算プログラムがインストールされたコンピュータからなり、算出された筋疲労度を出力する出力手段とは、画像として表示することのできる各種ディスプレイ及び/またはプリンターである。
【００３０】
本発明の運動支援システムは、(1)使用者を強制的に運動させる、コンピュータ制御可能なトレッドミルのような運動装置、(2)使用者に装着させ、所望箇所の関節角度θ及び床反力Fを測定する測定手段、(3)得られた前記θ及びFのデータを用いて筋疲労度を算出する手段とを有する運動支援システムである。特に、前記筋疲労度のデータを前記運動装置(1)のコンピュータ制御部に入力することにより、該運動装置の強制運動に使用者の筋疲労度が反映されて制御される。このように、本発明の運動支援システムは使用者の生理学的状況をリアルタイムで反映する運動支援システムである点が最大の特徴である。これによって、本発明の運動支援システムを使用することにより、使用者に無理のない範囲で、安全且つ効率よくリハビリ等の運動を支援することができる。また、安全性を更に高める上で、使用者に許容される疲労度を予め設定し、その値に達したときに自動的に運動システムが停止するように安全装置を備えることが好ましい。
以下本発明を実施例によって更に詳述するが、本発明はこれによって限定されるものではない。
【００３１】
実施例
一般式(1)で表される筋モデル、それを含む人体モデルを用い、下記の様にして歩行時の筋疲労の推定を行った。
トレッドミル上で4km/hの速度で10分間歩行した際の表面筋電位と、その際の歩行に周期をあわせた（T≒1秒）歩行データ（各関節角度、床反力）を、一般式(1)の演算プログラムをインストールしたコンピュータに入力し、シミュレーションによって疲労度を推定した。上記表面筋電位の計測を内側広筋について行った。計測のための電極の位置を図5に、計測時の電極や筋電計等の配置を図6に示す。推定結果と実験値の比較は図4に示した通りである。
【産業上の利用可能性】
【００３２】
本発明の筋疲労の評価方法は、使用者の運動部位と対応する関節角度及び床反力データから、運動による使用者の筋肉疲労度をコンピュータを用いて算出するので、筋電位を実測する必要がない上装置が小型になるので産業上極めて有意義である。また、本発明のシステムは使用者の状況がリアルタイムで反映されるので、安全である上支援効率が向上するという大きな利点を有する。
【図面の簡単な説明】
【００３３】
【図１】本発明の構成を説明するための説明図である。
【図２】筋疲労度を計算するためのルーチンフローである。
【図３】回転関節のモデル図である。
【図４】本発明を用いたシミュレーションの結果と、筋電位の実測値を比較したグラフである。
【図５】実施例において筋電位を測定したときの、電極位置を示す図である。
【図６】実施例における計測機器の配置を示す図である。
【符号の説明】
【００３４】
１ 運道具
２ 疲労度を計算するコンピュータ
３ 制御部
【特許請求の範囲】
【請求項１】
使用者に装着したセンサーを介して得られた関節角度θ及び床反力Fを用いて筋張力を計算し、次いで得られた筋張力を前記センサーを用いた測定開始から測定終了までの時間範囲で積分する筋疲労度の評価方法であって、前記計算が下記一般式(1)に基づいてなされることを特徴とする筋疲労度の評価方法；
一般式(1)
但し、一般式(1)中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉ごとに決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式である。
【請求項２】
前記関節角度が、足関節、膝関節及び股関節の各角度である、請求項1に記載された筋疲労度の評価方法。
【請求項３】
前記kが、収縮時には、
（但し、FACTORはMin(1,3.33・q)であり、Vceは収縮要素速度を収縮要素長さで割った値である。）で表され、伸張時においては、
と表され、特に高伸張速度部においては、
（ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータである。）と表される、請求項1又は2に記載された筋疲労度の評価方法。
【請求項４】
前記ffatigueが下記の式で表される、請求項1〜3の何れかに記載された筋疲労度の評価方法；
ffatigue=ｍ[(-0.5y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]+E
但し、ｍは筋肉ごとに決まるパラメータ、xは∫Fcedtを表し、yは運動の周期（秒）を表す。Eは周期が0付近での誤差を少なくするための補正項である。
【請求項５】
使用者に装着させることにより、装着した箇所の関節角度θ及び床反力を測定することのできる手段、前記θ、F及び下記一般式(1)を用いて使用者の筋疲労度を算出する手段、並びに算出された筋疲労度を出力する出力手段とからなることを特徴とする、筋疲労度評価装置；
一般式(1)
但し、一般式(1)中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉ごとに決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式であり、それぞれ下記のように表される。
（但し、FACTORはMin(1,3.33・q)であり、Vceは収縮要素速度を収縮要素長さで割った値である。）で表され、伸張時においては、
と表され、特に高伸張速度部においては、
（ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータである。）と表され、ffatigueは下記のように表される。
ffatigue=k[(-0.5y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]+E
ここでxは∫Fcedtを表し、yは運動の周期（秒）を表す。Eは周期が0付近での誤差を少なくするための補正項である。
【請求項６】
(1)使用者を強制的に運動させる、コンピュータ制御可能な運動装置、(2)使用者に装着させ、所望箇所の関節角度θ及び床反力Fを測定する測定手段、(3)得られた前記θ及びFのデータを用いて筋疲労度を算出する手段とを有する運動支援システムであって、前記筋疲労度のデータが前記運動装置(1)のコンピュータ制御部に入力されることにより、該運動装置の強制運動が使用者の筋疲労度を反映して制御されることを特徴とする、使用者の生理学的状況をリアルタイムで反映する運動支援システム。
【請求項７】
前記筋疲労度の算出が下記一般式(1)を用いて行われる、請求項4に記載された使用者の生理学的状況をリアルタイムで反映する運動支援システム。
一般式(1)
但し、一般式(1)中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉によって決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式であり、それぞれ下記のように表される。
（但し、FACTORはMin(1,3.33・q)であり、Vceは収縮要素速度を収縮要素長さで割った値である。）で表され、伸張時においては、
と表され、特に高伸張速度部においては、
（ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータである。）と表され、ffatigueは下記のように表される；
ffatigue=k[(-0.5y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]+E
ここでxは∫Fcedtを表し、yは運動の周期（秒）を表す。Eは周期が0付近での誤差を少なくするための補正項である。
【請求項８】
設定した筋疲労度に達したときに運動装置が停止する如く安全装置が設けられた、請求項6または7に記載された使用者の生理学的状況をリアルタイムで反映する運動支援システム。
【請求項１】
使用者に装着したセンサーを介して得られた関節角度θ及び床反力Fを用いて筋張力を計算し、次いで得られた筋張力を前記センサーを用いた測定開始から測定終了までの時間範囲で積分する筋疲労度の評価方法であって、前記計算が下記一般式(1)に基づいてなされることを特徴とする筋疲労度の評価方法；
一般式(1)
但し、一般式(1)中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉ごとに決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式である。
【請求項２】
前記関節角度が、足関節、膝関節及び股関節の各角度である、請求項1に記載された筋疲労度の評価方法。
【請求項３】
前記kが、収縮時には、
（但し、FACTORはMin(1,3.33・q)であり、Vceは収縮要素速度を収縮要素長さで割った値である。）で表され、伸張時においては、
と表され、特に高伸張速度部においては、
（ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータである。）と表される、請求項1又は2に記載された筋疲労度の評価方法。
【請求項４】
前記ffatigueが下記の式で表される、請求項1〜3の何れかに記載された筋疲労度の評価方法；
ffatigue=ｍ[(-0.5y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]+E
但し、ｍは筋肉ごとに決まるパラメータ、xは∫Fcedtを表し、yは運動の周期（秒）を表す。Eは周期が0付近での誤差を少なくするための補正項である。
【請求項５】
使用者に装着させることにより、装着した箇所の関節角度θ及び床反力を測定することのできる手段、前記θ、F及び下記一般式(1)を用いて使用者の筋疲労度を算出する手段、並びに算出された筋疲労度を出力する出力手段とからなることを特徴とする、筋疲労度評価装置；
一般式(1)
但し、一般式(1)中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉ごとに決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式であり、それぞれ下記のように表される。
（但し、FACTORはMin(1,3.33・q)であり、Vceは収縮要素速度を収縮要素長さで割った値である。）で表され、伸張時においては、
と表され、特に高伸張速度部においては、
（ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータである。）と表され、ffatigueは下記のように表される。
ffatigue=k[(-0.5y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]+E
ここでxは∫Fcedtを表し、yは運動の周期（秒）を表す。Eは周期が0付近での誤差を少なくするための補正項である。
【請求項６】
(1)使用者を強制的に運動させる、コンピュータ制御可能な運動装置、(2)使用者に装着させ、所望箇所の関節角度θ及び床反力Fを測定する測定手段、(3)得られた前記θ及びFのデータを用いて筋疲労度を算出する手段とを有する運動支援システムであって、前記筋疲労度のデータが前記運動装置(1)のコンピュータ制御部に入力されることにより、該運動装置の強制運動が使用者の筋疲労度を反映して制御されることを特徴とする、使用者の生理学的状況をリアルタイムで反映する運動支援システム。
【請求項７】
前記筋疲労度の算出が下記一般式(1)を用いて行われる、請求項4に記載された使用者の生理学的状況をリアルタイムで反映する運動支援システム。
一般式(1)
但し、一般式(1)中のFceは筋肉の収縮要素力、Vceは収縮要素の収縮速度、Lceは収縮要素長さ、qは興奮水準、αは各筋肉ごとのパラメータである。また、kは筋肉によって決まる筋肉の速さ−力関係式、ffatigueは積分筋電位疲労特性式であり、それぞれ下記のように表される。
（但し、FACTORはMin(1,3.33・q)であり、Vceは収縮要素速度を収縮要素長さで割った値である。）で表され、伸張時においては、
と表され、特に高伸張速度部においては、
（ここで、Slopelinは高伸張速度部における力と速さの関係を規定するパラメータである。）と表され、ffatigueは下記のように表される；
ffatigue=k[(-0.5y+9.0)×10-8・x3+(2.0y-4.0)×10-5・x2+(-2.2y+5.9)×10-3・x+0.3]+E
ここでxは∫Fcedtを表し、yは運動の周期（秒）を表す。Eは周期が0付近での誤差を少なくするための補正項である。
【請求項８】
設定した筋疲労度に達したときに運動装置が停止する如く安全装置が設けられた、請求項6または7に記載された使用者の生理学的状況をリアルタイムで反映する運動支援システム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２００７−２３６６６３（Ｐ２００７−２３６６６３Ａ）
【公開日】平成１９年９月２０日（２００７．９．２０）
【国際特許分類】
【出願番号】特願２００６−６３７６６（Ｐ２００６−６３７６６）
【出願日】平成１８年３月９日（２００６．３．９）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り 平成１８年３月１日 社団法人精密工学会発行の「２００６年度精密工学会春季大会学術講演会講演論文集」に発表
【出願人】（５９３０２７５５４）
【Ｆターム（参考）】
【公開日】平成１９年９月２０日（２００７．９．２０）
【国際特許分類】
【出願日】平成１８年３月９日（２００６．３．９）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り 平成１８年３月１日 社団法人精密工学会発行の「２００６年度精密工学会春季大会学術講演会講演論文集」に発表
【出願人】（５９３０２７５５４）
【Ｆターム（参考）】
[ Back to top ]