JP4045487B2 - 筋骨格系機械インピーダンス計測方法及び装置 - Google Patents
筋骨格系機械インピーダンス計測方法及び装置 Download PDFInfo
- Publication number
- JP4045487B2 JP4045487B2 JP2002165150A JP2002165150A JP4045487B2 JP 4045487 B2 JP4045487 B2 JP 4045487B2 JP 2002165150 A JP2002165150 A JP 2002165150A JP 2002165150 A JP2002165150 A JP 2002165150A JP 4045487 B2 JP4045487 B2 JP 4045487B2
- Authority
- JP
- Japan
- Prior art keywords
- musculoskeletal
- mechanical impedance
- module
- external force
- uncertainty
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、筋骨格系機械インピーダンス計測方法及びそれを実現する装置に関するものであり、一方においては、医学、更には運動医学における各種の診断、訓練装置に組込まれることにより利用可能であり、他方においては、メカトロニクス機器・ロボットなどを人間的に動作させる為のシステム開発にも応用されるものであり、具体的には、計測中に含まれるノイズ及びモデル化誤差に対して、運動などにより変化する腕、足などの筋骨格系の弾性係数、粘性係数といった機械インピーダンスを定量的に計測できる筋骨格系機械インピーダンス計測方法及び装置に関する。
【０００２】
【従来の技術】
腕の機械インピーダンス計測の応用としては、ロボット制御、遠隔操作技術、仮想現実感などにおける人‐機械インタフェースの設計、筋電動力義手の応用などが考えられる。
【０００３】
【発明が解決しようとする課題】
しかしながら、実際問題として、運動中又は作業中にノイズ、モデル化誤差、その他の悪要素が不可避である為に、例えば腕の機械インピーダンスを正確に測定することは困難であった。
【０００４】
ここにおいて、本発明の解決すべき主要な目的は以下の通りである。
本発明の第１の目的は、不可避的な悪要素が存在する状況下において、筋骨格系の機械インピーダンスを計測できる、筋骨格系機械インピーダンス計測方法及び装置を提供するものである。
【０００５】
本発明の第２の目的は、計測中に含まれるノイズ、モデル化に伴う不確かさの悪要素を考慮して、筋骨格系の機械インピーダンスを計測できる、筋骨格系機械インピーダンス計測方法及び装置を提供するものである。
【０００６】
本発明の第３の目的は、多関節の筋骨格系機械インピーダンスについても計測できる筋骨格系機械インピーダンス計測方法及び装置を提供するものである。
【０００７】
本発明の他の目的は、明細書、図面、特に特許請求の範囲における各請求項の記載から自ずと明らかとなろう。
【０００８】
【課題を解決するための手段】
本発明方法は、上記課題の解決に当たり、筋骨格系に対して逐次的に摂動を与え、位置及び外力を計測した結果から、摂動に対する変化成分を抽出し、その抽出した変化成分から、ノイズ又は力学的モデル構造に依存する係数の不確かさについて取除きを行って、安定性のよい機械インピーダンスの推定を行う、という特徴的構成手法を講じる。
【０００９】
本発明装置は、上記課題の解決に当たり、筋骨格系に対して逐次的に摂動を与え、位置・外力の変化成分につき抽出を行う計測系に、当該抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行う評価モジュールと、当該抽出した変化成分と当該評価モジュールにより評価された結果に基いて、ノイズ・不確かさを取り除いて、筋骨格系の機械インピーダンスの逐次計算を行う推定モジュールとを備える、という特徴的構成手段を講じる。
【００１０】
更に、具体的詳細に述べると、当該課題の解決では、本発明が次に列挙する新規な特徴的構成手法又は手段を採用することにより、上記目的を達成するようになされる。
【００１１】
本発明方法の第１の特徴は、ノイズ又は筋骨格系のモデル化の不確かさの存在の下において、当該筋骨格系の機械インピーダンスを計測系により計測する方法であって、評価モジュール及び推定モジュールを備えた前記計測系が、位置センサを取り付けたトルクモータにより、当該トルクモータに連結され枢結部に力センサ及びグリッパを設けたリンク機構の当該グリッパに誘導外力を発生させ、前記位置センサと前記力センサとを用いて、当該グリッパの位置・外力の変化成分を抽出し、前記評価モジュールが、前記抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行うとともに、前記推定モジュールが、前記抽出を行う抽出モジュールから入力された、今回の抽出した値と、今回の関節角ベクトルに関する行列と前回求めた機械インピーダンス行列との積との差を計算し、その各成分の大きさが、前記評価モジュールにより評価された評価値に基いて、それぞれ前記ノイズ又は不確かさの前記評価値の上限の大きさを越えた場合に、その越えた程度に応じて加減することにより、筋骨格系の機械インピーダンスの逐次計算を行ってなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１２】
本発明方法の第２の特徴は、上記本発明方法の第１の特徴における前記推定モジュールが、前記筋骨格系の機械インピーダンスとして、関節毎、関節相互間の弾性係数及び粘性係数を逐次計算してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１３】
本発明方法の第３の特徴は、上記本発明方法の第１又は第２の特徴における前記評価モジュールが、前記グリッパの位置・外力の変化成分を抽出した結果に含まれる、前記筋骨格系の力学的モデルに依存する係数の不確かさ、摂動の計測誤差につき一以上評価してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１４】
本発明方法の第４の特徴は、上記本発明方法の第１の特徴における前記抽出モジュールが、今回の関節角ベクトルに関する行列の一成分として、各関節角の微分計算を行ってなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１５】
本発明方法の第５の特徴は、上記本発明方法の第１又は第４の特徴における前記推定モジュールが、前記加減すべき値を、ＵＤ分解法を用いて算出してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１６】
本発明方法の第６の特徴は、上記本発明方法の第１、第４又は第５の特徴における前記推定モジュールが、数値計算、予期し得ないノイズの累積作用の内一以上の点についても配慮して加算すべき量を算出してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１８】
本発明装置の第１の特徴は、筋骨格系に対して逐次的に摂動を与え、位置・外力の変化成分につき抽出を行う計測系に、前記抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行う評価モジュールと、前記抽出を行う抽出モジュールから入力された、今回の抽出した値と、今回の関節角ベクトルに関する行列と前回求めた機械インピーダンス行列との積との差を計算し、その各成分の大きさが、前記評価モジュールにより評価された評価値に基いて、それぞれ前記ノイズ又は不確かさの前記評価値の上限の大きさを越えた場合に、その越えた程度に応じて加減することにより、前記逐次計算をなす筋骨格系の機械インピーダンスの逐次計算を行う推定モジュールとを備えてなる、筋骨格系機械インピーダンス推定装置の構成採用にある。
【００１９】
本発明装置の第２の特徴は、上記本発明装置の第１の特徴における前記推定モジュールが、前記筋骨格系の機械インピーダンスとして、関節毎、関節相互間の弾性係数及び粘性係数を逐次計算するモジュールである、筋骨格系機械インピーダンス推定装置の構成採用にある。
【００２０】
本発明装置の第３の特徴は、上記本発明装置の第１又は第２の特徴における前記評価モジュールが、前記筋骨格系に対して与えた摂動による、位置・外力の変化成分を抽出した結果に含まれる、当該筋骨格系の力学的モデルに依存する係数の不確かさ、摂動の計測誤差につき一以上評価するモジュールである、筋骨格系機械インピーダンス推定装置の構成採用にある。
【００２２】
本発明装置の第４の特徴は、上記本発明装置の第１の特徴における前記抽出モジュールが、今回の関節角ベクトルに関する行列の一成分として、各関節角の微分計算をも行うモジュールでもある、筋骨格系機械インピーダンス計測装置の構成採用にある。
【００２３】
本発明装置の第５の特徴は、上記本発明装置の第１又は第４の特徴における前記推定モジュールが、前記加減すべき値を、ＵＤ分解法を用いて算出するモジュールでもある、筋骨格系機械インピーダンス計測装置の構成採用にある。
【００２４】
本発明装置の第６の特徴は、上記本発明装置の第１、第４又は第５の特徴における前記推定モジュールが、数値計算、予期し得ないノイズの累積作用の内一以上の点についても配慮して加算すべき量を算出するモジュールでもある、筋骨格系機械インピーダンス計測装置の構成採用にある。
【００２５】
本発明装置の第７の特徴は、上記本発明装置の第１、第２、第３、第４、第５又は第６の特徴における前記計測系が、前記筋骨格系に摂動を与える誘導外力手段と、前記筋骨格系の位置を計測する誘導位置計測手段と、前記筋骨格系に与えられた外力を計測する誘導外力計測手段と、前記筋骨格系の力学的モデルを用いて前記誘導位置計測手段、前記誘導外力測定手段によりそれぞれ計測された位置、筋発生力の摂動による変化成分を抽出する誘導変化分抽出手段とを具備してなる、節筋骨格系機械インピーダンス推定装置の構成採用にある。
【００２６】
【発明の実施の形態】
以下、本発明の実施の形態を説明するに際して、先ず、概要を説明し、その後、それを具現化した装置例及び方法例につき該当図面を参照して説明する。
【００２７】
（概要）
本発明は、ノイズやモデル化誤差を含む場合に対して、ノイズやモデル化誤差の不確かさを評価するモジュール（後述の「ノイズ・不確かさ評価モジュール」に該当）とロバストな係数逐次推定モジュールを用いることにより、関節の粘性・弾性係数をオンラインで計測する、というものである。
【００２８】
これを、既に出願している平成８特許願第２５３６２０号（特開平１０−９４５２４号公報）、平成９年特許願第６０９１９号（特開平１０‐２４８８１４号公報）の場合（以下では、関連出願として、本発明と対比する。尚、前者は、運動中、作業中での筋骨格系インピーダンスを計測できるようにしたものであり、後者は、前者がオンライン的なこと、いわゆる実計測上のことを考えなかったのに対して、実時間での逐次計算できるようにしたものである。）に適用する様に具現化した装置例及び方法例を以下順に説明する。
【００２９】
（装置例）
図１は本発明の一実施形態である筋骨格系機械インピーダンス計測装置のブロック構成図であり、図２はその装置の態様模式図である。
【００３０】
図１中、１は多関節筋骨格系、２は誘導外力手段、３は誘導位置計測手段、４は誘導外力測定手段、５は誘導変化分抽出モジュール、６はノイズ・不確かさ評価モジュール、７はロバストな係数逐次推定モジュール、８は計測系である。
【００３１】
筋骨格系機械インピーダンス計測装置は、誘導外力手段２、誘導位置計測手段３、誘導外力測定手段４、誘導変化分抽出モジュール５を具備した計測系８に、ノイズ・不確かさ評価モジュール６と、単なる係数逐次推定モジュールではない、ロバストな係数逐次推定モジュール７とを備えたものであり、各構成要素につき説明すると次の通りである。
【００３２】
誘導外力手段２は、構造的不確かな（即ち、パラメータ的な）多関節筋骨格系１に摂動を与える手段である。
【００３３】
誘導位置計測手段３は多関節筋骨格系１の位置を、誘導外力測定手段４は与えられた外力を、それぞれ計測し、誘導変化分抽出モジュール５は腕の力学モデルを用いて、位置及び筋発生力の摂動による変化成分のみを抽出する。
【００３４】
ノイズ・不確かさ評価モジュール６は予め推定した力学的構造に依存する係数の不確かさ成分及び摂動の計測誤差についての評価を行う。
【００３５】
ロバストな係数逐次推定モジュール７は、ノイズ・不確かさ評価モジュール６により行われた評価に基いて、多関節筋骨格系の機械インピーダンスを逐次的に計算する。尚、「ロバストな」とは、「ｒｏｂｕｓｔな」（堅固な：外乱に対して安定性の強い）であり、機械インピーダンスの計算を、ノイズ・不確かさに配慮しそれを取り除いて（除外して）行われるという意味である。
【００３６】
かかる誘導位置計測手段３、誘導外力測定手段４及び誘導変化分抽出モジュール５は、前述の関連出願においても構成要素として採用していたものであり、図２を参照して、誘導外力手段１と共に本発明の実施の形態を説明する上で必要な点につき説明する。
【００３７】
図２中、１１は伝導リンク機構、１１ａ乃至１１ｄはリンク、１２はグリッパ、１３は固結部、１４は力センサ、１５はテーブル、ＰＣは誘導位置計測手段３、誘導外力測定手段４、誘導変化分抽出モジュール５、ノイズ・不確かさ評価モジュール６及びロバストな係数逐次推定モジュール７の内の計算とか制御を司るコンピュータ、βは腕、β１は前腕、β２は上腕、β３は腕関節、γは肩である。
【００３８】
誘導外力手段１は、具体的態様としては、テーブル１５面上方に設けられ伝導リンク機構１１を連結構成したリンク１１ａ乃至１１ｄと、支持固定回転側の二つのリンク１１ａ、１１ｂの端の固結部１３を駆動する（図示しない）トルクモータと、先端自由可動側の二つのリンク１１ｃ、１１ｄの端の枢結部に設けられ、テーブル１５の面上方を非接触にて水平方向回転でき、テーブル１５面の法線方向の力とモーメント力とを支持するグリッパ１２からなる。
【００３９】
そして、誘導外力手段２は、コンピュータＰＣに制御され、伝導リンク機構１１を解してテーブル１５下方に配置され（図示しない）位置センサが取り付けられ（図示しない）トルクモータから手先αに誘導外力が与えられる。手先αはグリッパ１２を握持して固定され、誘導外力手段１から腕βに与えられた外力は、誘導外力測定手段４（具体的には、枢結部に設置された（図示しない）力センサ）により計測される。
【００４０】
また、手先αの位置は、誘導位置計測装置３（具体的には、トルクモータに取り付けられた（図示しない）位置センサにより計測されたモータ回転角からｘｙ座標変換を行うこと）により、求められる。腕βの節角度及び手先αに加った力を関節トルクに換算した値は、予め求めた肩関節γの位置と前腕β１・上腕β２の長さからｘｙ座標変換を行って求められる。
【００４１】
（方法例）
ここで、本発明の一実施形態である、筋骨格系機械インピーダンス計測方法について、上記装置例に適用した場合について説明する。
図３は、その筋骨格系機械インピーダンス計測方法のフロー図である。
【００４２】
ＳＴ１として、誘導外力手段１により誘導された筋骨格系１の誘導位置、誘導外力を、それぞれ誘導位置計測手段２、誘導外力測定手段３により、計測し、それに基いて、式１より、筋により発生したトルクを、トルクベクトルτｉｎとして計算する。
【００４３】
但し、ｑ［θ１ θ２］Ｔ は関節角ベクトル（θ１、θ２はそれぞれ肩関節γ、肘関節β３の角度）、Ｉは慣性行列、Ｈは遠心力、コリオリ力などの非線形干渉力ベクトルで式２、式３より算出する。またτｅｘｔは計測により得られた外部トルクベクトルである。
【００４４】
尚、Ｚ1、Ｚ２、Ｚ３は、前腕β１、上腕β２の重さ、長さ、慣性などを含んだ、腕βの構造（即ち、力学的モデル）によるパラメータ（この点からも分かる様に、このパラメータには「不確かさ」が内在している。）であり、かかる値は腕βの力学モデルを用いて予め推定しておき、式１の計算に用いられる。
【００４５】
【外１】
【００４６】
ＳＴ２として、誘導変化分抽出モジュール５により、前述のように計測（計算）した関節角ベクトルｑ及び算出した筋発生トルクベクトルτｉｎから、帯域フィルタなどを用いて、誘導による変化分ｑｆ（＝［θ1 ｆ、θ2 ｆ］Ｔ）、τｉｎ ｆ（＝［τｉｎ＿１ ｆ、τｉｎ＿２ ｆ］Ｔ）を抽出する。
【００４７】
以上のＳＴ１、ＳＴ２は、本発明を適用した前述の関連出願による。
このＳＴ２に続く手順が、本発明の要となる、ノイズ・不確かさ評価モジュール６とロバストな係数逐次推定モジュール７による機械インピーダンスの逐次算出である。
【００４８】
もっとも、このＳＴ２の結果により実際に誘動変化分抽出モジュール５により抽出されて出力された成分中には、各種の不確かな信号が含まれている。かかるノイズのような信号は、予め推定した力学的構造に依存する係数の不確かさ及び摂動の計測誤差で構成される。そして、かかる要素は、当然オンライン推定の精度に影響を及ぼす。
【００４９】
一方、そのノイズのような信号が筋発生トルクベクトルτｉｎ（［τｉｎ _ １ 、τｉｎ _ ２］Ｔ）に含まれ、多関節の機械インピーダンスとの不相関関係を求め難い。
【００５０】
そこで、ＳＴ３として、実際に含まれるノイズ、不確かさとして、予め推定した力学的構造に依存する係数（例えば、腕の長さ、重さなど）の不確かさと摂動の計測誤差とを合わせたΔを考え、ノイズ・不確かさ評価モジュール６はこのΔにつき評価した結果Δｆとして出力する手段である。即ち、その要素の上界を求める手段である。
【００５１】
更に、ＳＴ４として、既に述べような有界外乱の影響で、オンライン推定するときに、再帰的に変更される方程式の正定性を保つのが困難であり、パラメータ推定が不安定になってしまう場合がある為、多関節筋骨格系１の機械インピーダンス係数をロバストな係数逐次推定モジュール７を使って逐次に算出する。
【００５２】
以上のことを元に、更に詳細にＳＴ１乃至ＳＴ４について説明する。
【００５３】
既に述べた要領で、ＳＴ１では、誘導位置計測手段３により計測された関節角ベクトルｑが出力され、誘導外力測定手段２により計測された外力に基いて計算した筋発生トルクベクトルτｉｎが出力されて、誘導変化分抽出モジュール５に入力されることになる。
【００５４】
次にＳＴ２では、ＳＴ１により誘導変化分抽出モジュール５に入力された関節角ベクトルｑ、筋発生トルクベクトルτｉｎを基にしてその誘導変化成分をフィルタなどを用いて抽出する。また、関節角ベクトルの微分の計算（式５参照）も行い、ロバストな係数逐次推定モジュール７に、その抽出された筋発生トルクベクトルτｉｎ ｆと後述の関節角θ１、θ２に関する値を要素とする行列Ｘとが入力される。
【００５５】
即ち、筋発生トルクベクトル（信号）τｉｎが誘動変化分抽出モジュール５を通った信号τｉｎ ｆは、式４のように記述される。
τｉｎ ｆ＝ＸＵ＋Δｆ （式４）
ここで、Ｘ、Ｕは、式５、式６となる。尚、Ｒ１１、Ｄ１１はそれぞれ肩関節γの弾性と粘性、Ｒ２２、Ｄ２２はそれぞれ肘関節β３の弾性と粘性、Ｒ１２、Ｒ２１、Ｄ１２、Ｄ２１は二関節に連動した弾性と粘性の各係数であり、Δｆはノイズ・不確かさとして想定された値Δからの誘導値である。
【００５６】
【外２】
【００５７】
次に、ＳＴ３についてであるが、予め推定した力学的モデル、力学的構造に依存する係数の不確かさ（尚、摂動の計測誤差を合わせても構わない。）Δが入力され、この不確かさΔに誘導変化分抽出モジュール５でなされる抽出誤差などのノイズ・不確かさを作用させてΔｆとし、ロバストな係数逐次推定モジュール７に出力する。
【００５８】
そして、本発明の要となるＳＴ４、即ちロバストな係数逐次推定モジュール７による、ノイズ・不確かさ評価モジュール６で得られた情報（具体的には後述する、Δｆの上界ε）に基いたロバストなオンライン推定については、式７（乃至式１３）を用いて行う。
【００５９】
尚、Ｗは対角１の上三角行列、Ｚは対角行列である。Ｌｉは定数行列であり、また、正定値である。Ｖｉは正定数である。尚、ＷとＺはＰが正定値対称行列となる様に採用した分解法（ＵＤ分解法：因みに、変数としてＵは既に使用したので、ＵＤ分解法にいう、行列Ｕ、Ｄは、本明細書での説明では、行列Ｗ、Ｚに該当する。）の要素である。ｉ＝１，２（以下、同様。）である。
【００６０】
【外３】
【００６１】
【外４】
【００６２】
かかる一連の式について説明すると、▲１▼Δｆの影響を受けて、従来の推定モジュールでは不安定となるので、当該一連の式からロバストなオンライン推定とすると共に、▲２▼理論的には安定条件を満足したとしても、数値計算上及び予期し得ないノイズなどの累積作用により推定モジュールが不安定となる（式９のようにＭ行列の正定性が崩れる）場合もあり得ることから、かかる場合を回避すべく式１０の様にＵＤ分解法に基き、しかも、そのＰ（ｔ）は正定数の行列に収束しても、各成分がゼロに近い正の定数になることもあるので、実計測上で、式９の様に、Ｐ（ｔ）に正則な行列であるＬｉを加えてある。式８のＶｉについても同様の趣旨で任意に値が設定される。
【００６３】
ところで、式１０のＷｉ（ｔ）、Ｚｉ（ｔ）の決定については、Ｗｉ（ｔ−１）、Ｚｉ（ｔ−１）、設定された初期条件などから、逐次的に計算されるものであり、以下の逐次計算により決定される。尚、この逐次計算での説明において、Ｘｉは式５による。
【００６４】
時刻ｔでのＺｉの第ｊｊ成分であるＺｉ，ｊｊ（ｔ）は、ｊ＝１，２，・・，８として、式１４から求まる。
Ｚｉ，ｊｊ（ｔ）＝
（αｉ，ｊｊ（ｔ−１）／αｉ，ｊｊ（ｔ））Ｚｉ，ｊｊ（ｔ−１）（式１４）
但し、αｉ，ｊ（ｔ）＝αｉ，ｊ（ｔ−１）＋ｆｉ，ｊｂｉ，ｊ
αｉ，ｊ（０）は初期条件として経験則で定めた正の定数値である。
ｆｉ，ｊ，ｂｉ，ｊは、次式のｆｉ，ｂｉの第ｊ要素である。
ｆｉ＝Ｗｉ Ｔ（ｔ−１）Ｘｉ Ｔ（ｔ） ｂｉ＝Ｚｉ（ｔ−１）ｆｉ
【００６５】
時刻ｔでのＷｉの第ｊｍ成分Ｗｉ，ｊｍ（ｔ）は、ｊについてｊ＝１，・・，ｍ−１で、ｍにつきｍ＝１，２，・・，８として、式１５から求まる。
Ｗｉ，ｊｍ（ｔ）＝Ｗｉ，ｊｍ（ｔ−１）＋ｖｉ，ｊｋｉ，ｍ （式１５）
但し、ｔが更新する毎に ｖｉ，ｊ←ｖｉ，ｊ＋Ｗｉ，ｊｍｖｉ，ｍとする。
尚、ｖｉ，ｊ＝ｂｉ，ｊ 、 ｋｉ，ｍ＝‐ｆｉ，ｍ／αｉ，ｍ（ｔ−１）
【００６６】
【外５】
【００６７】
【実施例】
【外６】
【００６８】
図４は、腕βを動かし、更に腕βの機械インピーダンス係数を変化させた場合における、従来の手法（関連出願）によるシミュレーション結果であり、（ａ）は腕βの肩・肘関節の角度θ１、θ２、（ｂ）は算出した筋発生トルクτ１、τ２、（ｃ）と（ｄ）は係数逐次推定モジュールにより算出した弾性係数Ｒ１１、Ｒ１２とＲ２１、Ｒ２２、（ｅ）と（ｆ）は、同様に算出した粘性係数Ｄ１１、Ｄ１２とＤ２１、Ｄ２２の時間変化である。尚、（ｃ）（ｄ）のグラフにおいて、実線は弾性係数の従来法での推定値、点線は弾性係数の真値（設定値）を示し、同様、（ｅ）（ｆ）のグラフにおいて、実線は粘性係数の従来法での推定値、点線は粘性係数の真値（設定値）を示す。
【００６９】
図５は、図４と同じ条件の下、腕βを動かし、更に腕βの機械インピーダンス係数を変化させた場合における、本発明でのシミュレーション結果であって、（ａ）は腕βの肩・肘関節の角度θ１、θ２、（ｂ）は算出した筋発生トルクτ１、τ２、（ｃ）と（ｄ）はノイズ・不確かさ評価モジュール６及びロバストな係数逐次推定モジュール７により算出した弾性係数Ｒ１１、Ｒ１２とＲ２１、Ｒ２２、（ｅ）と（ｆ）は、同様に算出した粘性係数Ｄ１１、Ｄ１２とＤ２１、Ｄ２２の時間変化である。尚、（ｃ）（ｄ）のグラフにおいて、実線は弾性係数の本発明による推定値、点線は弾性係数の真値（設定値）を示し、同様、（ｅ）（ｆ）のグラフにおいて、実線は粘性係数の本発明による推定値、点線は粘性係数の真値（設定値）を示す。
【００７０】
両図から一瞥して、従来ではノイズの影響で時間と共に変化する各弾性と粘性係数を推定できなかったところ、本発明により良好な結果が得られることが、分かる。
【００７１】
以上、本発明の実施の形態である装置例及び方法例を説明し、更に実施例にてシミュレーションの結果を示したものの、本発明の目的を達し、下記する効果を奏する範囲において、適宜変更して実施可能である。
【００７２】
【発明の効果】
本発明によれば、実時計測で起こるノイズ及びモデル化誤差の存在の下において、運動、作業などにより変化する腕、足などの筋骨格系の粘性弾性を定量的に計測できる。これにより、医学、更には運動医学における各種の診断、訓練装置にも組込まれ、また、メカトロニクス機器・ロボットなどを人間的に動作させる為のシステム開発にも応用されるなどの極めて優れた効果を発揮する。
【図面の簡単な説明】
【図１】本発明の一実施形態である筋骨格系機械インピーダンス計測装置のブロック構成図である。
【図２】図１の筋骨格系機械インピーダンス計測装置の態様模式図である
【図３】本発明の一実施形態である筋骨格系機械インピーダンス計測方法のフロー図である。
【図４】従来でのシミュレーション結果で、（ａ）は関節角、（ｂ）は関節トルク、（ｃ）（ｄ）は弾性係数、（ｅ）（ｆ）は粘性係数についての時間に対する結果である。
【図５】本発明に基くシミュレーション結果で、（ａ）は関節角、（ｂ）は関節トルク、（ｃ）（ｄ）は弾性係数、（ｅ）（ｆ）は粘性係数についての時間に対する結果である。
【符号の説明】
１…多関節筋骨格系
１１…伝導リンク機構
１１ａ、１１ｂ、１１ｃ、１１ｄ…リンク
１２…グリッパ
１３…固結部
１４…力センサ
１５…テーブル
２…誘導外力手段
３…誘導位置計測手段
４…誘導外力測定手段
５…誘導変化分抽出モジュール
６…ノイズ・不確かさ評価モジュール
７…ロバストな係数逐次推定モジュール
８…計測系
ＰＣ…コンピュータ
α…手先
β…腕
β１…前腕
β２…上腕
β３…肘関節
γ…肩関節
【発明の属する技術分野】
本発明は、筋骨格系機械インピーダンス計測方法及びそれを実現する装置に関するものであり、一方においては、医学、更には運動医学における各種の診断、訓練装置に組込まれることにより利用可能であり、他方においては、メカトロニクス機器・ロボットなどを人間的に動作させる為のシステム開発にも応用されるものであり、具体的には、計測中に含まれるノイズ及びモデル化誤差に対して、運動などにより変化する腕、足などの筋骨格系の弾性係数、粘性係数といった機械インピーダンスを定量的に計測できる筋骨格系機械インピーダンス計測方法及び装置に関する。
【０００２】
【従来の技術】
腕の機械インピーダンス計測の応用としては、ロボット制御、遠隔操作技術、仮想現実感などにおける人‐機械インタフェースの設計、筋電動力義手の応用などが考えられる。
【０００３】
【発明が解決しようとする課題】
しかしながら、実際問題として、運動中又は作業中にノイズ、モデル化誤差、その他の悪要素が不可避である為に、例えば腕の機械インピーダンスを正確に測定することは困難であった。
【０００４】
ここにおいて、本発明の解決すべき主要な目的は以下の通りである。
本発明の第１の目的は、不可避的な悪要素が存在する状況下において、筋骨格系の機械インピーダンスを計測できる、筋骨格系機械インピーダンス計測方法及び装置を提供するものである。
【０００５】
本発明の第２の目的は、計測中に含まれるノイズ、モデル化に伴う不確かさの悪要素を考慮して、筋骨格系の機械インピーダンスを計測できる、筋骨格系機械インピーダンス計測方法及び装置を提供するものである。
【０００６】
本発明の第３の目的は、多関節の筋骨格系機械インピーダンスについても計測できる筋骨格系機械インピーダンス計測方法及び装置を提供するものである。
【０００７】
本発明の他の目的は、明細書、図面、特に特許請求の範囲における各請求項の記載から自ずと明らかとなろう。
【０００８】
【課題を解決するための手段】
本発明方法は、上記課題の解決に当たり、筋骨格系に対して逐次的に摂動を与え、位置及び外力を計測した結果から、摂動に対する変化成分を抽出し、その抽出した変化成分から、ノイズ又は力学的モデル構造に依存する係数の不確かさについて取除きを行って、安定性のよい機械インピーダンスの推定を行う、という特徴的構成手法を講じる。
【０００９】
本発明装置は、上記課題の解決に当たり、筋骨格系に対して逐次的に摂動を与え、位置・外力の変化成分につき抽出を行う計測系に、当該抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行う評価モジュールと、当該抽出した変化成分と当該評価モジュールにより評価された結果に基いて、ノイズ・不確かさを取り除いて、筋骨格系の機械インピーダンスの逐次計算を行う推定モジュールとを備える、という特徴的構成手段を講じる。
【００１０】
更に、具体的詳細に述べると、当該課題の解決では、本発明が次に列挙する新規な特徴的構成手法又は手段を採用することにより、上記目的を達成するようになされる。
【００１１】
本発明方法の第１の特徴は、ノイズ又は筋骨格系のモデル化の不確かさの存在の下において、当該筋骨格系の機械インピーダンスを計測系により計測する方法であって、評価モジュール及び推定モジュールを備えた前記計測系が、位置センサを取り付けたトルクモータにより、当該トルクモータに連結され枢結部に力センサ及びグリッパを設けたリンク機構の当該グリッパに誘導外力を発生させ、前記位置センサと前記力センサとを用いて、当該グリッパの位置・外力の変化成分を抽出し、前記評価モジュールが、前記抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行うとともに、前記推定モジュールが、前記抽出を行う抽出モジュールから入力された、今回の抽出した値と、今回の関節角ベクトルに関する行列と前回求めた機械インピーダンス行列との積との差を計算し、その各成分の大きさが、前記評価モジュールにより評価された評価値に基いて、それぞれ前記ノイズ又は不確かさの前記評価値の上限の大きさを越えた場合に、その越えた程度に応じて加減することにより、筋骨格系の機械インピーダンスの逐次計算を行ってなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１２】
本発明方法の第２の特徴は、上記本発明方法の第１の特徴における前記推定モジュールが、前記筋骨格系の機械インピーダンスとして、関節毎、関節相互間の弾性係数及び粘性係数を逐次計算してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１３】
本発明方法の第３の特徴は、上記本発明方法の第１又は第２の特徴における前記評価モジュールが、前記グリッパの位置・外力の変化成分を抽出した結果に含まれる、前記筋骨格系の力学的モデルに依存する係数の不確かさ、摂動の計測誤差につき一以上評価してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１４】
本発明方法の第４の特徴は、上記本発明方法の第１の特徴における前記抽出モジュールが、今回の関節角ベクトルに関する行列の一成分として、各関節角の微分計算を行ってなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１５】
本発明方法の第５の特徴は、上記本発明方法の第１又は第４の特徴における前記推定モジュールが、前記加減すべき値を、ＵＤ分解法を用いて算出してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１６】
本発明方法の第６の特徴は、上記本発明方法の第１、第４又は第５の特徴における前記推定モジュールが、数値計算、予期し得ないノイズの累積作用の内一以上の点についても配慮して加算すべき量を算出してなる、筋骨格系機械インピーダンス計測方法の構成採用にある。
【００１８】
本発明装置の第１の特徴は、筋骨格系に対して逐次的に摂動を与え、位置・外力の変化成分につき抽出を行う計測系に、前記抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行う評価モジュールと、前記抽出を行う抽出モジュールから入力された、今回の抽出した値と、今回の関節角ベクトルに関する行列と前回求めた機械インピーダンス行列との積との差を計算し、その各成分の大きさが、前記評価モジュールにより評価された評価値に基いて、それぞれ前記ノイズ又は不確かさの前記評価値の上限の大きさを越えた場合に、その越えた程度に応じて加減することにより、前記逐次計算をなす筋骨格系の機械インピーダンスの逐次計算を行う推定モジュールとを備えてなる、筋骨格系機械インピーダンス推定装置の構成採用にある。
【００１９】
本発明装置の第２の特徴は、上記本発明装置の第１の特徴における前記推定モジュールが、前記筋骨格系の機械インピーダンスとして、関節毎、関節相互間の弾性係数及び粘性係数を逐次計算するモジュールである、筋骨格系機械インピーダンス推定装置の構成採用にある。
【００２０】
本発明装置の第３の特徴は、上記本発明装置の第１又は第２の特徴における前記評価モジュールが、前記筋骨格系に対して与えた摂動による、位置・外力の変化成分を抽出した結果に含まれる、当該筋骨格系の力学的モデルに依存する係数の不確かさ、摂動の計測誤差につき一以上評価するモジュールである、筋骨格系機械インピーダンス推定装置の構成採用にある。
【００２２】
本発明装置の第４の特徴は、上記本発明装置の第１の特徴における前記抽出モジュールが、今回の関節角ベクトルに関する行列の一成分として、各関節角の微分計算をも行うモジュールでもある、筋骨格系機械インピーダンス計測装置の構成採用にある。
【００２３】
本発明装置の第５の特徴は、上記本発明装置の第１又は第４の特徴における前記推定モジュールが、前記加減すべき値を、ＵＤ分解法を用いて算出するモジュールでもある、筋骨格系機械インピーダンス計測装置の構成採用にある。
【００２４】
本発明装置の第６の特徴は、上記本発明装置の第１、第４又は第５の特徴における前記推定モジュールが、数値計算、予期し得ないノイズの累積作用の内一以上の点についても配慮して加算すべき量を算出するモジュールでもある、筋骨格系機械インピーダンス計測装置の構成採用にある。
【００２５】
本発明装置の第７の特徴は、上記本発明装置の第１、第２、第３、第４、第５又は第６の特徴における前記計測系が、前記筋骨格系に摂動を与える誘導外力手段と、前記筋骨格系の位置を計測する誘導位置計測手段と、前記筋骨格系に与えられた外力を計測する誘導外力計測手段と、前記筋骨格系の力学的モデルを用いて前記誘導位置計測手段、前記誘導外力測定手段によりそれぞれ計測された位置、筋発生力の摂動による変化成分を抽出する誘導変化分抽出手段とを具備してなる、節筋骨格系機械インピーダンス推定装置の構成採用にある。
【００２６】
【発明の実施の形態】
以下、本発明の実施の形態を説明するに際して、先ず、概要を説明し、その後、それを具現化した装置例及び方法例につき該当図面を参照して説明する。
【００２７】
（概要）
本発明は、ノイズやモデル化誤差を含む場合に対して、ノイズやモデル化誤差の不確かさを評価するモジュール（後述の「ノイズ・不確かさ評価モジュール」に該当）とロバストな係数逐次推定モジュールを用いることにより、関節の粘性・弾性係数をオンラインで計測する、というものである。
【００２８】
これを、既に出願している平成８特許願第２５３６２０号（特開平１０−９４５２４号公報）、平成９年特許願第６０９１９号（特開平１０‐２４８８１４号公報）の場合（以下では、関連出願として、本発明と対比する。尚、前者は、運動中、作業中での筋骨格系インピーダンスを計測できるようにしたものであり、後者は、前者がオンライン的なこと、いわゆる実計測上のことを考えなかったのに対して、実時間での逐次計算できるようにしたものである。）に適用する様に具現化した装置例及び方法例を以下順に説明する。
【００２９】
（装置例）
図１は本発明の一実施形態である筋骨格系機械インピーダンス計測装置のブロック構成図であり、図２はその装置の態様模式図である。
【００３０】
図１中、１は多関節筋骨格系、２は誘導外力手段、３は誘導位置計測手段、４は誘導外力測定手段、５は誘導変化分抽出モジュール、６はノイズ・不確かさ評価モジュール、７はロバストな係数逐次推定モジュール、８は計測系である。
【００３１】
筋骨格系機械インピーダンス計測装置は、誘導外力手段２、誘導位置計測手段３、誘導外力測定手段４、誘導変化分抽出モジュール５を具備した計測系８に、ノイズ・不確かさ評価モジュール６と、単なる係数逐次推定モジュールではない、ロバストな係数逐次推定モジュール７とを備えたものであり、各構成要素につき説明すると次の通りである。
【００３２】
誘導外力手段２は、構造的不確かな（即ち、パラメータ的な）多関節筋骨格系１に摂動を与える手段である。
【００３３】
誘導位置計測手段３は多関節筋骨格系１の位置を、誘導外力測定手段４は与えられた外力を、それぞれ計測し、誘導変化分抽出モジュール５は腕の力学モデルを用いて、位置及び筋発生力の摂動による変化成分のみを抽出する。
【００３４】
ノイズ・不確かさ評価モジュール６は予め推定した力学的構造に依存する係数の不確かさ成分及び摂動の計測誤差についての評価を行う。
【００３５】
ロバストな係数逐次推定モジュール７は、ノイズ・不確かさ評価モジュール６により行われた評価に基いて、多関節筋骨格系の機械インピーダンスを逐次的に計算する。尚、「ロバストな」とは、「ｒｏｂｕｓｔな」（堅固な：外乱に対して安定性の強い）であり、機械インピーダンスの計算を、ノイズ・不確かさに配慮しそれを取り除いて（除外して）行われるという意味である。
【００３６】
かかる誘導位置計測手段３、誘導外力測定手段４及び誘導変化分抽出モジュール５は、前述の関連出願においても構成要素として採用していたものであり、図２を参照して、誘導外力手段１と共に本発明の実施の形態を説明する上で必要な点につき説明する。
【００３７】
図２中、１１は伝導リンク機構、１１ａ乃至１１ｄはリンク、１２はグリッパ、１３は固結部、１４は力センサ、１５はテーブル、ＰＣは誘導位置計測手段３、誘導外力測定手段４、誘導変化分抽出モジュール５、ノイズ・不確かさ評価モジュール６及びロバストな係数逐次推定モジュール７の内の計算とか制御を司るコンピュータ、βは腕、β１は前腕、β２は上腕、β３は腕関節、γは肩である。
【００３８】
誘導外力手段１は、具体的態様としては、テーブル１５面上方に設けられ伝導リンク機構１１を連結構成したリンク１１ａ乃至１１ｄと、支持固定回転側の二つのリンク１１ａ、１１ｂの端の固結部１３を駆動する（図示しない）トルクモータと、先端自由可動側の二つのリンク１１ｃ、１１ｄの端の枢結部に設けられ、テーブル１５の面上方を非接触にて水平方向回転でき、テーブル１５面の法線方向の力とモーメント力とを支持するグリッパ１２からなる。
【００３９】
そして、誘導外力手段２は、コンピュータＰＣに制御され、伝導リンク機構１１を解してテーブル１５下方に配置され（図示しない）位置センサが取り付けられ（図示しない）トルクモータから手先αに誘導外力が与えられる。手先αはグリッパ１２を握持して固定され、誘導外力手段１から腕βに与えられた外力は、誘導外力測定手段４（具体的には、枢結部に設置された（図示しない）力センサ）により計測される。
【００４０】
また、手先αの位置は、誘導位置計測装置３（具体的には、トルクモータに取り付けられた（図示しない）位置センサにより計測されたモータ回転角からｘｙ座標変換を行うこと）により、求められる。腕βの節角度及び手先αに加った力を関節トルクに換算した値は、予め求めた肩関節γの位置と前腕β１・上腕β２の長さからｘｙ座標変換を行って求められる。
【００４１】
（方法例）
ここで、本発明の一実施形態である、筋骨格系機械インピーダンス計測方法について、上記装置例に適用した場合について説明する。
図３は、その筋骨格系機械インピーダンス計測方法のフロー図である。
【００４２】
ＳＴ１として、誘導外力手段１により誘導された筋骨格系１の誘導位置、誘導外力を、それぞれ誘導位置計測手段２、誘導外力測定手段３により、計測し、それに基いて、式１より、筋により発生したトルクを、トルクベクトルτｉｎとして計算する。
【００４３】
但し、ｑ［θ１ θ２］Ｔ は関節角ベクトル（θ１、θ２はそれぞれ肩関節γ、肘関節β３の角度）、Ｉは慣性行列、Ｈは遠心力、コリオリ力などの非線形干渉力ベクトルで式２、式３より算出する。またτｅｘｔは計測により得られた外部トルクベクトルである。
【００４４】
尚、Ｚ1、Ｚ２、Ｚ３は、前腕β１、上腕β２の重さ、長さ、慣性などを含んだ、腕βの構造（即ち、力学的モデル）によるパラメータ（この点からも分かる様に、このパラメータには「不確かさ」が内在している。）であり、かかる値は腕βの力学モデルを用いて予め推定しておき、式１の計算に用いられる。
【００４５】
【外１】
ＳＴ２として、誘導変化分抽出モジュール５により、前述のように計測（計算）した関節角ベクトルｑ及び算出した筋発生トルクベクトルτｉｎから、帯域フィルタなどを用いて、誘導による変化分ｑｆ（＝［θ1 ｆ、θ2 ｆ］Ｔ）、τｉｎ ｆ（＝［τｉｎ＿１ ｆ、τｉｎ＿２ ｆ］Ｔ）を抽出する。
【００４７】
以上のＳＴ１、ＳＴ２は、本発明を適用した前述の関連出願による。
このＳＴ２に続く手順が、本発明の要となる、ノイズ・不確かさ評価モジュール６とロバストな係数逐次推定モジュール７による機械インピーダンスの逐次算出である。
【００４８】
もっとも、このＳＴ２の結果により実際に誘動変化分抽出モジュール５により抽出されて出力された成分中には、各種の不確かな信号が含まれている。かかるノイズのような信号は、予め推定した力学的構造に依存する係数の不確かさ及び摂動の計測誤差で構成される。そして、かかる要素は、当然オンライン推定の精度に影響を及ぼす。
【００４９】
一方、そのノイズのような信号が筋発生トルクベクトルτｉｎ（［τｉｎ _ １ 、τｉｎ _ ２］Ｔ）に含まれ、多関節の機械インピーダンスとの不相関関係を求め難い。
【００５０】
そこで、ＳＴ３として、実際に含まれるノイズ、不確かさとして、予め推定した力学的構造に依存する係数（例えば、腕の長さ、重さなど）の不確かさと摂動の計測誤差とを合わせたΔを考え、ノイズ・不確かさ評価モジュール６はこのΔにつき評価した結果Δｆとして出力する手段である。即ち、その要素の上界を求める手段である。
【００５１】
更に、ＳＴ４として、既に述べような有界外乱の影響で、オンライン推定するときに、再帰的に変更される方程式の正定性を保つのが困難であり、パラメータ推定が不安定になってしまう場合がある為、多関節筋骨格系１の機械インピーダンス係数をロバストな係数逐次推定モジュール７を使って逐次に算出する。
【００５２】
以上のことを元に、更に詳細にＳＴ１乃至ＳＴ４について説明する。
【００５３】
既に述べた要領で、ＳＴ１では、誘導位置計測手段３により計測された関節角ベクトルｑが出力され、誘導外力測定手段２により計測された外力に基いて計算した筋発生トルクベクトルτｉｎが出力されて、誘導変化分抽出モジュール５に入力されることになる。
【００５４】
次にＳＴ２では、ＳＴ１により誘導変化分抽出モジュール５に入力された関節角ベクトルｑ、筋発生トルクベクトルτｉｎを基にしてその誘導変化成分をフィルタなどを用いて抽出する。また、関節角ベクトルの微分の計算（式５参照）も行い、ロバストな係数逐次推定モジュール７に、その抽出された筋発生トルクベクトルτｉｎ ｆと後述の関節角θ１、θ２に関する値を要素とする行列Ｘとが入力される。
【００５５】
即ち、筋発生トルクベクトル（信号）τｉｎが誘動変化分抽出モジュール５を通った信号τｉｎ ｆは、式４のように記述される。
τｉｎ ｆ＝ＸＵ＋Δｆ （式４）
ここで、Ｘ、Ｕは、式５、式６となる。尚、Ｒ１１、Ｄ１１はそれぞれ肩関節γの弾性と粘性、Ｒ２２、Ｄ２２はそれぞれ肘関節β３の弾性と粘性、Ｒ１２、Ｒ２１、Ｄ１２、Ｄ２１は二関節に連動した弾性と粘性の各係数であり、Δｆはノイズ・不確かさとして想定された値Δからの誘導値である。
【００５６】
【外２】
次に、ＳＴ３についてであるが、予め推定した力学的モデル、力学的構造に依存する係数の不確かさ（尚、摂動の計測誤差を合わせても構わない。）Δが入力され、この不確かさΔに誘導変化分抽出モジュール５でなされる抽出誤差などのノイズ・不確かさを作用させてΔｆとし、ロバストな係数逐次推定モジュール７に出力する。
【００５８】
そして、本発明の要となるＳＴ４、即ちロバストな係数逐次推定モジュール７による、ノイズ・不確かさ評価モジュール６で得られた情報（具体的には後述する、Δｆの上界ε）に基いたロバストなオンライン推定については、式７（乃至式１３）を用いて行う。
【００５９】
尚、Ｗは対角１の上三角行列、Ｚは対角行列である。Ｌｉは定数行列であり、また、正定値である。Ｖｉは正定数である。尚、ＷとＺはＰが正定値対称行列となる様に採用した分解法（ＵＤ分解法：因みに、変数としてＵは既に使用したので、ＵＤ分解法にいう、行列Ｕ、Ｄは、本明細書での説明では、行列Ｗ、Ｚに該当する。）の要素である。ｉ＝１，２（以下、同様。）である。
【００６０】
【外３】
【外４】
かかる一連の式について説明すると、▲１▼Δｆの影響を受けて、従来の推定モジュールでは不安定となるので、当該一連の式からロバストなオンライン推定とすると共に、▲２▼理論的には安定条件を満足したとしても、数値計算上及び予期し得ないノイズなどの累積作用により推定モジュールが不安定となる（式９のようにＭ行列の正定性が崩れる）場合もあり得ることから、かかる場合を回避すべく式１０の様にＵＤ分解法に基き、しかも、そのＰ（ｔ）は正定数の行列に収束しても、各成分がゼロに近い正の定数になることもあるので、実計測上で、式９の様に、Ｐ（ｔ）に正則な行列であるＬｉを加えてある。式８のＶｉについても同様の趣旨で任意に値が設定される。
【００６３】
ところで、式１０のＷｉ（ｔ）、Ｚｉ（ｔ）の決定については、Ｗｉ（ｔ−１）、Ｚｉ（ｔ−１）、設定された初期条件などから、逐次的に計算されるものであり、以下の逐次計算により決定される。尚、この逐次計算での説明において、Ｘｉは式５による。
【００６４】
時刻ｔでのＺｉの第ｊｊ成分であるＺｉ，ｊｊ（ｔ）は、ｊ＝１，２，・・，８として、式１４から求まる。
Ｚｉ，ｊｊ（ｔ）＝
（αｉ，ｊｊ（ｔ−１）／αｉ，ｊｊ（ｔ））Ｚｉ，ｊｊ（ｔ−１）（式１４）
但し、αｉ，ｊ（ｔ）＝αｉ，ｊ（ｔ−１）＋ｆｉ，ｊｂｉ，ｊ
αｉ，ｊ（０）は初期条件として経験則で定めた正の定数値である。
ｆｉ，ｊ，ｂｉ，ｊは、次式のｆｉ，ｂｉの第ｊ要素である。
ｆｉ＝Ｗｉ Ｔ（ｔ−１）Ｘｉ Ｔ（ｔ） ｂｉ＝Ｚｉ（ｔ−１）ｆｉ
【００６５】
時刻ｔでのＷｉの第ｊｍ成分Ｗｉ，ｊｍ（ｔ）は、ｊについてｊ＝１，・・，ｍ−１で、ｍにつきｍ＝１，２，・・，８として、式１５から求まる。
Ｗｉ，ｊｍ（ｔ）＝Ｗｉ，ｊｍ（ｔ−１）＋ｖｉ，ｊｋｉ，ｍ （式１５）
但し、ｔが更新する毎に ｖｉ，ｊ←ｖｉ，ｊ＋Ｗｉ，ｊｍｖｉ，ｍとする。
尚、ｖｉ，ｊ＝ｂｉ，ｊ 、 ｋｉ，ｍ＝‐ｆｉ，ｍ／αｉ，ｍ（ｔ−１）
【００６６】
【外５】
【実施例】
【外６】
図４は、腕βを動かし、更に腕βの機械インピーダンス係数を変化させた場合における、従来の手法（関連出願）によるシミュレーション結果であり、（ａ）は腕βの肩・肘関節の角度θ１、θ２、（ｂ）は算出した筋発生トルクτ１、τ２、（ｃ）と（ｄ）は係数逐次推定モジュールにより算出した弾性係数Ｒ１１、Ｒ１２とＲ２１、Ｒ２２、（ｅ）と（ｆ）は、同様に算出した粘性係数Ｄ１１、Ｄ１２とＤ２１、Ｄ２２の時間変化である。尚、（ｃ）（ｄ）のグラフにおいて、実線は弾性係数の従来法での推定値、点線は弾性係数の真値（設定値）を示し、同様、（ｅ）（ｆ）のグラフにおいて、実線は粘性係数の従来法での推定値、点線は粘性係数の真値（設定値）を示す。
【００６９】
図５は、図４と同じ条件の下、腕βを動かし、更に腕βの機械インピーダンス係数を変化させた場合における、本発明でのシミュレーション結果であって、（ａ）は腕βの肩・肘関節の角度θ１、θ２、（ｂ）は算出した筋発生トルクτ１、τ２、（ｃ）と（ｄ）はノイズ・不確かさ評価モジュール６及びロバストな係数逐次推定モジュール７により算出した弾性係数Ｒ１１、Ｒ１２とＲ２１、Ｒ２２、（ｅ）と（ｆ）は、同様に算出した粘性係数Ｄ１１、Ｄ１２とＤ２１、Ｄ２２の時間変化である。尚、（ｃ）（ｄ）のグラフにおいて、実線は弾性係数の本発明による推定値、点線は弾性係数の真値（設定値）を示し、同様、（ｅ）（ｆ）のグラフにおいて、実線は粘性係数の本発明による推定値、点線は粘性係数の真値（設定値）を示す。
【００７０】
両図から一瞥して、従来ではノイズの影響で時間と共に変化する各弾性と粘性係数を推定できなかったところ、本発明により良好な結果が得られることが、分かる。
【００７１】
以上、本発明の実施の形態である装置例及び方法例を説明し、更に実施例にてシミュレーションの結果を示したものの、本発明の目的を達し、下記する効果を奏する範囲において、適宜変更して実施可能である。
【００７２】
【発明の効果】
本発明によれば、実時計測で起こるノイズ及びモデル化誤差の存在の下において、運動、作業などにより変化する腕、足などの筋骨格系の粘性弾性を定量的に計測できる。これにより、医学、更には運動医学における各種の診断、訓練装置にも組込まれ、また、メカトロニクス機器・ロボットなどを人間的に動作させる為のシステム開発にも応用されるなどの極めて優れた効果を発揮する。
【図面の簡単な説明】
【図１】本発明の一実施形態である筋骨格系機械インピーダンス計測装置のブロック構成図である。
【図２】図１の筋骨格系機械インピーダンス計測装置の態様模式図である
【図３】本発明の一実施形態である筋骨格系機械インピーダンス計測方法のフロー図である。
【図４】従来でのシミュレーション結果で、（ａ）は関節角、（ｂ）は関節トルク、（ｃ）（ｄ）は弾性係数、（ｅ）（ｆ）は粘性係数についての時間に対する結果である。
【図５】本発明に基くシミュレーション結果で、（ａ）は関節角、（ｂ）は関節トルク、（ｃ）（ｄ）は弾性係数、（ｅ）（ｆ）は粘性係数についての時間に対する結果である。
【符号の説明】
１…多関節筋骨格系
１１…伝導リンク機構
１１ａ、１１ｂ、１１ｃ、１１ｄ…リンク
１２…グリッパ
１３…固結部
１４…力センサ
１５…テーブル
２…誘導外力手段
３…誘導位置計測手段
４…誘導外力測定手段
５…誘導変化分抽出モジュール
６…ノイズ・不確かさ評価モジュール
７…ロバストな係数逐次推定モジュール
８…計測系
ＰＣ…コンピュータ
α…手先
β…腕
β１…前腕
β２…上腕
β３…肘関節
γ…肩関節
Claims (13)
- ノイズ又は筋骨格系のモデル化の不確かさの存在の下において、当該筋骨格系の機械インピーダンスを計測系により計測する方法であって、
評価モジュール及び推定モジュールを備えた前記計測系が、位置センサを取り付けたトルクモータにより、当該トルクモータに連結され枢結部に力センサ及びグリッパを設けたリンク機構の当該グリッパに誘導外力を発生させ、前記位置センサと前記力センサとを用いて、当該グリッパの位置・外力の変化成分を抽出し、
前記評価モジュールが、前記抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行うとともに、
前記推定モジュールが、前記抽出を行う抽出モジュールから入力された、今回の抽出した値と、今回の関節角ベクトルに関する行列と前回求めた機械インピーダンス行列との積との差を計算し、その各成分の大きさが、前記評価モジュールにより評価された評価値に基いて、それぞれ前記ノイズ又は不確かさの前記評価値の上限の大きさを越えた場合に、その越えた程度に応じて加減することにより、筋骨格系の機械インピーダンスの逐次計算を行う、
ことを特徴とする筋骨格系機械インピーダンス計測方法。 - 前記推定モジュールは、
前記筋骨格系の機械インピーダンスとして、関節毎、関節相互間の弾性係数及び粘性係数を逐次計算する、
ことを特徴とする請求項１に記載の筋骨格系機械インピーダンス計測方法。 - 前記評価モジュールは、
前記グリッパの位置・外力の変化成分を抽出した結果に含まれる、前記筋骨格系の力学的モデルに依存する係数の不確かさ、摂動の計測誤差につき一以上評価する、
ことを特徴とする請求項１又は２に記載の筋骨格系機械インピーダンス計測方法。 - 前記抽出モジュールは、
今回の関節角ベクトルに関する行列の一成分として、各関節角の微分計算を行う、
ことを特徴とする請求項１に記載の筋骨格系機械インピーダンス計測方法。 - 前記推定モジュールは、
前記加減すべき値を、ＵＤ分解法を用いて算出する、
ことを特徴とする請求項１又は４に記載の筋骨格系機械インピーダンス計測方法。 - 前記推定モジュールは、
数値計算、予期し得ないノイズの累積作用の内一以上の点についても配慮して加算すべき量を算出するモジュールでもある、
ことを特徴とする請求項１、４又は５のいずれか１項に記載の筋骨格系機械インピーダンス計測方法。 - 筋骨格系に対して逐次的に摂動を与え、位置・外力の変化成分につき抽出を行う計測系に、
前記抽出の際に用いた力学的モデルに依存する係数の不確かさ成分、摂動の計測誤差の内一以上について評価を行う評価モジュールと、
前記抽出を行う抽出モジュールから入力された、今回の抽出した値と、今回の関節角ベクトルに関する行列と前回求めた機械インピーダンス行列との積との差を計算し、その各成分の大きさが、前記評価モジュールにより評価された評価値に基いて、それぞれ前記ノイズ又は不確かさの前記評価値の上限の大きさを越えた場合に、その越えた程度に応じて加減することにより、前記逐次計算をなす筋骨格系の機械インピーダンスの逐次計算を行う推定モジュールとを備えた、
ことを特徴とする筋骨格系機械インピーダンス推定装置。 - 前記推定モジュールは、
前記筋骨格系の機械インピーダンスとして、関節毎、関節相互間の弾性係数及び粘性係数を逐次計算するモジュールである、
ことを特徴とする請求項７に記載の筋骨格系機械インピーダンス推定装置。 - 前記評価モジュールは、
前記筋骨格系に対して与えた摂動による、位置・外力の変化成分を抽出した結果に含まれる、当該筋骨格系の力学的モデルに依存する係数の不確かさ、摂動の計測誤差につき一以上評価するモジュールである、
ことを特徴とする請求項７又は８に記載の筋骨格系機械インピーダンス推定装置。 - 前記抽出モジュールは、
今回の関節角ベクトルに関する行列の一成分として、各関節角の微分計算をも行うモジュールでもある、
ことを特徴とする請求項７に記載の筋骨格系機械インピーダンス計測装置。 - 前記推定モジュールは、
前記加減すべき値を、ＵＤ分解法を用いて算出するモジュールでもある、
ことを特徴とする請求項７又は１０に記載の筋骨格系機械インピーダンス計測装置。 - 前記推定モジュールは、
数値計算、予期し得ないノイズの累積作用の内一以上の点についても配慮して加算すべき量を算出するモジュールでもある、
ことを特徴とする請求項７、１０又は１１のいずれか１項に記載の筋骨格系機械インピーダンス計測装置。 - 前記計測系は、
前記筋骨格系に摂動を与える誘導外力手段と、
前記筋骨格系の位置を計測する誘導位置計測手段と、
前記筋骨格系に与えられた外力を計測する誘導外力計測手段と、
前記筋骨格系の力学的モデルを用いて前記誘導位置計測手段、前記誘導外力測定手段によりそれぞれ計測された位置、筋発生力の摂動による変化成分を抽出する誘導変化分抽出手段とを具備する、
ことを特徴とする請求項７、８、９、１０、１１又は１２のいずれか１項に記載の節筋骨格系機械インピーダンス推定装置。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002165150A JP4045487B2 (ja) | 2002-06-06 | 2002-06-06 | 筋骨格系機械インピーダンス計測方法及び装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002165150A JP4045487B2 (ja) | 2002-06-06 | 2002-06-06 | 筋骨格系機械インピーダンス計測方法及び装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2004008423A JP2004008423A (ja) | 2004-01-15 |
| JP4045487B2 true JP4045487B2 (ja) | 2008-02-13 |
Family
ID=30433057
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2002165150A Expired - Fee Related JP4045487B2 (ja) | 2002-06-06 | 2002-06-06 | 筋骨格系機械インピーダンス計測方法及び装置 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4045487B2 (ja) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR102200644B1 (ko) * | 2019-11-18 | 2021-01-12 | 울산과학기술원 | 임피던스 추정 로봇용 3차원 운동 구동 장치 |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5623759B2 (ja) * | 2010-02-26 | 2014-11-12 | 公益財団法人東京都医学総合研究所 | 筋電図信号に基づいた脳内の並列運動制御機能の同定及び評価法 |
| CN107080547A (zh) * | 2017-06-01 | 2017-08-22 | 中国科学院宁波材料技术与工程研究所 | 一种人体上肢运动特性与机械阻抗的测量系统及方法 |
-
2002
- 2002-06-06 JP JP2002165150A patent/JP4045487B2/ja not_active Expired - Fee Related
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR102200644B1 (ko) * | 2019-11-18 | 2021-01-12 | 울산과학기술원 | 임피던스 추정 로봇용 3차원 운동 구동 장치 |
Also Published As
| Publication number | Publication date |
|---|---|
| JP2004008423A (ja) | 2004-01-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Kry et al. | Interaction capture and synthesis | |
| Good et al. | Dynamic models for control system design of integrated robot and drive systems | |
| JP5761832B2 (ja) | 動作補助装置、及び動作補助装置の同調制御方法 | |
| CN105437209B (zh) | 一种具有人机交互动力学模型的外骨骼式上肢康复机器人 | |
| Zheng et al. | Modeling and simulation of planar multibody systems with revolute clearance joints considering stiction based on an LCP method | |
| CN110303521B (zh) | 机器人关节力矩信号采集系统及方法 | |
| Haghighipanah et al. | Improving position precision of a servo-controlled elastic cable driven surgical robot using unscented kalman filter | |
| Tahmasebi et al. | Dynamic parameter identification and analysis of a PHANToM haptic device | |
| Barjuei et al. | Bond graph modeling of an exoskeleton actuator | |
| Ramadurai et al. | Application of unscented kalman filter to a cable driven surgical robot: A simulation study | |
| Chang et al. | Stochastic estimation of human arm impedance under nonlinear friction in robot joints: A model study | |
| Jia et al. | Improved dynamic parameter identification method relying on proprioception for manipulators | |
| JP4045487B2 (ja) | 筋骨格系機械インピーダンス計測方法及び装置 | |
| Valera et al. | Controller–observer design and dynamic parameter identification for model-based control of an electromechanical lower-limb rehabilitation system | |
| Driels | Using passive end-point motion constraints to calibrate robot manipulators | |
| Trumić et al. | Force/torque-sensorless joint stiffness estimation in articulated soft robots | |
| Burkus et al. | A validation procedure to identify joint friction, reductor self-locking and gear backlash parameters | |
| CN114867584A (zh) | 获得机器人臂的振动特性的方法 | |
| Giraud et al. | Modeling and compensation of the internal friction torque of a travelling wave ultrasonic motor | |
| Farooq et al. | Hybrid force/position control scheme for flexible joint robot with friction between and the end-effector and the environment | |
| JP2002257673A (ja) | 機械構造系の履歴を有する非線形復元力特性の解析方法 | |
| Penzlin et al. | Torque estimation in variable stiffness actuators | |
| Wen et al. | Haptic force control based on impedance/admittance control aided by visual feedback | |
| Silawatchananai et al. | Force control of an upper limb exoskeleton for virtual reality using impedance control | |
| Hoshino et al. | Vibration control using a state observer that considers disturbances of a golf swing robot |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20040722 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20061201 |
|
| RD01 | Notification of change of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7426 Effective date: 20070625 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20070731 |
|
| RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20070815 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20071001 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20071030 |
|
| RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20071109 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20071107 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101130 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| RD05 | Notification of revocation of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7425 Effective date: 20080418 |
|
| A072 | Dismissal of procedure |
Free format text: JAPANESE INTERMEDIATE CODE: A072 Effective date: 20080812 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101130 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111130 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111130 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121130 Year of fee payment: 5 |
|
| LAPS | Cancellation because of no payment of annual fees |