JP2014089113A - 姿勢推定装置及びプログラム - Google Patents
姿勢推定装置及びプログラム Download PDFInfo
- Publication number
- JP2014089113A JP2014089113A JP2012239131A JP2012239131A JP2014089113A JP 2014089113 A JP2014089113 A JP 2014089113A JP 2012239131 A JP2012239131 A JP 2012239131A JP 2012239131 A JP2012239131 A JP 2012239131A JP 2014089113 A JP2014089113 A JP 2014089113A
- Authority
- JP
- Japan
- Prior art keywords
- posture
- unit
- value
- estimated
- vector
- 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.)
- Pending
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
- G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning
- G01C21/16—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation
- G01C21/165—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments
- G01C21/1654—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning by integrating acceleration or speed, i.e. inertial navigation combined with non-inertial navigation instruments with electromagnetic compass
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/04—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by terrestrial means
- G01C21/08—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by terrestrial means involving use of the magnetic field of the earth
Landscapes
- Remote Sensing (AREA)
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Electromagnetism (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Navigation (AREA)
- Measuring Magnetic Variables (AREA)
Abstract
【課題】真の姿勢から乖離した姿勢の推定を防止する。
【解決手段】3次元磁気センサ70と3次元加速度センサ80を備える携帯機器1の姿勢qを推定する姿勢推定装置は、携帯機器1の動きが安定しているか否かを判定する判定部220と、3次元磁気センサ70及び3次元加速度センサ80からの出力値に基づいて観測値ベクトルyを生成する観測値ベクトル生成部240と、単位時間経過後の携帯機器1の姿勢qの推定値である推定姿勢q+ kを算出するカルマンフィルタ部260とを備え、カルマンフィルタ部260は、単位時間経過後の姿勢を予測する状態遷移モデル演算部264と、前記予測結果に基づいて観測値ベクトルyの値を推定する観測モデル部281と、前記予測結果と真の姿勢との差分を推定する誤差推定部282と、判定部220の判定結果と誤差推定部282の推定結果とに基づいて推定姿勢q+ kを算出する姿勢更新部275とを具備する。
【選択図】図4
Description
本発明は姿勢推定装置及びプログラムに関する。
物体の姿勢等、動的システムの状態を推定する演算を実行する場合、1種類のセンサの出力結果に基づいて算出するよりも、地磁気センサ、加速度センサ、角速度センサ等の異種の物理量を測定する複数のセンサの出力結果を統合して算出した方が、短時間に正確な推定をすることができる。
異種の物理量を測定する複数のセンサの出力を統合し動的システムの状態を推定する方法として、カルマンフィルタを用いる方法が知られている。例えば、非特許文献1には、非線形カルマンフィルタの1種であるシグマポイントカルマンフィルタを用いて、3軸の地磁気センサ、3軸の角速度センサ、及び、3軸の加速度センサからの出力を統合することで、姿勢を推定する姿勢推定装置が開示されている。
J. L. Crassidis, and F. L. Markley, "Unscented Filtering for Spacecraft Attitude Estimation," Journal of Guidance, Control and Dynamics, 26 (2003), pp.536-542
地磁気センサや加速度センサ等のセンサは、検出対象となる成分以外の成分(ノイズ)を検出する場合がある。例えば、加速度センサは、重力加速度を検出対象とする場合であっても、当該加速度センサが振動していれば重力加速度以外の振動に係る加速度成分を検出する。また、地磁気を検出対象とする地磁気センサは、当該地磁気センサの近傍に磁界を発生させる物体が存在すれば地磁気以外の磁気成分を検出する。
姿勢推定装置が利用するセンサが、検出対象となる成分以外のノイズを検出する場合、姿勢推定装置は、当該ノイズの影響を受け、真の姿勢とは大きくかけ離れた姿勢を推定することがある。
姿勢推定装置が利用するセンサが、検出対象となる成分以外のノイズを検出する場合、姿勢推定装置は、当該ノイズの影響を受け、真の姿勢とは大きくかけ離れた姿勢を推定することがある。
本発明は、上述した点に鑑みてなされたものであり、姿勢推定装置が利用するセンサが検出対象となる成分以外のノイズ成分を検出する場合であっても、真の姿勢から乖離した姿勢の推定を防止することのできる姿勢推定装置の提供を、解決課題とする。
上述した課題を解決するため、本発明に係る姿勢推定装置は、3方向の加速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元加速度センサを含む複数のセンサを備える機器に設けられ、当該機器の姿勢を推定する姿勢推定装置であって、前記3次元加速度センサからの出力値に基づいて、前記機器の動きが安定しているか否かを判定する判定部と、前記判定部の判定結果と、前記複数のセンサのうち少なくとも一部のセンサの出力値の各々を要素とする観測値ベクトルとに基づいて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢の推定値を算出するカルマンフィルタ部と、を備え、前記カルマンフィルタ部は、前記機器の姿勢の経時的な変化を表す状態遷移モデルを用いて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢を予測する状態遷移モデル部と、前記状態遷移モデル部の予測結果と前記観測値ベクトルとに基づいて、前記状態遷移モデル部の予測結果と前記ある時刻よりも単位時間経過後の時刻における真の姿勢との差分を推定した値である推定姿勢誤差を生成する予測値補正部と、前記判定部の判定結果が肯定である場合には、前記推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出し、前記判定部の判定結果が否定である場合には、前記推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出する姿勢更新部と、を具備する、ことを特徴とする。
状態推定装置が設けられる機器の動きが安定していない場合、3次元加速度センサからの出力値には機器の振動に起因するノイズが重畳し、3次元加速度センサからの出力値の示すベクトルと重力加速度を表すベクトルとが大きく乖離する場合がある。この場合、3次元加速度センサからの出力値に基づいて機器の姿勢を推定すると、推定された姿勢のうち重力加速度方向からの傾きを表す成分は、真の姿勢の傾きから大きく乖離する可能性が高い。
この発明によれば、姿勢推定装置が設けられる機器の動きが安定していないと判定される場合に、姿勢更新部は、観測値ベクトルに基づいて生成された推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差により、状態遷移モデル部の予測結果を補正することで、単位時間経過後の姿勢の推定値を算出する。すなわち、姿勢更新部は、状態推定装置が設けられる機器の動きが安定しない場合には、推定姿勢誤差から、重力加速度方向に対する傾きを表す成分を取り除く。よって、観測値ベクトルが3次元加速度センサからの出力値を要素とする場合であって、且つ、3次元加速度センサからの出力値の示すベクトルが重力加速度とは異なる方向を示す場合であっても、3次元加速度センサからの出力値の示すベクトルと重力加速度方向とが一致していると看做して機器の姿勢の重力加速度方向からの傾きを推定することを防止することができる。
これにより、本発明に係る姿勢推定装置は、機器の姿勢の推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きから大きく乖離した傾きを表すように更新することなく、正確な姿勢を推定することが可能となる。
この発明によれば、姿勢推定装置が設けられる機器の動きが安定していないと判定される場合に、姿勢更新部は、観測値ベクトルに基づいて生成された推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差により、状態遷移モデル部の予測結果を補正することで、単位時間経過後の姿勢の推定値を算出する。すなわち、姿勢更新部は、状態推定装置が設けられる機器の動きが安定しない場合には、推定姿勢誤差から、重力加速度方向に対する傾きを表す成分を取り除く。よって、観測値ベクトルが3次元加速度センサからの出力値を要素とする場合であって、且つ、3次元加速度センサからの出力値の示すベクトルが重力加速度とは異なる方向を示す場合であっても、3次元加速度センサからの出力値の示すベクトルと重力加速度方向とが一致していると看做して機器の姿勢の重力加速度方向からの傾きを推定することを防止することができる。
これにより、本発明に係る姿勢推定装置は、機器の姿勢の推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きから大きく乖離した傾きを表すように更新することなく、正確な姿勢を推定することが可能となる。
また、本発明に係る姿勢推定装置は、3方向の磁気成分をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元磁気センサと、3方向の加速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元加速度センサと、3方向にそれぞれ延在する3軸の周りの角速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元角速度センサと、を備える機器に設けられ、当該機器の姿勢を推定する姿勢推定装置であって、前記3次元加速度センサからの出力値に基づいて、前記機器の動きが安定しているか否かを判定する判定部と、前記判定部の判定結果が肯定である場合には、前記3次元磁気センサからの出力値と前記3次元加速度センサからの出力値とを要素とする観測値ベクトルを生成し、前記判定部の判定結果が否定である場合には、前記3次元磁気センサからの出力値を要素とする観測値ベクトルを生成する観測値ベクトル生成部と、前記判定部の判定結果と、前記観測値ベクトルとに基づいて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢の推定値を算出するカルマンフィルタ部と、を備え、前記カルマンフィルタ部は、前記機器の姿勢の経時的な変化を表す状態遷移モデルに対して、ある時刻の姿勢の推定値と前記3次元角速度センサからの出力値とを適用することで、当該ある時刻よりも単位時間経過後の時刻における姿勢を予測して、これを姿勢予測値として出力する状態遷移モデル部と、前記観測値ベクトルの各要素の値と前記機器の姿勢との関係を表す観測モデルに対して、前記姿勢予測値を適用することで、当該ある時刻よりも単位時間経過後の時刻における観測値ベクトルの各要素の値を推定して、これを推定観測値ベクトルとして出力する観測モデル部と、前記推定観測値ベクトルと、前記観測値ベクトルとの差分である観測残差を算出するとともに、算出した観測残差と、前記姿勢予測値とに基づいて、前記ある時刻よりも単位時間経過後の時刻における真の姿勢と前記姿勢予測値との差分を推定した値である推定姿勢誤差を生成する誤差推定部と、前記判定部の判定結果が肯定である場合には、前記推定姿勢誤差と、前記姿勢予測値とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出し、前記判定部の判定結果が否定である場合には、前記推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差と、前記姿勢予測値とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出する姿勢更新部と、を具備する、ことを特徴とする。
この発明に係る姿勢推定装置は、姿勢推定装置が設けられる機器の動きが安定していないと判定される場合に、3次元加速度センサからの出力値を含まない観測値ベクトルを生成する。これにより、姿勢推定装置は、3次元加速度センサからの出力値にノイズが重畳している場合に、3次元加速度センサからの出力値に基づいて機器の姿勢を推定することを防止することができる。
なお、通常、地磁気の方向は地磁気の伏角をパラメータとして表現される。そのため、3次元磁気センサからの出力値の示すベクトルと地磁気を表すベクトルとが一致していると看做すことができる場合、3次元磁気センサからの出力値は地磁気の伏角を正確に反映した値となる。また、この場合、3次元磁気センサからの出力値から求められる地磁気の重力加速度成分の方向と、真の重力加速度方向とは、一致していると看做すことができる。よって、この場合、3次元磁気センサからの出力値に基づいて正確に求められる重力加速度方向を参照することで、機器の姿勢の重力加速度方向からの傾きを正確に推定することができる。
しかし、機器の近傍に磁界を発生させる物体が存在する場合、当該物体の発する磁界に起因するノイズが3次元磁気センサからの出力値に重畳する。この場合、3次元磁気センサからの出力値は地磁気の伏角を正確に反映した値とはならず、3次元磁気センサからの出力値から求められる地磁気の重力加速度成分の方向と、真の重力加速度方向とは、異なることになる。よって、この場合、3次元磁気センサからの出力値に基づいて、機器の姿勢の重力加速度方向からの傾きを正確に推定することができない。
この発明によれば、姿勢推定装置が設けられる機器の動きが安定していないと判定され、3次元磁気センサからの出力値のみに基づいて観測値ベクトルが生成される場合に、姿勢更新部は、当該観測値ベクトル(3次元磁気センサからの出力値)に基づいて生成された推定姿勢誤差から重力加速度方向に対する傾きを表す成分を取り除いた抽出推定姿勢誤差により、状態遷移モデル部の予測結果を補正し、単位時間経過後の姿勢の推定値を算出する。よって、3次元磁気センサからの出力値から求められる地磁気の重力加速度成分の方向と、真の重力加速度方向とが、異なる方向を示す場合であっても、これらを一致していると看做して機器の姿勢の重力加速度方向からの傾きを推定することを防止することができる。これにより、本発明に係る姿勢推定装置は、機器の姿勢の推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きから大きく乖離した傾きを表すように更新することなく、正確な姿勢を推定することが可能となる。
しかし、機器の近傍に磁界を発生させる物体が存在する場合、当該物体の発する磁界に起因するノイズが3次元磁気センサからの出力値に重畳する。この場合、3次元磁気センサからの出力値は地磁気の伏角を正確に反映した値とはならず、3次元磁気センサからの出力値から求められる地磁気の重力加速度成分の方向と、真の重力加速度方向とは、異なることになる。よって、この場合、3次元磁気センサからの出力値に基づいて、機器の姿勢の重力加速度方向からの傾きを正確に推定することができない。
この発明によれば、姿勢推定装置が設けられる機器の動きが安定していないと判定され、3次元磁気センサからの出力値のみに基づいて観測値ベクトルが生成される場合に、姿勢更新部は、当該観測値ベクトル(3次元磁気センサからの出力値)に基づいて生成された推定姿勢誤差から重力加速度方向に対する傾きを表す成分を取り除いた抽出推定姿勢誤差により、状態遷移モデル部の予測結果を補正し、単位時間経過後の姿勢の推定値を算出する。よって、3次元磁気センサからの出力値から求められる地磁気の重力加速度成分の方向と、真の重力加速度方向とが、異なる方向を示す場合であっても、これらを一致していると看做して機器の姿勢の重力加速度方向からの傾きを推定することを防止することができる。これにより、本発明に係る姿勢推定装置は、機器の姿勢の推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きから大きく乖離した傾きを表すように更新することなく、正確な姿勢を推定することが可能となる。
また、上述した態様において、前記判定部は、前記3次元加速度センサからの出力値を表すベクトルの大きさと、重力加速度の大きさとの差分値の絶対値が、所定の値以下である場合に、前記機器の動きが安定していると判定する、ことを特徴とすることが好ましい。
機器の動きが安定していない場合、機器に加わる加速度がノイズとして3次元加速度センサの出力値に重畳する。この態様によれば、判定部が、3次元加速度センサの出力値の大きさと重力加速度の大きさとを比較するため、機器の動きが安定しているか否かを把握することが可能となる。
機器の動きが安定していない場合、機器に加わる加速度がノイズとして3次元加速度センサの出力値に重畳する。この態様によれば、判定部が、3次元加速度センサの出力値の大きさと重力加速度の大きさとを比較するため、機器の動きが安定しているか否かを把握することが可能となる。
また、本発明に係るプログラムは、3方向の加速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元加速度センサを含む複数のセンサを備える機器に設けられ、当該機器の姿勢を推定する姿勢推定装置を制御するコンピュータのプログラムであって、前記コンピュータを、前記3次元加速度センサからの出力値に基づいて、前記機器の動きが安定しているか否かを判定する判定部と、前記判定部の判定結果と、前記複数のセンサのうち少なくとも一部のセンサの出力値の各々を要素とする観測値ベクトルとに基づいて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢の推定値を算出するカルマンフィルタ部と、して機能させ、更に、前記カルマンフィルタ部を、前記機器の姿勢の経時的な変化を表す状態遷移モデルを用いて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢を予測する状態遷移モデル部と、前記状態遷移モデル部の予測結果と前記観測値ベクトルとに基づいて、前記状態遷移モデル部の予測結果と前記ある時刻よりも単位時間経過後の時刻における真の姿勢との差分を推定した値である推定姿勢誤差を生成する予測値補正部と、前記判定部の判定結果が肯定である場合には、前記推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出し、前記判定部の判定結果が否定である場合には、前記推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出する姿勢更新部と、して機能させる、ことを特徴とする
<A.実施形態>
本発明の実施の形態について図面を参照して説明する。
本発明の実施の形態について図面を参照して説明する。
<1. 機器構成及びソフトウェア構成について>
図1は、本発明の実施形態に係る携帯機器のブロック図であり、図2は携帯機器の外観を示す斜視図である。携帯機器1は、携帯機器1の姿勢に応じて携帯機器1が備える画面(後述する表示部50)に表示されている地図などの画像を回転させることによって、画像によって表されている方位を現実空間の方位に追従させる機能を備える。この機能は、各種センサの出力に基づいてカルマンフィルタの演算を行い、携帯機器1の姿勢を推定することによって実現される。
図1は、本発明の実施形態に係る携帯機器のブロック図であり、図2は携帯機器の外観を示す斜視図である。携帯機器1は、携帯機器1の姿勢に応じて携帯機器1が備える画面(後述する表示部50)に表示されている地図などの画像を回転させることによって、画像によって表されている方位を現実空間の方位に追従させる機能を備える。この機能は、各種センサの出力に基づいてカルマンフィルタの演算を行い、携帯機器1の姿勢を推定することによって実現される。
携帯機器1は、各種の構成要素とバスを介して接続され装置全体を制御するCPU10、CPU10の作業領域として機能するRAM(蓄積部)20、姿勢推定プログラム100等の各種のプログラムやデータを記憶したROM30、通信を実行する通信部40、画像を表示する表示部50、及びGPS部60を備える。
また、携帯機器1は、地磁気等の磁気を検出して磁気データmを出力する3次元磁気センサ70、加速度を検出して加速度データaを出力する3次元加速度センサ80、及び角速度を検出して角速度データωを出力する3次元角速度センサ90を備える。
また、携帯機器1は、地磁気等の磁気を検出して磁気データmを出力する3次元磁気センサ70、加速度を検出して加速度データaを出力する3次元加速度センサ80、及び角速度を検出して角速度データωを出力する3次元角速度センサ90を備える。
表示部50は、CPU10が姿勢推定プログラム100を実行することにより推定した携帯機器1の姿勢qの推定値に基づいて、方位などを表す矢印などの画像を表示する。
GPS部60は、GPS衛星からの信号を受信して携帯機器1の位置情報(緯度、経度)を生成する。
GPS部60は、GPS衛星からの信号を受信して携帯機器1の位置情報(緯度、経度)を生成する。
3次元磁気センサ70は、X軸磁気センサ71、Y軸磁気センサ72、及びZ軸磁気センサ73を備える。各センサは、MI素子(磁気インピーダンス素子)、MR素子(磁気抵抗効果素子)などを用いて構成することができる。磁気センサI/F74は、X軸磁気センサ71、Y軸磁気センサ72、及び、Z軸磁気センサ73からの出力信号をAD変換して磁気データmを出力する。この磁気データmは、x軸、y軸およびz軸の3成分によって表されるベクトルデータである。より具体的には、磁気データmは、携帯機器1に固定された座標系において、X軸磁気センサ71からの出力値をx軸成分として表し、Y軸磁気センサ72からの出力値をy軸成分として表し、Z軸磁気センサ73からの出力値をz軸成分として表すベクトルデータである。
ところで、3次元磁気センサ70が検出する磁気データmには、地磁気Bgが含まれる。
地磁気Bgは、一般的に、磁極北に向かう水平成分と伏角により定められる鉛直成分を有する磁界であり、地上に固定された座標系において、一定の向き及び大きさを有するベクトルGBgとして表される(なお、ベクトルの符号の左上に示す添え字「G」は、当該ベクトルが地上に固定された座標系で表されたベクトルであることを意味する)。つまり、携帯機器1に固定された座標系において、地磁気Bgは、携帯機器1の姿勢qに連動して向きを変える、一定の大きさのベクトルSBg(q)として表される(なお、ベクトルの符号の左上に示す添え字「S」は、当該ベクトルが携帯機器1に固定された座標系で表されたベクトルであることを意味する)。従って、携帯機器1に固定された座標系において地磁気Bgを表すベクトルSBg(q)を算出することで、携帯機器1の姿勢qを求めることができる。
しかし、携帯機器1の近傍に、スピーカー等の磁界を発生させる物品や、各種金属等の着磁性を有する物品が存在する場合、これらの物品からの磁界(ノイズ)の影響を受けることがある。この場合、磁気データmの示す値は、地磁気Bgと、携帯機器1の近傍に存在する物品が発する磁界とが重畳した値を示す。すなわち、携帯機器1の近傍に磁界を発生させる物体が存在するときには、磁気データmに基づいて地磁気Bgの示す方向(つまり、ベクトルSBg(q)の方向)を正確に求めることは困難である場合がある。この場合、磁気データmに基づいて携帯機器1の姿勢qを正確に求めることはできない。
地磁気Bgは、一般的に、磁極北に向かう水平成分と伏角により定められる鉛直成分を有する磁界であり、地上に固定された座標系において、一定の向き及び大きさを有するベクトルGBgとして表される(なお、ベクトルの符号の左上に示す添え字「G」は、当該ベクトルが地上に固定された座標系で表されたベクトルであることを意味する)。つまり、携帯機器1に固定された座標系において、地磁気Bgは、携帯機器1の姿勢qに連動して向きを変える、一定の大きさのベクトルSBg(q)として表される(なお、ベクトルの符号の左上に示す添え字「S」は、当該ベクトルが携帯機器1に固定された座標系で表されたベクトルであることを意味する)。従って、携帯機器1に固定された座標系において地磁気Bgを表すベクトルSBg(q)を算出することで、携帯機器1の姿勢qを求めることができる。
しかし、携帯機器1の近傍に、スピーカー等の磁界を発生させる物品や、各種金属等の着磁性を有する物品が存在する場合、これらの物品からの磁界(ノイズ)の影響を受けることがある。この場合、磁気データmの示す値は、地磁気Bgと、携帯機器1の近傍に存在する物品が発する磁界とが重畳した値を示す。すなわち、携帯機器1の近傍に磁界を発生させる物体が存在するときには、磁気データmに基づいて地磁気Bgの示す方向(つまり、ベクトルSBg(q)の方向)を正確に求めることは困難である場合がある。この場合、磁気データmに基づいて携帯機器1の姿勢qを正確に求めることはできない。
3次元加速度センサ80は、X軸加速度センサ81、Y軸加速度センサ82、及びZ軸加速度センサ83を備える。各センサは、ピエゾ抵抗型、静電容量型、熱検知型などのような検知方式であってもよい。加速度センサI/F84は、各センサからの出力信号をAD変換して加速度データaを出力する。この加速度データaは、3次元加速度センサ80と一体となって運動する携帯機器1に固定された座標系における慣性力と重力との合力を、x軸、y軸及びz軸の3成分を有するベクトルとして示すデータである。
従って、携帯機器1が静止状態または等速直線運動状態にあれば、加速度データaは、携帯機器1に固定された座標系において重力加速度gの大きさと方向とを示すベクトルデータとなる。
従って、携帯機器1が静止状態または等速直線運動状態にあれば、加速度データaは、携帯機器1に固定された座標系において重力加速度gの大きさと方向とを示すベクトルデータとなる。
3次元角速度センサ90は、X軸角速度センサ91、Y軸角速度センサ92、及びZ軸角速度センサ93を備える。角速度センサI/F94は、各センサからの出力信号をAD変換して角速度データωを出力する。角速度データωは、携帯機器1に固定された座標系において、3方向に延在する各軸の周りの角速度を示すベクトルデータである。
CPU10は、ROM30に格納されている姿勢推定プログラム100を実行することによって、携帯機器1の姿勢qを推定する。すなわち、携帯機器1は、CPU10が姿勢推定プログラム100を実行することにより、姿勢推定装置として機能する。
図3は、姿勢推定装置の機能(すなわち、CPU10が姿勢推定プログラム100を実行することにより実現される機能)を表した機能ブロック図である。図3に示すように、姿勢推定装置は、加速度データa、磁気データm、及び、角速度データωに基づいて、携帯機器1の姿勢qを推定する姿勢推定部200を備える。姿勢推定部200は、判定部220、観測値ベクトル生成部240、及び、カルマンフィルタ部260を備える。
判定部220は、加速度データaに基づいて、携帯機器1の動きが安定しているか否か、すなわち、携帯機器1に大きな振動が生じていないか否かを判定する。
具体的には、判定部220は、以下の式(1)に示す条件、すなわち、「加速度データaが示すベクトルの大きさと重力加速度の大きさgとの差分値の絶対値が、所定の閾値εGよりも小さい」という条件を充足するか否かを判定する。そして、判定部220は、当該条件を充足する場合には、携帯機器1の動きが安定していると判定し、充足しない場合には、携帯機器1の動きが不安定である(携帯機器1が振動している)と判定する。
なお、各種変数、ベクトル、行列等の右下に付された添え字「k」は、当該各種変数、ベクトル、行列等の時刻T=kにおける値であることを意味する。以下では、時刻T=kにおける加速度データaの示す3次元ベクトルを加速度ベクトルakと称し、時刻T=kにおける磁気データmの示す3次元ベクトルを磁気ベクトルmkと称し、時刻T=kにおける角速度データωの示す3次元ベクトルを角速度ベクトルωkと称する。
具体的には、判定部220は、以下の式(1)に示す条件、すなわち、「加速度データaが示すベクトルの大きさと重力加速度の大きさgとの差分値の絶対値が、所定の閾値εGよりも小さい」という条件を充足するか否かを判定する。そして、判定部220は、当該条件を充足する場合には、携帯機器1の動きが安定していると判定し、充足しない場合には、携帯機器1の動きが不安定である(携帯機器1が振動している)と判定する。
なお、各種変数、ベクトル、行列等の右下に付された添え字「k」は、当該各種変数、ベクトル、行列等の時刻T=kにおける値であることを意味する。以下では、時刻T=kにおける加速度データaの示す3次元ベクトルを加速度ベクトルakと称し、時刻T=kにおける磁気データmの示す3次元ベクトルを磁気ベクトルmkと称し、時刻T=kにおける角速度データωの示す3次元ベクトルを角速度ベクトルωkと称する。
観測値ベクトル生成部240は、判定部220の判定結果に基づいて、加速度データaの各成分及び磁気データmの各成分のうちの少なくとも一部を要素とする観測値ベクトルyを生成し、これをカルマンフィルタ部260に対して出力する。
具体的には、観測値ベクトル生成部240は、判定部220の判定結果が肯定である場合(つまり、携帯機器1の動きが安定していると判定された場合)には、以下の式(2)に示す、加速度ベクトルakと磁気ベクトルmkとを要素とする6次元のベクトルである観測値ベクトルykを生成する。
一方、観測値ベクトル生成部240は、判定部220の判定結果が否定である場合(つまり、携帯機器1が振動していると判定された場合)には、以下の式(3)に示す、磁気ベクトルmkを要素とする3次元のベクトルである観測値ベクトルykを生成する。
具体的には、観測値ベクトル生成部240は、判定部220の判定結果が肯定である場合(つまり、携帯機器1の動きが安定していると判定された場合)には、以下の式(2)に示す、加速度ベクトルakと磁気ベクトルmkとを要素とする6次元のベクトルである観測値ベクトルykを生成する。
一方、観測値ベクトル生成部240は、判定部220の判定結果が否定である場合(つまり、携帯機器1が振動していると判定された場合)には、以下の式(3)に示す、磁気ベクトルmkを要素とする3次元のベクトルである観測値ベクトルykを生成する。
本実施形態において、観測値ベクトル生成部240は、3次元加速度センサ80とカルマンフィルタ部260との間に設けられたスイッチSWを含んで構成される。
スイッチSWは、判定部220の判定結果が肯定である場合にオンする。これにより、観測値ベクトル生成部240は、3次元加速度センサ80が出力する加速度データa、及び、3次元磁気センサ70が出力する磁気データmの双方を、カルマンフィルタ部260に対して伝達させる。
一方、スイッチSWは、判定部220の判定結果が否定である場合にオフする。これにより、観測値ベクトル生成部240は、3次元磁気センサ70が出力する磁気データmのみを、カルマンフィルタ部260に対して伝達させる。
スイッチSWは、判定部220の判定結果が肯定である場合にオンする。これにより、観測値ベクトル生成部240は、3次元加速度センサ80が出力する加速度データa、及び、3次元磁気センサ70が出力する磁気データmの双方を、カルマンフィルタ部260に対して伝達させる。
一方、スイッチSWは、判定部220の判定結果が否定である場合にオフする。これにより、観測値ベクトル生成部240は、3次元磁気センサ70が出力する磁気データmのみを、カルマンフィルタ部260に対して伝達させる。
カルマンフィルタ部260は、3次元角速度センサ90が出力する角速度データωと、観測値ベクトルyとに基づいて、カルマンフィルタの演算を実行し、携帯機器1の姿勢qを推定する。
上述のとおり、観測値ベクトルyは、判定部220の判定結果が肯定である場合には6次元のベクトルであり、判定部220の判定結果が否定である場合には3次元のベクトルである。そのため、本実施形態に係るカルマンフィルタ部260は、判定部220の判定結果に基づいて、観測値ベクトルyの次元を考慮しつつ姿勢qを推定する演算を実行する。
以下、カルマンフィルタ部260の動作について詳述する。
上述のとおり、観測値ベクトルyは、判定部220の判定結果が肯定である場合には6次元のベクトルであり、判定部220の判定結果が否定である場合には3次元のベクトルである。そのため、本実施形態に係るカルマンフィルタ部260は、判定部220の判定結果に基づいて、観測値ベクトルyの次元を考慮しつつ姿勢qを推定する演算を実行する。
以下、カルマンフィルタ部260の動作について詳述する。
<2. カルマンフィルタの演算について>
一般的にカルマンフィルタは、動的システムの状態の経時的な変化を表す論理上のモデルである状態遷移モデルを用いて、ある時刻(時刻T=k−1)のシステムの状態から、当該ある時刻よりも単位時間経過した後(時刻T=k)のシステムの状態を予測する。
また、カルマンフィルタは、状態遷移モデルを用いて予測したシステムの状態の予測値と観測値ベクトルykとの関係を表す論理上のモデルである観測モデルを用いて、システムの状態の予測値から、観測値ベクトルykを推定する。以下では、観測値ベクトルykの推定値を推定観測値ベクトルy− kと称する。
次に、カルマンフィルタは、推定観測値ベクトルy− kと、実際の観測値である観測値ベクトルykとの差分である観測残差ekを算出する。
そして、カルマンフィルタは、状態遷移モデルを用いて予測したシステムの状態の予測値を、観測残差ekに基づいて補正(更新)することで、ある時刻から単位時間が経過した後のシステムの状態の推定値を算出する。
一般的にカルマンフィルタは、動的システムの状態の経時的な変化を表す論理上のモデルである状態遷移モデルを用いて、ある時刻(時刻T=k−1)のシステムの状態から、当該ある時刻よりも単位時間経過した後(時刻T=k)のシステムの状態を予測する。
また、カルマンフィルタは、状態遷移モデルを用いて予測したシステムの状態の予測値と観測値ベクトルykとの関係を表す論理上のモデルである観測モデルを用いて、システムの状態の予測値から、観測値ベクトルykを推定する。以下では、観測値ベクトルykの推定値を推定観測値ベクトルy− kと称する。
次に、カルマンフィルタは、推定観測値ベクトルy− kと、実際の観測値である観測値ベクトルykとの差分である観測残差ekを算出する。
そして、カルマンフィルタは、状態遷移モデルを用いて予測したシステムの状態の予測値を、観測残差ekに基づいて補正(更新)することで、ある時刻から単位時間が経過した後のシステムの状態の推定値を算出する。
本実施形態におけるカルマンフィルタ(カルマンフィルタ部260)が推定するシステムの状態には、携帯機器1の姿勢qが含まれる。すなわち、カルマンフィルタ部260は、時刻T=k−1における携帯機器1の姿勢qの推定値である推定姿勢q+ k−1に基づいて、時刻T=kにおける姿勢qの推定値である推定姿勢q+ kを算出(推定)する。
本実施形態では、姿勢qを、クォータニオンを用いて表現する。クォータニオンとは、物体の姿勢(回転状態)を表す4次元数である。例えば、携帯機器1に固定された座標系の各軸と、地上に固定された座標系の各軸とが一致する姿勢を基準姿勢と定義し、基準姿勢をq=[0、0、0、1]Tと表現する。このとき、携帯機器1の任意の姿勢qは、基準姿勢に対して3次元の単位ベクトルρの方向に延在する軸を回転軸として角度ψだけ回転した姿勢として表現できる。すなわち、姿勢qは、クォータニオンを用いて、以下の式(4)で表される。
本実施形態では、姿勢qを、クォータニオンを用いて表現する。クォータニオンとは、物体の姿勢(回転状態)を表す4次元数である。例えば、携帯機器1に固定された座標系の各軸と、地上に固定された座標系の各軸とが一致する姿勢を基準姿勢と定義し、基準姿勢をq=[0、0、0、1]Tと表現する。このとき、携帯機器1の任意の姿勢qは、基準姿勢に対して3次元の単位ベクトルρの方向に延在する軸を回転軸として角度ψだけ回転した姿勢として表現できる。すなわち、姿勢qは、クォータニオンを用いて、以下の式(4)で表される。
また、本実施形態におけるカルマンフィルタ(カルマンフィルタ部260)が推定するシステムの状態には、状態遷移モデルを用いて求めた単位時間経過後(時刻T=k)の姿勢qの予測値である姿勢予測値q− kと、単位時間経過後の真の姿勢qとの差分(以下、単に「真の姿勢との差分」と称する場合がある)を推定した値である推定姿勢誤差δq+ kが含まれる。本実施形態では、推定姿勢誤差δq+ kは、クォータニオンを用いて表現される。
本実施形態における、カルマンフィルタの状態遷移モデルは、非線形関数fを用いて、以下に示す式(5)で与えられる。式(5)に示される状態ベクトルxkは、姿勢の変化量である推定姿勢誤差δq+ kをMRPs (modified Rodrigues parameters)を用いて表現するための3次元のベクトルであり、後述する更新後の状態ベクトルx+ kは、推定姿勢誤差δq+ kと一致する(厳密には、表現形式のみが相違する)。また、式(5)に示されるプロセスノイズwkは、3次元のベクトルであり、0を中心とするガウスノイズである。なお、プロセスノイズwkの共分散行列をQkと表す。共分散行列Qkは、3行3列の正方行列である。
式(5)に示すように、状態遷移モデルは、時刻T=k−1における角速度ベクトルωk−1に基づいて、単位時間経過後の時刻T=kにおける状態ベクトルxkを予測するものである。以下では、状態遷移モデルを用いて予測した状態ベクトルxkを、予測状態ベクトルx− kと称する。
式(5)に示すように、状態遷移モデルは、時刻T=k−1における角速度ベクトルωk−1に基づいて、単位時間経過後の時刻T=kにおける状態ベクトルxkを予測するものである。以下では、状態遷移モデルを用いて予測した状態ベクトルxkを、予測状態ベクトルx− kと称する。
状態遷移モデルは、式(5)のように、姿勢の変化量である状態ベクトルxkを予測するものとして表現されるが、以下の式(6)のように、姿勢qの経時的な変化を予測するものとしてより具体的に表現することもできる。
式(6)の右辺の演算子Ωは、式(7)で定義される。また、式(7)のI3×3は、3行3列の単位行列を表す。式(7)の演算子Γは、3次元ベクトルl=(l1,l2,l3)を用いて、式(8)で定義される。Δtは、単位時間(時刻T=k−1から時刻T=kまでの時間)である。式(7)の成分Ψkは、式(9)で定義される。
式(6)に示すように、状態遷移モデルは、時刻T=k−1における姿勢qk−1と、時刻T=k−1における角速度ベクトルωk−1に基づいて、単位時間経過後の時刻T=kにおける姿勢qkを予測するものである。上述のとおり、状態遷移モデルを用いて予測した姿勢qkが、姿勢予測値q− kである。
式(6)の右辺の演算子Ωは、式(7)で定義される。また、式(7)のI3×3は、3行3列の単位行列を表す。式(7)の演算子Γは、3次元ベクトルl=(l1,l2,l3)を用いて、式(8)で定義される。Δtは、単位時間(時刻T=k−1から時刻T=kまでの時間)である。式(7)の成分Ψkは、式(9)で定義される。
式(6)に示すように、状態遷移モデルは、時刻T=k−1における姿勢qk−1と、時刻T=k−1における角速度ベクトルωk−1に基づいて、単位時間経過後の時刻T=kにおける姿勢qkを予測するものである。上述のとおり、状態遷移モデルを用いて予測した姿勢qkが、姿勢予測値q− kである。
本実施形態における、カルマンフィルタの観測モデルは、非線形関数hを用いて、以下に示す式(10)で与えられる。式(10)に示される観測ノイズvkは、0を中心とするガウスノイズである。なお、観測ノイズvkの共分散行列をRkと表す。
上述のとおり、観測値ベクトルykは、判定部220の判定結果が肯定である場合、加速度ベクトルakと磁気ベクトルmkとを要素とする6次元のベクトルであり、判定部220の判定結果が否定の場合、磁気ベクトルmkを要素とする3次元のベクトルである。従って、判定部220の判定結果が肯定である場合、観測ノイズvkは、3次元磁気センサ70のノイズと3次元加速度センサ80のノイズとを表す6次元のベクトルであり、共分散行列Rkは6行6列の正方行列である。一方、判定部220の判定結果が否定である場合、観測ノイズvkは、3次元磁気センサ70のノイズを表す3次元のベクトルであり、共分散行列Rkは3行3列の正方行列である。
式(10)に示すように、観測モデルは、時刻T=kにおける姿勢qkを用いて、時刻T=kにおける観測値ベクトルykを推定するものである。上述のとおり、観測モデルを用いて推定した観測値ベクトルを、推定観測値ベクトルy− kと称する。
なお、非線形関数hの詳細については、後述する。
上述のとおり、観測値ベクトルykは、判定部220の判定結果が肯定である場合、加速度ベクトルakと磁気ベクトルmkとを要素とする6次元のベクトルであり、判定部220の判定結果が否定の場合、磁気ベクトルmkを要素とする3次元のベクトルである。従って、判定部220の判定結果が肯定である場合、観測ノイズvkは、3次元磁気センサ70のノイズと3次元加速度センサ80のノイズとを表す6次元のベクトルであり、共分散行列Rkは6行6列の正方行列である。一方、判定部220の判定結果が否定である場合、観測ノイズvkは、3次元磁気センサ70のノイズを表す3次元のベクトルであり、共分散行列Rkは3行3列の正方行列である。
式(10)に示すように、観測モデルは、時刻T=kにおける姿勢qkを用いて、時刻T=kにおける観測値ベクトルykを推定するものである。上述のとおり、観測モデルを用いて推定した観測値ベクトルを、推定観測値ベクトルy− kと称する。
なお、非線形関数hの詳細については、後述する。
時刻T=kにおける観測残差ekは、観測値ベクトルykと推定観測値ベクトルy− kとの差分を表すベクトルであり、以下に示す式(11)で表される。式(10)及び式(11)から明らかなように、観測残差ekは、判定部220の判定結果が肯定である場合は6次元のベクトルであり、判定部220の判定結果が否定の場合は3次元のベクトルである。
カルマンフィルタは、以下の式(12)に示すように、観測残差ek、及び、式(13)に示すカルマンゲインKkを用いて、状態ベクトルxkを、予測状態ベクトルx− kから状態ベクトルx+ kへと更新する。また、カルマンフィルタは、以下の式(14)に示すように、状態ベクトルxkの推定誤差の共分散行列Pkを更新する。ここで、P− kは、予測状態ベクトルx− kの推定誤差の共分散行列であり、P+ kは、更新後の状態ベクトルx+ kの推定誤差の共分散行列である。また、Pyy kは、観測残差ekの共分散行列であり、Pxy kは、予測状態ベクトルx− kと推定観測値ベクトルy− kとの相互共分散行列である。
カルマンフィルタは、以下の式(12)に示すように、観測残差ek、及び、式(13)に示すカルマンゲインKkを用いて、状態ベクトルxkを、予測状態ベクトルx− kから状態ベクトルx+ kへと更新する。また、カルマンフィルタは、以下の式(14)に示すように、状態ベクトルxkの推定誤差の共分散行列Pkを更新する。ここで、P− kは、予測状態ベクトルx− kの推定誤差の共分散行列であり、P+ kは、更新後の状態ベクトルx+ kの推定誤差の共分散行列である。また、Pyy kは、観測残差ekの共分散行列であり、Pxy kは、予測状態ベクトルx− kと推定観測値ベクトルy− kとの相互共分散行列である。
図4に、カルマンフィルタ部260の機能ブロック図を示す。カルマンフィルタ部260は、アンセンテッド変換を用いたシグマポイントカルマンフィルタにより、式(6)〜式(14)に示した非線形カルマンフィルタの演算を実行する。
図4に示すように、初期化部261は、加算器272が出力する更新後の状態ベクトルx+ kを単位時間だけ遅延させた状態ベクトルx+ k−1の各要素を「0」に設定して、これを出力する。すなわち、初期化部261は、実質的には、ゼロベクトルO3=[0,0,0]Tを生成してこれを出力する。
また、初期化部261は、更新後の状態ベクトルx+ kの推定誤差の共分散行列P+ kを単位時間遅延させることで、状態ベクトルx+ k−1の推定誤差の共分散行列P+ k−1を生成し、これらを出力する。
また、初期化部261は、更新後の状態ベクトルx+ kの推定誤差の共分散行列P+ kを単位時間遅延させることで、状態ベクトルx+ k−1の推定誤差の共分散行列P+ k−1を生成し、これらを出力する。
シグマポイント生成部262は、共分散行列P+ k−1、及び、ゼロベクトルO3(つまり、各要素が「0」に設定された状態ベクトルx+ k−1)に基づいて、「2dim(x)+1」個のシグマポイントχk−1(i)(但し、i=1,2,3,…,2dim(x))を生成する。ここで、シグマポイントχk−1(i)の各々は3次元ベクトルであり、シグマポイントχk−1(0)は各要素が「0」に設定された状態ベクトルx+ k−1(すなわち、ゼロベクトルO3)である。また、「dim(x)」は、状態ベクトルxkの次元を表す整数、すなわち「3」である。つまり、シグマポイント生成部262は、7個のシグマポイントを生成する。このようなシグマポイントの生成は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。
また、シグマポイント生成部262は、3次元のベクトルである「2dim(x)+1」個のシグマポイントχk−1(i)の各々の表現形式を、クォータニオンに変換することで、「2dim(x)+1」個のシグマポイントδq+ χ,k−1(i)に変換する。シグマポイントχk−1(i)からシグマポイントδq+ χ,k−1(i)へと変換する演算、すなわち、姿勢の表現形式をMRPsからクォータニオンに変換する演算は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。
シグマポイントχk−1(0)から変換されたシグマポイントδq+ χ,k−1(0)は、基準姿勢[0、0、0、1]Tを表す。また、シグマポイントδq+ χ,k−1(i)は、基準姿勢からの姿勢変化量を表す。
また、シグマポイント生成部262は、3次元のベクトルである「2dim(x)+1」個のシグマポイントχk−1(i)の各々の表現形式を、クォータニオンに変換することで、「2dim(x)+1」個のシグマポイントδq+ χ,k−1(i)に変換する。シグマポイントχk−1(i)からシグマポイントδq+ χ,k−1(i)へと変換する演算、すなわち、姿勢の表現形式をMRPsからクォータニオンに変換する演算は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。
シグマポイントχk−1(0)から変換されたシグマポイントδq+ χ,k−1(0)は、基準姿勢[0、0、0、1]Tを表す。また、シグマポイントδq+ χ,k−1(i)は、基準姿勢からの姿勢変化量を表す。
遅延部276は、姿勢更新部275が出力した時刻T=kにおける推定姿勢q+ kを単位時間だけ遅延させた推定姿勢q+ k−1を出力する。
姿勢演算部263は、以下の式(15)に示すように、推定姿勢q+ k−1とシグマポイントδq+ χ,k−1(i)とに基づいて、「2dim(x)+1」個のシグマポイントδq+ χ,k−1(i)の各々対応する、「2dim(x)+1」個の姿勢q+ χ,k−1(i)を生成する。
ここで、式(15)に表される演算子[×]は、クォータニオン積を意味する。クォータニオンは、式(4)からも明らかなように、姿勢そのものを表す以外に、姿勢の変化量(任意の軸を中心とする回転の大きさ)を表現するものでもある。そして、クォータニオン積は、一方のクォータニオンを姿勢と看做し、他方のクォータニオンを姿勢変化量と看做すことで、一方のクォータニオンで表された姿勢を、他方のクォータニオンで表された姿勢変化量だけ姿勢変化させる、公知の演算である。
上述のとおり、シグマポイントδq+ χ,k−1(i)は基準姿勢からの姿勢変化量を表す。よって、式(15)の、姿勢q+ χ,k−1(i)は、推定姿勢q+ k−1を、シグマポイントδq+ χ,k−1(i)で表される姿勢変化量だけ姿勢変化させた姿勢である。
基準姿勢を表すシグマポイントδq+ χ,k−1(0)は、姿勢変化が無いこと(姿勢変化量がゼロであること)を表す。よって、式(16)に示すように、推定姿勢q+ k−1をシグマポイントδq+ χ,k−1(0)により姿勢変化させた姿勢q+ χ,k−1(0)は推定姿勢q+ k−1に等しくなる。
ここで、式(15)に表される演算子[×]は、クォータニオン積を意味する。クォータニオンは、式(4)からも明らかなように、姿勢そのものを表す以外に、姿勢の変化量(任意の軸を中心とする回転の大きさ)を表現するものでもある。そして、クォータニオン積は、一方のクォータニオンを姿勢と看做し、他方のクォータニオンを姿勢変化量と看做すことで、一方のクォータニオンで表された姿勢を、他方のクォータニオンで表された姿勢変化量だけ姿勢変化させる、公知の演算である。
上述のとおり、シグマポイントδq+ χ,k−1(i)は基準姿勢からの姿勢変化量を表す。よって、式(15)の、姿勢q+ χ,k−1(i)は、推定姿勢q+ k−1を、シグマポイントδq+ χ,k−1(i)で表される姿勢変化量だけ姿勢変化させた姿勢である。
基準姿勢を表すシグマポイントδq+ χ,k−1(0)は、姿勢変化が無いこと(姿勢変化量がゼロであること)を表す。よって、式(16)に示すように、推定姿勢q+ k−1をシグマポイントδq+ χ,k−1(0)により姿勢変化させた姿勢q+ χ,k−1(0)は推定姿勢q+ k−1に等しくなる。
状態遷移モデル演算部264(以下、単に「状態遷移モデル部」と称する場合がある)は、姿勢演算部263が生成した「2dim(x)+1」個の姿勢q+ χ,k−1(i)のそれぞれと、角速度ベクトルωk−1とを、式(6)に示す状態遷移モデルに適用して、「2dim(x)+1」個の姿勢予測値q− χ,k(i)を生成する。
なお、姿勢q+ χ,k−1(0)を式(6)に適用して得られる姿勢予測値q− χ,k(0)は、遅延部276が出力した推定姿勢q+ k−1を式(6)に適用して得られる姿勢予測値q− kと等しい。
なお、姿勢q+ χ,k−1(0)を式(6)に適用して得られる姿勢予測値q− χ,k(0)は、遅延部276が出力した推定姿勢q+ k−1を式(6)に適用して得られる姿勢予測値q− kと等しい。
差分姿勢演算部265は、状態遷移モデル演算部264が生成した姿勢予測値q− χ,k(i)と、姿勢予測値q− χ,k(0)と、に基づいて、姿勢予測値q− χ,k(0)の示す姿勢から姿勢予測値q− χ,k(i)の示す姿勢に姿勢変化させるときの姿勢変化量を演算する。当該姿勢変化量を求める演算は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。これにより、差分姿勢演算部265は、「2dim(x)+1」個の姿勢予測値q− χ,k(i)のそれぞれを、姿勢予測値q− χ,k(0)からの姿勢の変化量である、「2dim(x)+1」個のシグマポイントδq− χ,k(i)に変換する。
そして、差分姿勢演算部265は、クォータニオンで表された「2dim(x)+1」個のシグマポイントδq− χ,k(i)のそれぞれの表現形式をMRPsに変換することで、「2dim(x)+1」個のシグマポイントχk(i)を生成する。姿勢の表現形式をクォータニオンからMRPsに変換する演算は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。
そして、差分姿勢演算部265は、クォータニオンで表された「2dim(x)+1」個のシグマポイントδq− χ,k(i)のそれぞれの表現形式をMRPsに変換することで、「2dim(x)+1」個のシグマポイントχk(i)を生成する。姿勢の表現形式をクォータニオンからMRPsに変換する演算は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。
予測状態ベクトル生成部266は、差分姿勢演算部265が生成した「2dim(x)+1」個のシグマポイントχk(i)の加重平均である予測状態ベクトルx− kを生成する。なお、当該加重平均の演算は、公知の方法を適宜用いて行えばよい。
共分散生成部267は、予測状態ベクトルx− kの推定誤差の共分散行列P− kを、予測状態ベクトルx− k、及び、シグマポイントχk(i)に基づいて生成する。この演算は、公知の方法を適宜用いて行えばよい。
観測モデル演算部268は、以下の式(17)に示すように、状態遷移モデル演算部264が生成した「2dim(x)+1」個の姿勢予測値q− χ,k(i)のそれぞれを、観測モデルに適用することで、「2dim(x)+1」個の推定観測値γk(i)を生成する。
ここで、観測モデルで用いられる非線形関数hの詳細について、説明する。
非線形関数hは、判定部220の判定結果が肯定の場合には、以下の式(18)で表され、判定部220の判定結果が否定の場合には、以下の式(19)で表される。
ここで、γmは、磁気ベクトルmkの推定値であり、以下の式(20)で表される。式(20)に現れるベクトルGBgは、地上に固定された座標系において地磁気Bgの向き及び大きさを表すベクトルであり、以下の式(21)で表される。式(21)に現れる、値rは、地磁気Bgの強さを表す値であり、値φは、地磁気Bgの伏角を表す値である。また、式(20)に現れる行列B(q)は、携帯機器1が姿勢qである場合に、地上に固定された座標系において表現されたベクトルを、携帯機器1に固定された座標系において表現するための座標変換を行う行列であり、式(22)により表される。
γaは、加速度ベクトルakの推定値であり、以下の式(23)で表される。式(23)に現れるベクトルGGRVは、地上に固定された座標系において重力加速度の向き及び大きさを表したベクトルを、重力加速度の大きさgで正規化した3次元のベクトルである。
非線形関数hは、判定部220の判定結果が肯定の場合には、以下の式(18)で表され、判定部220の判定結果が否定の場合には、以下の式(19)で表される。
ここで、γmは、磁気ベクトルmkの推定値であり、以下の式(20)で表される。式(20)に現れるベクトルGBgは、地上に固定された座標系において地磁気Bgの向き及び大きさを表すベクトルであり、以下の式(21)で表される。式(21)に現れる、値rは、地磁気Bgの強さを表す値であり、値φは、地磁気Bgの伏角を表す値である。また、式(20)に現れる行列B(q)は、携帯機器1が姿勢qである場合に、地上に固定された座標系において表現されたベクトルを、携帯機器1に固定された座標系において表現するための座標変換を行う行列であり、式(22)により表される。
γaは、加速度ベクトルakの推定値であり、以下の式(23)で表される。式(23)に現れるベクトルGGRVは、地上に固定された座標系において重力加速度の向き及び大きさを表したベクトルを、重力加速度の大きさgで正規化した3次元のベクトルである。
推定観測値ベクトル生成部269は、観測モデル演算部268が生成した「2dim(x)+1」個の推定観測値γk(i)の加重平均である推定観測値ベクトルy− kを生成する。なお、当該加重平均の演算は、公知の方法を適宜用いて行えばよい。
観測モデル演算部268、及び、推定観測値ベクトル生成部269は、観測モデルに対して姿勢予測値q− χ,k(i)を適用することで、推定観測値ベクトルy− kを生成する観測モデル部281として機能する。
減算器270は、観測値ベクトルykから、推定観測値ベクトルy− kを減算することで、観測残差ekを生成する。
カルマンゲイン生成部271は、推定観測値ベクトルy− k、「2dim(x)+1」個の推定観測値γk(i)、及び、観測ノイズvkの共分散行列Rkに基づいて、観測残差ekの共分散行列Pyy kを生成する。また、カルマンゲイン生成部271は、予測状態ベクトルx− k、推定観測値ベクトルy− k、「2dim(x)+1」個のシグマポイントχk(i)、及び、「2dim(x)+1」個の推定観測値γk(i)に基づいて、予測状態ベクトルx− kと推定観測値ベクトルy− kとの相互共分散行列Pxy kを生成する。なお、カルマンゲイン生成部271が実行する、共分散行列Pyy k、及び、相互共分散行列Pxy kを生成する具体的な演算は、公知の方法、例えば非特許文献1に記載された方法等を適宜用いて行えばよい。
次に、カルマンゲイン生成部271は、式(13)に示すように、共分散行列Pyy kと相互共分散行列Pxy kとに基づいて、カルマンゲインKkを生成する。カルマンゲインKkは、判定部220の判定結果が肯定の場合には、3行6列の行列であり、判定部220の判定結果が否定の場合には、3行3列の行列である。
そして、カルマンゲイン生成部271は、カルマンゲインKk及び観測残差ekに基づいて、式(12)の右辺第2項に示される3次元ベクトル「Kkek」を生成するとともに、
カルマンゲインKk及び共分散行列Pyy kに基づいて、式(14)の右辺第2項に示される3行3列の行列「KkPyy kKk T」を生成する。
次に、カルマンゲイン生成部271は、式(13)に示すように、共分散行列Pyy kと相互共分散行列Pxy kとに基づいて、カルマンゲインKkを生成する。カルマンゲインKkは、判定部220の判定結果が肯定の場合には、3行6列の行列であり、判定部220の判定結果が否定の場合には、3行3列の行列である。
そして、カルマンゲイン生成部271は、カルマンゲインKk及び観測残差ekに基づいて、式(12)の右辺第2項に示される3次元ベクトル「Kkek」を生成するとともに、
カルマンゲインKk及び共分散行列Pyy kに基づいて、式(14)の右辺第2項に示される3行3列の行列「KkPyy kKk T」を生成する。
加算器272は、式(12)に示すように、予測状態ベクトルx− kに、カルマンゲイン生成部271が算出したベクトル「Kkek」を加算することで、更新後の状態ベクトルx+ kを生成する。
推定姿勢誤差演算部273は、MRPsを用いて表現された更新後の状態ベクトルx+ kの表現形式をクォータニオンに変換することで、推定姿勢誤差δq+ kを生成する。
このように、差分姿勢演算部265、予測状態ベクトル生成部266、減算器270、カルマンゲイン生成部271、加算器272、及び、推定姿勢誤差演算部273は、観測値ベクトルykに基づいて観測残差ekを生成するとともに、観測残差ekと姿勢予測値q− χ,k(i)とに基づいて、推定姿勢誤差δq+ kを生成する誤差推定部282として機能する。
なお、以下では、観測モデル部281、及び、誤差推定部282を、「予測値補正部」と総称する場合がある。
なお、以下では、観測モデル部281、及び、誤差推定部282を、「予測値補正部」と総称する場合がある。
減算器274は、式(14)に示すように、予測状態ベクトルx− kの推定誤差の共分散行列P− kから、カルマンゲイン生成部271が算出した行列「KkPyy kKk T」を減算することで、更新後の状態ベクトルx+ kの推定誤差の共分散行列P+ kを生成する。
姿勢更新部275は、状態遷移モデル演算部264で生成された姿勢予測値q− χ,k(0)と、推定姿勢誤差δq+ kとに基づいて、推定姿勢q+ kを生成する。姿勢更新部275が生成する推定姿勢q+ kは、携帯機器1の姿勢qの推定値として出力される。
推定姿勢q+ kを生成するための演算は、判定部220の判定結果により異なる手順となる。そこで、以下では、判定部220の判定結果が肯定である場合、及び、否定である場合のそれぞれについて、推定姿勢q+ kを生成する手順を個別に説明する。
推定姿勢q+ kを生成するための演算は、判定部220の判定結果により異なる手順となる。そこで、以下では、判定部220の判定結果が肯定である場合、及び、否定である場合のそれぞれについて、推定姿勢q+ kを生成する手順を個別に説明する。
まず、判定部220の判定結果が肯定である場合、姿勢更新部275は、以下の式(25)に示すように、推定姿勢誤差δq+ kと姿勢予測値q− χ,k(0)とのクォータニオン積として、推定姿勢q+ kを算出する。
上述のとおり、推定姿勢誤差δq+ kは、姿勢予測値q− χ,k(0)と真の姿勢qとの差分を、3次元磁気センサ70からの出力値と3次元加速度センサ80からの出力値とを用いて推定した値である。また、姿勢予測値q− χ,k(0)は、状態遷移モデルを用いて、3次元角速度センサ90からの出力値に基づいて求めた姿勢qの予測値である。すなわち、カルマンフィルタ部260は、3次元磁気センサ70、3次元加速度センサ80、及び、3次元角速度センサ90の全てのセンサからの出力値を用いて、推定姿勢q+ kを推定する。
このように、カルマンフィルタ部260は、判定部220の判定結果が肯定である場合、3つのセンサからの出力値を統合して姿勢を推定する演算を行うため、姿勢を正確に推定することができる。
このように、カルマンフィルタ部260は、判定部220の判定結果が肯定である場合、3つのセンサからの出力値を統合して姿勢を推定する演算を行うため、姿勢を正確に推定することができる。
他方、判定部220の判定結果が否定である場合、姿勢更新部275は、以下の手順に従い推定姿勢q+ kを算出する。なお、上述のとおり、クォータニオンは、姿勢そのものを表すことも、姿勢の変化量を表すこともできる。以下の演算で登場するクォータニオンqA〜qCについても、その両方の意味で使用される。
第1に、姿勢更新部275は、以下の式(26)に従い、クォータニオンqAを生成する。クォータニオンqAは、姿勢予測値q− χ,k(0)を、推定姿勢誤差δq+ kにより姿勢変化させた姿勢を表すクォータニオンである。
第2に、姿勢更新部275は、以下の式(27)に従い、クォータニオンqBを生成する。式(27)に現れる[q− χ,k(0)]−1は、姿勢予測値q− χ,k(0)の逆クォータニオンである。
逆クォータニオンとは、クォータニオンが表す姿勢変化とは逆の姿勢変化を表す4次元数である。具体的には、例えば、クォータニオンq0が、姿勢q1から姿勢q2への姿勢変化を表している場合、逆クォータニオン[q0]−1は、姿勢q2から姿勢q1への姿勢変化を表す。
すなわち、クォータニオンqBは、クォータニオンqAとして表される姿勢を、姿勢予測値q− χ,k(0)により表される姿勢変化とは逆向きに姿勢変化させた姿勢を表す。換言すれば、クォータニオンqAは、姿勢予測値q− χ,k(0)で表される姿勢を、推定姿勢誤差δq+ kにより姿勢変化させた姿勢であるのに対して、クォータニオンqBは、基準姿勢=[0、0、0、1]Tを、推定姿勢誤差δq+ kにより姿勢変化させた姿勢(つまり、推定姿勢誤差δq+ kそのもの)を表す。
第3に、姿勢更新部275は、クォータニオンqC(以下、「抽出推定姿勢誤差」と称する場合がある)を生成する。クォータニオンqCは、図5に示すように、クォータニオンqBにより示される姿勢変化のうち、地上に固定された座標系ΣGのZ軸(つまり、重力加速度方向に延在する軸)からの傾きを表す成分を取り除き、当該Z軸を回転軸とする回転成分を抽出したクォータニオンである。
より具体的には、クォータニオンqBの各成分を、以下の式(28)で表したときに、クォータニオンqBの成分qB4が0以上の場合には、クォータニオンqCは、式(29)に従って算出され、クォータニオンqBの成分qB4が0よりも小さい場合には、クォータニオンqCは、式(30)に従って算出される。
第4に、姿勢更新部275は、式(31)に従い、姿勢予測値q− χ,k(0)とクォータニオンqCとのクォータニオン積として、推定姿勢q+ kを算出する。
第1に、姿勢更新部275は、以下の式(26)に従い、クォータニオンqAを生成する。クォータニオンqAは、姿勢予測値q− χ,k(0)を、推定姿勢誤差δq+ kにより姿勢変化させた姿勢を表すクォータニオンである。
第2に、姿勢更新部275は、以下の式(27)に従い、クォータニオンqBを生成する。式(27)に現れる[q− χ,k(0)]−1は、姿勢予測値q− χ,k(0)の逆クォータニオンである。
逆クォータニオンとは、クォータニオンが表す姿勢変化とは逆の姿勢変化を表す4次元数である。具体的には、例えば、クォータニオンq0が、姿勢q1から姿勢q2への姿勢変化を表している場合、逆クォータニオン[q0]−1は、姿勢q2から姿勢q1への姿勢変化を表す。
すなわち、クォータニオンqBは、クォータニオンqAとして表される姿勢を、姿勢予測値q− χ,k(0)により表される姿勢変化とは逆向きに姿勢変化させた姿勢を表す。換言すれば、クォータニオンqAは、姿勢予測値q− χ,k(0)で表される姿勢を、推定姿勢誤差δq+ kにより姿勢変化させた姿勢であるのに対して、クォータニオンqBは、基準姿勢=[0、0、0、1]Tを、推定姿勢誤差δq+ kにより姿勢変化させた姿勢(つまり、推定姿勢誤差δq+ kそのもの)を表す。
第3に、姿勢更新部275は、クォータニオンqC(以下、「抽出推定姿勢誤差」と称する場合がある)を生成する。クォータニオンqCは、図5に示すように、クォータニオンqBにより示される姿勢変化のうち、地上に固定された座標系ΣGのZ軸(つまり、重力加速度方向に延在する軸)からの傾きを表す成分を取り除き、当該Z軸を回転軸とする回転成分を抽出したクォータニオンである。
より具体的には、クォータニオンqBの各成分を、以下の式(28)で表したときに、クォータニオンqBの成分qB4が0以上の場合には、クォータニオンqCは、式(29)に従って算出され、クォータニオンqBの成分qB4が0よりも小さい場合には、クォータニオンqCは、式(30)に従って算出される。
第4に、姿勢更新部275は、式(31)に従い、姿勢予測値q− χ,k(0)とクォータニオンqCとのクォータニオン積として、推定姿勢q+ kを算出する。
判定部220の判定結果が否定である場合、観測値ベクトルyには、加速度ベクトルakは含まれず、磁気ベクトルmkのみが含まれる。すなわち、判定部220の判定結果が否定である場合、観測モデルにおいて、加速度ベクトルakの推定値γaは算出されず、磁気ベクトルmkの推定値γmのみが推定される。
式(23)に示す加速度ベクトルakの推定値γaは、地上に固定された座標系において重力加速度方向を表すベクトルGGRVを携帯機器1に固定された座標系において表したベクトルSGRVの推定値である。
加速度ベクトルakとベクトルSGRVとが一致していると看做すことができる場合、加速度ベクトルakより、重力加速度方向を知ることができる。よって、この場合、観測残差ek(より具体的には、加速度ベクトルakの推定値γaと加速度ベクトルakの差分を表す成分)に基づいて、カルマンフィルタ部260が推定する携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きに近づくように更新することができる。すなわち、判定部220の判定結果が肯定であり、加速度ベクトルakとベクトルSGRVとが一致していると看做すことができれば、加速度ベクトルakは、携帯機器1の姿勢qのうち重力加速度方向からの傾きを表す成分の推定に大きく貢献する。
加速度ベクトルakとベクトルSGRVとが一致していると看做すことができる場合、加速度ベクトルakより、重力加速度方向を知ることができる。よって、この場合、観測残差ek(より具体的には、加速度ベクトルakの推定値γaと加速度ベクトルakの差分を表す成分)に基づいて、カルマンフィルタ部260が推定する携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きに近づくように更新することができる。すなわち、判定部220の判定結果が肯定であり、加速度ベクトルakとベクトルSGRVとが一致していると看做すことができれば、加速度ベクトルakは、携帯機器1の姿勢qのうち重力加速度方向からの傾きを表す成分の推定に大きく貢献する。
しかし、判定部220の判定結果が否定である場合、加速度ベクトルakとベクトルSGRVとは大きく異なることがある。この場合、加速度ベクトルakに基づいて、重力加速度方向を知ることはできない。よって、この場合、カルマンフィルタ部260が推定する携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分は、真の傾きからは大きく異なる傾きを表すように更新されることになる。
これに対して本実施形態に係るカルマンフィルタ部260は、判定部220の判定結果が否定である場合、加速度ベクトルakを算出せずに、携帯機器1の姿勢qを推定する。従って、加速度ベクトルakとベクトルSGRVとは大きく異なる場合であっても、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きからは大きく異なる傾きを表すように更新されることを防止することが可能となる。
これに対して本実施形態に係るカルマンフィルタ部260は、判定部220の判定結果が否定である場合、加速度ベクトルakを算出せずに、携帯機器1の姿勢qを推定する。従って、加速度ベクトルakとベクトルSGRVとは大きく異なる場合であっても、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きからは大きく異なる傾きを表すように更新されることを防止することが可能となる。
一方、磁気ベクトルmkの推定値γmは、地上に固定された座標系において地磁気Bgの方向及び大きさを表すベクトルGBgを、携帯機器1に固定された座標系において表したベクトルSBg(q)の推定値である。そして、本実施形態では、式(21)に示すように、地磁気Bgの方向及び大きさを表すベクトルGBgは、地磁気Bgの伏角φをパラメータとして用いて表現される。そのため、磁気ベクトルmkとベクトルSBg(q)とが一致していると看做すことができる場合、磁気ベクトルmkは、地磁気Bgの伏角φを正確に反映した値となる。
また、上述のとおり、地磁気Bgの鉛直成分(重力加速度方向の成分)は、伏角φにより定められる。よって、磁気ベクトルmkが、地磁気Bgの伏角φを正確に反映した値となる場合には、磁気ベクトルmkから求められる地磁気Bgの鉛直成分の方向と重力加速度方向とが一致する。すなわち、伏角φを正確に知ることができれば、地磁気Bgの向きから、重力加速度方向を知ることができる。この場合、磁気ベクトルmkと磁気ベクトルmkの推定値γmとに基づいて、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きに近づくように更新することができる。
しかし、上述のとおり、携帯機器1の近傍に磁界を発生させる物体が存在し、磁気データmの示す値(磁気ベクトルmk)に、当該物体が発する磁界がノイズとして重畳する場合、磁気ベクトルmkとベクトルSBg(q)とが大きく乖離することがある。この場合、磁気ベクトルmkは、地磁気Bgの伏角φを正確に反映した値とはならず、磁気ベクトルmkから求められる地磁気Bgの鉛直成分の方向と重力加速度方向とが大きく乖離する。よって、この場合、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きから異なる傾きを表すように更新されることになる。すなわち、磁気ベクトルmkのみに基づいて携帯機器1の姿勢qを推定しても、姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が正しく推定されない可能性がある。
また、上述のとおり、地磁気Bgの鉛直成分(重力加速度方向の成分)は、伏角φにより定められる。よって、磁気ベクトルmkが、地磁気Bgの伏角φを正確に反映した値となる場合には、磁気ベクトルmkから求められる地磁気Bgの鉛直成分の方向と重力加速度方向とが一致する。すなわち、伏角φを正確に知ることができれば、地磁気Bgの向きから、重力加速度方向を知ることができる。この場合、磁気ベクトルmkと磁気ベクトルmkの推定値γmとに基づいて、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分を、真の傾きに近づくように更新することができる。
しかし、上述のとおり、携帯機器1の近傍に磁界を発生させる物体が存在し、磁気データmの示す値(磁気ベクトルmk)に、当該物体が発する磁界がノイズとして重畳する場合、磁気ベクトルmkとベクトルSBg(q)とが大きく乖離することがある。この場合、磁気ベクトルmkは、地磁気Bgの伏角φを正確に反映した値とはならず、磁気ベクトルmkから求められる地磁気Bgの鉛直成分の方向と重力加速度方向とが大きく乖離する。よって、この場合、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きから異なる傾きを表すように更新されることになる。すなわち、磁気ベクトルmkのみに基づいて携帯機器1の姿勢qを推定しても、姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が正しく推定されない可能性がある。
これに対して本実施形態に係る姿勢更新部275は、判定部220の判定結果が否定である場合、磁気ベクトルmkと磁気ベクトルmkの推定値γmとに基づいて生成される推定姿勢誤差δq+ kから、重力加速度方向に延在する軸からの傾きを表す成分を取り除くことで当該軸を回転軸とする回転成分を抽出したクォータニオンqCを用いて、推定姿勢q+ kを算出する。これにより、姿勢更新部275は、判定部220の判定結果が否定である場合に、磁気ベクトルmkに基づいて、携帯機器1の姿勢qのうち重力加速度方向からの傾きを表す成分を推定することを防止する。
従って、磁気ベクトルmkとベクトルSBg(q)とは大きく異なり、磁気ベクトルmkが地磁気Bgの伏角φを正確に反映した値とならない場合であっても、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きからは大きく異なる傾きを表すように更新されることを防止することが可能となる。
よって、本実施形態に係る姿勢推定装置は、センサが検出対象となる成分以外のノイズ成分を検出する場合であっても、真の姿勢とは大きくかけ離れた姿勢の推定を防止することができる。
従って、磁気ベクトルmkとベクトルSBg(q)とは大きく異なり、磁気ベクトルmkが地磁気Bgの伏角φを正確に反映した値とならない場合であっても、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きからは大きく異なる傾きを表すように更新されることを防止することが可能となる。
よって、本実施形態に係る姿勢推定装置は、センサが検出対象となる成分以外のノイズ成分を検出する場合であっても、真の姿勢とは大きくかけ離れた姿勢の推定を防止することができる。
<B.変形例>
本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を組み合わせることもできる。
本発明は上述した実施形態に限定されるものではなく、例えば、以下の変形が可能である。また、以下に示す変形例のうちの2以上の変形例を組み合わせることもできる。
(1)変形例1
上述した実施形態に係る携帯機器1は、シグマポイントカルマンフィルタを用いて非線形カルマンフィルタの演算を実行するものであったが、本発明はこのような形態に限定されるものでは無く、拡張カルマンフィルタ等、公知の非線形カルマンフィルタを適宜適用して演算を行うものでもよい。
上述した実施形態に係る携帯機器1は、シグマポイントカルマンフィルタを用いて非線形カルマンフィルタの演算を実行するものであったが、本発明はこのような形態に限定されるものでは無く、拡張カルマンフィルタ等、公知の非線形カルマンフィルタを適宜適用して演算を行うものでもよい。
(2)変形例2
上述した実施形態及び変形例において、姿勢推定部200は、観測値ベクトル生成部240を備えるが、これを備えないものであってもよい。すなわち、観測値ベクトルykは、常に、加速度ベクトルakと磁気ベクトルmkとを要素とするものであってもよい。
本発明に係るカルマンフィルタ部260は、判定部220の判定結果に基づいて、推定姿勢q+ kを算出する。すなわち、判定部220の判定結果が否定である場合には、推定姿勢誤差δq+ kのうち、重力加速度方向に延在する軸を回転軸とする回転成分を抽出したクォータニオンqCを用いて、推定姿勢q+ kを算出する。従って、仮に、判定部220の判定結果が否定であり、加速度ベクトルakとベクトルSGRVとは大きく異なる場合であっても、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きからは大きく異なる傾きを表すように更新されることを防止することができる。
上述した実施形態及び変形例において、姿勢推定部200は、観測値ベクトル生成部240を備えるが、これを備えないものであってもよい。すなわち、観測値ベクトルykは、常に、加速度ベクトルakと磁気ベクトルmkとを要素とするものであってもよい。
本発明に係るカルマンフィルタ部260は、判定部220の判定結果に基づいて、推定姿勢q+ kを算出する。すなわち、判定部220の判定結果が否定である場合には、推定姿勢誤差δq+ kのうち、重力加速度方向に延在する軸を回転軸とする回転成分を抽出したクォータニオンqCを用いて、推定姿勢q+ kを算出する。従って、仮に、判定部220の判定結果が否定であり、加速度ベクトルakとベクトルSGRVとは大きく異なる場合であっても、携帯機器1の姿勢qの推定値のうち重力加速度方向からの傾きを表す成分が、真の傾きからは大きく異なる傾きを表すように更新されることを防止することができる。
(3)変形例3
上述した実施形態及び変形例において、携帯機器1は、3次元磁気センサ70、3次元加速度センサ80、及び3次元角速度センサ90を備えるが、本発明はこのような形態に限定されるものではなく、これら3種類のセンサ以外のセンサを備えるものであってもよい。また、これら3種類のセンサのうち3次元加速度センサ80を少なくとも含む複数のセンサを備えるものであってもよい。
つまり、本発明に係る姿勢推定装置は、異種の物理量を測定する複数のセンサを備える携帯機器1に設けられ、これら複数のセンサの出力を統合しシステムの状態を推定するカルマンフィルタの演算を行うものであればよい。
上述した実施形態及び変形例において、携帯機器1は、3次元磁気センサ70、3次元加速度センサ80、及び3次元角速度センサ90を備えるが、本発明はこのような形態に限定されるものではなく、これら3種類のセンサ以外のセンサを備えるものであってもよい。また、これら3種類のセンサのうち3次元加速度センサ80を少なくとも含む複数のセンサを備えるものであってもよい。
つまり、本発明に係る姿勢推定装置は、異種の物理量を測定する複数のセンサを備える携帯機器1に設けられ、これら複数のセンサの出力を統合しシステムの状態を推定するカルマンフィルタの演算を行うものであればよい。
1…携帯機器、10…CPU、70…3次元磁気センサ、80…3次元加速度センサ、90…3次元角速度センサ、100…姿勢推定プログラム、200…姿勢推定部、220…判定部、240…観測値ベクトル生成部、260…カルマンフィルタ部、264…状態遷移モデル演算部、268…観測モデル演算部、269…推定観測値ベクトル生成部、275…姿勢更新部、281…観測モデル部、282…誤差推定部。
Claims (4)
- 3方向の加速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元加速度センサを含む複数のセンサを備える機器に設けられ、当該機器の姿勢を推定する姿勢推定装置であって、
前記3次元加速度センサからの出力値に基づいて、前記機器の動きが安定しているか否かを判定する判定部と、
前記判定部の判定結果と、前記複数のセンサのうち少なくとも一部のセンサの出力値の各々を要素とする観測値ベクトルとに基づいて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢の推定値を算出するカルマンフィルタ部と、
を備え、
前記カルマンフィルタ部は、
前記機器の姿勢の経時的な変化を表す状態遷移モデルを用いて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢を予測する状態遷移モデル部と、
前記状態遷移モデル部の予測結果と前記観測値ベクトルとに基づいて、前記状態遷移モデル部の予測結果と前記ある時刻よりも単位時間経過後の時刻における真の姿勢との差分を推定した値である推定姿勢誤差を生成する予測値補正部と、
前記判定部の判定結果が肯定である場合には、前記推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出し、前記判定部の判定結果が否定である場合には、前記推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出する姿勢更新部と、
を具備する、ことを特徴とする姿勢推定装置。 - 3方向の磁気成分をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元磁気センサと、
3方向の加速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元加速度センサと、
3方向にそれぞれ延在する3軸の周りの角速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元角速度センサと、
を備える機器に設けられ、
当該機器の姿勢を推定する姿勢推定装置であって、
前記3次元加速度センサからの出力値に基づいて、前記機器の動きが安定しているか否かを判定する判定部と、
前記判定部の判定結果が肯定である場合には、前記3次元磁気センサからの出力値と前記3次元加速度センサからの出力値とを要素とする観測値ベクトルを生成し、前記判定部の判定結果が否定である場合には、前記3次元磁気センサからの出力値を要素とする観測値ベクトルを生成する観測値ベクトル生成部と、
前記判定部の判定結果と、前記観測値ベクトルとに基づいて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢の推定値を算出するカルマンフィルタ部と、
を備え、
前記カルマンフィルタ部は、
前記機器の姿勢の経時的な変化を表す状態遷移モデルに対して、ある時刻の姿勢の推定値と前記3次元角速度センサからの出力値とを適用することで、当該ある時刻よりも単位時間経過後の時刻における姿勢を予測して、これを姿勢予測値として出力する状態遷移モデル部と、
前記観測値ベクトルの各要素の値と前記機器の姿勢との関係を表す観測モデルに対して、前記姿勢予測値を適用することで、当該ある時刻よりも単位時間経過後の時刻における観測値ベクトルの各要素の値を推定して、これを推定観測値ベクトルとして出力する観測モデル部と、
前記推定観測値ベクトルと、前記観測値ベクトルとの差分である観測残差を算出するとともに、算出した観測残差と、前記姿勢予測値とに基づいて、前記ある時刻よりも単位時間経過後の時刻における真の姿勢と前記姿勢予測値との差分を推定した値である推定姿勢誤差を生成する誤差推定部と、
前記判定部の判定結果が肯定である場合には、前記推定姿勢誤差と、前記姿勢予測値とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出し、前記判定部の判定結果が否定である場合には、前記推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差と、前記姿勢予測値とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出する姿勢更新部と、
を具備する、ことを特徴とする姿勢推定装置。 - 前記判定部は、
前記3次元加速度センサからの出力値を表すベクトルの大きさと、重力加速度の大きさとの差分値の絶対値が、所定の値以下である場合に、前記機器の動きが安定していると判定する、
ことを特徴とする、請求項1または2に記載の姿勢推定装置。 - 3方向の加速度をそれぞれ検出して、検出結果を3軸の座標系におけるベクトルデータとして順次出力する3次元加速度センサを含む複数のセンサを備える機器に設けられ、当該機器の姿勢を推定する姿勢推定装置を制御するコンピュータのプログラムであって、
前記コンピュータを、
前記3次元加速度センサからの出力値に基づいて、前記機器の動きが安定しているか否かを判定する判定部と、
前記判定部の判定結果と、前記複数のセンサのうち少なくとも一部のセンサの出力値の各々を要素とする観測値ベクトルとに基づいて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢の推定値を算出するカルマンフィルタ部と、
して機能させ、
更に、
前記カルマンフィルタ部を、
前記機器の姿勢の経時的な変化を表す状態遷移モデルを用いて、ある時刻の姿勢の推定値から、当該ある時刻よりも単位時間経過後の時刻における姿勢を予測する状態遷移モデル部と、
前記状態遷移モデル部の予測結果と前記観測値ベクトルとに基づいて、前記状態遷移モデル部の予測結果と前記ある時刻よりも単位時間経過後の時刻における真の姿勢との差分を推定した値である推定姿勢誤差を生成する予測値補正部と、
前記判定部の判定結果が肯定である場合には、前記推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出し、前記判定部の判定結果が否定である場合には、前記推定姿勢誤差から重力加速度方向に延在する軸を回転軸とする回転成分を抽出した抽出推定姿勢誤差と、前記状態遷移モデル部の予測結果とに基づいて、前記ある時刻から単位時間経過後の姿勢の推定値を算出する姿勢更新部と、
して機能させる、
ことを特徴とするプログラム。
Priority Applications (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2012239131A JP2014089113A (ja) | 2012-10-30 | 2012-10-30 | 姿勢推定装置及びプログラム |
US14/066,148 US20140122015A1 (en) | 2012-10-30 | 2013-10-29 | Attitude estimation method and apparatus |
CN201310526022.7A CN103791905B (zh) | 2012-10-30 | 2013-10-30 | 姿态估计方法和装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2012239131A JP2014089113A (ja) | 2012-10-30 | 2012-10-30 | 姿勢推定装置及びプログラム |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2014089113A true JP2014089113A (ja) | 2014-05-15 |
Family
ID=50548119
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2012239131A Pending JP2014089113A (ja) | 2012-10-30 | 2012-10-30 | 姿勢推定装置及びプログラム |
Country Status (3)
Country | Link |
---|---|
US (1) | US20140122015A1 (ja) |
JP (1) | JP2014089113A (ja) |
CN (1) | CN103791905B (ja) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6029124B1 (ja) * | 2015-05-13 | 2016-11-24 | 防衛装備庁長官 | 音源位置推定装置、方法、及びプログラム |
JP2017134515A (ja) * | 2016-01-26 | 2017-08-03 | トヨタ自動車株式会社 | 状態推定装置 |
JP2020529016A (ja) * | 2017-07-28 | 2020-10-01 | シスナヴ | 磁気センサにより測定される磁場からの方位の決定 |
JP7404227B2 (ja) | 2017-07-28 | 2023-12-25 | シスナヴ | 磁場の測定から判定された方位を特徴付ける方法および装置 |
Families Citing this family (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR20140070722A (ko) * | 2012-11-26 | 2014-06-11 | 한국전자통신연구원 | 다중속도시스템 결합 장치 및 운용 방법 |
JP2014209091A (ja) | 2013-03-25 | 2014-11-06 | ローム株式会社 | 半導体装置 |
KR101698682B1 (ko) * | 2015-08-26 | 2017-01-23 | 매그나칩 반도체 유한회사 | 지자기 센서의 출력값을 보정하는 방법 및 장치 |
CN105438499B (zh) * | 2015-11-17 | 2017-06-06 | 上海新跃仪表厂 | 绕空间轴的偏流角跟踪控制方法 |
CN109477716A (zh) * | 2016-06-02 | 2019-03-15 | 三菱电机株式会社 | 姿态估计装置、姿态估计方法和观测系统 |
CN107543546B (zh) * | 2016-06-28 | 2021-03-05 | 沈阳新松机器人自动化股份有限公司 | 一种六轴运动传感器的姿态解算方法及装置 |
CN108731675B (zh) * | 2017-04-18 | 2021-10-22 | 富士通株式会社 | 待定位物航向变化量的测量方法、测量装置和电子设备 |
CN106979780B (zh) * | 2017-05-22 | 2019-06-14 | 江苏亘德科技有限公司 | 一种无人车实时姿态测量方法 |
CN107422305B (zh) | 2017-06-06 | 2020-03-13 | 歌尔股份有限公司 | 一种麦克风阵列声源定位方法和装置 |
WO2019084804A1 (zh) * | 2017-10-31 | 2019-05-09 | 深圳市大疆创新科技有限公司 | 一种视觉里程计及其实现方法 |
CN108682059B (zh) * | 2018-06-07 | 2020-03-03 | 青岛迈金智能科技有限公司 | 一种基于三轴地磁传感器的设备姿态识别方法 |
FR3082611B1 (fr) * | 2018-06-13 | 2020-10-16 | Sysnav | Procede de calibration de magnetometres equipant un objet |
CN111352506A (zh) * | 2020-02-07 | 2020-06-30 | 联想(北京)有限公司 | 图像处理方法、装置、设备及计算机可读存储介质 |
CN111811506B (zh) * | 2020-09-15 | 2020-12-01 | 中国人民解放军国防科技大学 | 视觉/惯性里程计组合导航方法、电子设备及存储介质 |
CN113137983B (zh) * | 2021-04-30 | 2023-08-22 | 深圳市恒星物联科技有限公司 | 一种自学习的井盖姿态监测方法及监测系统 |
CN114668362B (zh) * | 2022-03-18 | 2022-11-11 | 元化智能科技(深圳)有限公司 | 无线胶囊内窥镜的定位系统、装置及计算机设备 |
CN115855072B (zh) * | 2023-03-03 | 2023-05-09 | 北京千种幻影科技有限公司 | 驾驶模拟平台的姿态估算方法、装置、设备及存储介质 |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JPS491439B1 (ja) * | 1968-09-27 | 1974-01-14 | ||
US6493631B1 (en) * | 2001-05-31 | 2002-12-10 | Mlho, Inc. | Geophysical inertial navigation system |
CN101652629A (zh) * | 2007-04-02 | 2010-02-17 | Nxp股份有限公司 | 用于定向感测的方法和系统 |
CN101782391A (zh) * | 2009-06-22 | 2010-07-21 | 北京航空航天大学 | 机动加速度辅助的扩展卡尔曼滤波航姿系统姿态估计方法 |
CN101839719A (zh) * | 2010-05-16 | 2010-09-22 | 中北大学 | 一种基于陀螺、地磁传感器的惯性测量装置 |
JP6243586B2 (ja) * | 2010-08-06 | 2017-12-06 | 任天堂株式会社 | ゲームシステム、ゲーム装置、ゲームプログラム、および、ゲーム処理方法 |
CN102252676B (zh) * | 2011-05-06 | 2014-03-12 | 微迈森惯性技术开发(北京)有限公司 | 运动姿态数据获取、人体运动姿态追踪方法及相关设备 |
-
2012
- 2012-10-30 JP JP2012239131A patent/JP2014089113A/ja active Pending
-
2013
- 2013-10-29 US US14/066,148 patent/US20140122015A1/en not_active Abandoned
- 2013-10-30 CN CN201310526022.7A patent/CN103791905B/zh not_active Expired - Fee Related
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6029124B1 (ja) * | 2015-05-13 | 2016-11-24 | 防衛装備庁長官 | 音源位置推定装置、方法、及びプログラム |
JP2017134515A (ja) * | 2016-01-26 | 2017-08-03 | トヨタ自動車株式会社 | 状態推定装置 |
JP2020529016A (ja) * | 2017-07-28 | 2020-10-01 | シスナヴ | 磁気センサにより測定される磁場からの方位の決定 |
JP7404227B2 (ja) | 2017-07-28 | 2023-12-25 | シスナヴ | 磁場の測定から判定された方位を特徴付ける方法および装置 |
Also Published As
Publication number | Publication date |
---|---|
US20140122015A1 (en) | 2014-05-01 |
CN103791905A (zh) | 2014-05-14 |
CN103791905B (zh) | 2016-08-31 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP2014089113A (ja) | 姿勢推定装置及びプログラム | |
LaValle et al. | Head tracking for the Oculus Rift | |
JP5061264B1 (ja) | 小型姿勢センサ | |
Ghobadi et al. | Robust attitude estimation from uncertain observations of inertial sensors using covariance inflated multiplicative extended Kalman filter | |
Del Rosario et al. | Quaternion-based complementary filter for attitude determination of a smartphone | |
EP2951530B1 (en) | Inertial device, method, and program | |
Zhang et al. | A novel adaptive Kalman filtering approach to human motion tracking with magnetic-inertial sensors | |
US20100250177A1 (en) | Orientation measurement of an object | |
Wu et al. | A novel approach for attitude estimation based on MEMS inertial sensors using nonlinear complementary filters | |
JP2013064695A (ja) | 状態推定装置、オフセット更新方法およびオフセット更新プログラム | |
KR101922700B1 (ko) | 가속도 센서와 지자기 센서 기반의 각속도 산출 방법 및 장치 | |
CN105841695B (zh) | 信息处理装置、信息处理方法以及记录介质 | |
Suh et al. | Quaternion-based indirect Kalman filter discarding pitch and roll information contained in magnetic sensors | |
JP2017207456A (ja) | 姿勢推定装置、姿勢推定方法、制御プログラム、および記録媒体 | |
Troni et al. | Adaptive Estimation of Measurement Bias in Three-Dimensional Field Sensors with Angular Rate Sensors: Theory and Comparative Experimental Evaluation. | |
JP7025215B2 (ja) | 測位システム及び測位方法 | |
CN111895988A (zh) | 无人机导航信息更新方法及装置 | |
JP2015179002A (ja) | 姿勢推定方法、姿勢推定装置及びプログラム | |
JP2013096724A (ja) | 状態推定装置 | |
JP2013122384A (ja) | カルマンフィルタ、及び、状態推定装置 | |
JP2013061309A (ja) | カルマンフィルタ、状態推定装置、カルマンフィルタの制御方法、及びカルマンフィルタの制御プログラム | |
Munguia et al. | An attitude and heading reference system (AHRS) based in a dual filter | |
JP2013185898A (ja) | 状態推定装置 | |
JP2015094631A (ja) | 位置算出装置及び位置算出方法 | |
Gao et al. | A novel robust Kalman filter on AHRS in the magnetic distortion environment |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
RD04 | Notification of resignation of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7424 Effective date: 20150410 |