以下、本発明を実施するための形態について図面を参照しながら説明する。
まず、図1を用いて、慣性航法により生じる2点間の距離の誤差と方位誤差について説明する。図1は、慣性航法により生じる2点間の距離の誤差と方位誤差について説明する図の一例である。
受信エリアAREA1と受信エリアAREA2は慣性装置が進入した場合に絶対位置情報を電波や音波などの無線で受信可能なエリアを示す。P1、P2は絶対位置情報に含まれる絶対位置である。絶対位置情報とは、何らかの基準(座標系)において正確であることが確認されている位置(緯度・経度・標高など)であり、天井や通路などに設置された送信機により定期的に送信されている。
P´1は絶対位置情報の受信エリアAREA1においてパラメトリックな手法(カルマンフィルタなど)により補正された位置を示し、P´2は絶対位置情報の受信エリアAREA2においてパラメトリックな手法により補正された位置を示す。円E1とE2は、パラメトリックな手法により推定された誤差の推定値(誤差情報)である。
慣性装置1を携帯した歩行者が受信エリアAREA1から受信エリアAREA2に移動する場合を考える。受信エリアAREA1で歩行者の位置はP´1に補正されておりこれを始点とする。歩行者が受信エリアAREA2に到達した際に、慣性航法のみによる移動の軌跡1012と推定された位置1001がP2と大きく乖離している。乖離に起因するのは距離と方位の推定誤差である。つまり、歩行速度を推定するための歩行速度パラメータが正しくない場合や、慣性装置1が推定した方位に誤差が含まれている場合、慣性航法のみにより推定された位置1001と、パラメトリックな手法で補正された位置P´2とは乖離する。図1では、方位の推定誤差(ヨー角補正値)Δψと距離の推定誤差ΔLが示されている。
これらの推定誤差が発生する原因は、歩行者の歩行速度を推定するための歩行速度推定パラメータが該歩行者にとって適切でないために積算して求められる位置1001に誤差が重畳するためであること、及び、慣性装置1が推定する方位(ヨー角)の補正が行われないためである。よって、これらを精度良く補正し、より正確な位置推定を行なうことのできる慣性装置1が求められる。
<本実施形態の慣性装置による測位方法の概略>
続いて、本実施形態の慣性装置1による測位方法の概略について説明する。本実施形態の測位方法は、慣性航法によって得られる測位結果(少なくとも位置座標が含まれる。慣性航法にパラメトリックな手法が用いられる場合は誤差情報も含まれる。)を絶対位置情報を受信した際にパラメトリックな手法を用いて補正し、位置座標に対する誤差(図1のE1、E2)を推定する。
慣性航法で得られる位置座標は、この補正された位置座標と誤差でリセットされる。そして、後述する「歩行条件」が成立する場合、次の処理を実行する。
1.最新の補正後の位置座標(P´2)と前回の補正後の位置座標(P´1)の2点間の距離A及び進行方向(距離Aの直線で結ぶ第1方位)を計算する。
2.前回の補正後の位置座標(P´1)と慣性航法のみで得られた位置1001の距離Bと2点間の進行方向(距離Bの直線で結ぶ第2方位)を計算する。
3.距離Aを距離Bで除算した結果を、歩行速度補正パラメータとして以後の慣性航法による位置の推定に使用する。
4.1の第1方位と2の第2方位の差分をヨー角補正値(方向補正パラメータ)Δψとして姿勢を補正する。
このような測位方法によれば、パラメトリックな手法で歩行速度補正パラメータ及び方向補正パラメータを算出し、補正された進行方向に補正された歩行速度で以降の慣性航法による測位を行うため、慣性航法により位置を精度良く推定することができる。
なお、パラメトリックな手法とは、母集団の特性を規定する母数について仮説を設けて推測する手法であり、例えば、母集団の正規性や等分散性を仮定して真値や誤差を推定する手法を言う。これに対比される言葉としてノンパラメトリックな手法があり、こちらは
母集団の分布型(母数)について一切の仮定を設けないで真値や誤差を推定する手法である。本実施形態ではパラメトリックな手法の一例としてカルマンフィルタを使用する。
<慣性装置の概略>
図2は、慣性装置1の座標系を説明する図の一例である。図2(a)は後述する絶対座標系であり、図2(b)は慣性装置1に固定されたデバイス座標系を示す。本明細書では、絶対座標系として、ユーザの歩行する平面がX軸とY軸で表され、X軸及びY軸に対し鉛直方向がZ軸で表される。X軸とY軸は垂直である。また、デバイス座標系をX´軸とY´軸とZ´軸で表す。
慣性装置1は、情報処理装置と称される倍があり、例えば、スマートデバイスと呼ばれる端末、携帯電話、スマートフォン、タブレット端末、ウェアラブルPC、PDA(Personal Digital Assistant)又はノートブックPCのように、ユーザが持ち運ぶことのできる小型のデバイスである。慣性装置1は、市販されているスマートフォン等に搭載されているような、慣性センサ(三軸の加速度センサ、三軸の角速度センサ)及び三軸の地磁気センサを備えている。慣性装置1は、各センサを利用して、加速度及び角速度の変化並びに向きを検出することができる。地磁気センサは地磁気を検出して絶対値としての方位を検出する。加速度センサは加速度の変化を検出する。角速度センサは角速度の変化を検出する。また、慣性センサとは慣性の変化を検出して変化の大きさに応じた信号を出力するセンサの総称である。
なお、慣性装置1は、上で例示したデバイス以外のデバイス(例えば、音楽プレーヤ、活動量計、腕時計等)であってもよい。また、慣性装置1は、他のデバイス(例えば、歩行ロボットや、動物に装着する首輪等)に内蔵されて提供されてもよい。これにより、鉛直方向に一定の周期で運動する多様な人又は物体の進行方向を推定することが可能となる。
<ハードウェア構成>
次に、図3を用いて、本発明の一実施形態における慣性装置1のハードウェア構成を説明する。図3は、慣性装置1が、スマートフォンのようなスマートデバイスである場合の、ハードウェア構成例である。
慣性装置1は、CPU11、RAM12、ROM13、加速度センサ14、角速度センサ15、地磁気センサ16、マイク17、スピーカ18、通信モジュール19、Bluetooth(登録商標)通信モジュール20、GPS受信モジュール21、ディスプレイ22、タッチパネル23、電池24、気圧センサ25及びバス26を有する。
CPU11は、慣性装置1の動作制御を行うプログラムを実行する。RAM12は、CPU11のワークエリア等を構成する。ROM13は、CPU11が実行するプログラム
や、プログラムの実行に必要なデータを記憶する。加速度センサ14は、慣性装置1のデバイス座標系におけるX'、Y'、Z'軸方向の加速度(センサ出力)を検出する。角速度センサ15(又はジャイロセンサ)は、慣性装置1の、デバイス座標系におけるX'、Y'、Z'軸方向の角速度(センサ出力)を検出する。地磁気センサ16は、磁北を表す三次元のベクトル(センサ出力)を出力し、慣性装置1の向きを検出する。気圧センサ25は、気圧を測定し、当該慣性装置1の高度を検出する。
マイク17は、ユーザの声等の音声を電気信号に変換する。スピーカ18は、電気信号を音声として出力する。通信モジュール19は、3Gネットワーク及び/又は無線LANに接続された他の装置と通信するための装置である。Bluetooth(登録商標)通信モジュール20は、Bluetooth(登録商標)を用いて通信するための装置である。GPS受信モジュール21は、GPS衛星やIMES(Indoor Messaging Service)によって送信される測位信号を受信するための装置である。
ディスプレイ22は、ユーザに対して画面を提示するための装置である。タッチパネル23は、ユーザからの入力を受け付ける装置である。電池24は、慣性装置1を駆動するための電力を供給する装置である。バス26は、電池24を除く各装置を相互に接続する。
なお、マイク17、スピーカ18、通信モジュール19、Bluetooth(登録商標)通信モジュール20、GPS受信モジュール21、ディスプレイ22及びタッチパネル23は、任意の構成要素である。例えば、慣性装置1が、表示画面を有しない活動量計のようなデバイスである場合には、これらの装置を備えていなくてもよい。
また、慣性装置1は、Bluetooth(登録商標)通信モジュール20の代わりに、他の規格に従って無線通信する装置(例えば、ZigBee(登録商標)通信モジュール)を備えていてもよい。
<慣性装置の全体構成>
図4は、慣性装置1の機能ブロック図の一例である。慣性装置1は、姿勢情報変換処理部100、進行方向推定部200、補正パラメータ推定処理部700、位置推定部800、姿勢情報推定部110、デバイス座標系三軸加速度取得部120、及び、測位誤差・座標推定部330を有している。慣性装置1が有するこれら各部は、図3に示したCPU11がRAM12からROM13に展開されたプログラムを実行し、図3に示す各部を制御することで実現される機能又は手段である。
補正パラメータ推定処理部700は、補正パラメータ(ヨー角補正値704と歩行速度推定パラメータ705)を算出する。詳しくは後述される。
姿勢情報推定部110は、加速度、角速度、及び地磁気の検出値を取得して後述するクォータニオン(後述する方向余弦行列又は回転行列DCM1と等価の情報を有している)を作成する。クォータニオンは姿勢情報変換処理部100に通知される。
姿勢情報変換処理部100は、方向余弦行列(又は回転行列DCM1)によりデバイス座標系の加速度を絶対座標系に変換する。また、姿勢情報変換処理部100は、ヨー角補正値704を取得して、進行方向を補正する。このようにして進行方向が補正された絶対座標系加速度ベクトル107が進行方向推定部200に通知される。
進行方向推定部200は、絶対座標系に変換された絶対座標系加速度ベクトル107を用いて歩行者の進行方向213を作成する。詳細は後述される。
位置推定部800は、進行方向推定部200が推定した進行方向213から進行方向速度ベクトルを算出し、カルマンフィルタなどのパラメトリックな手法を用いて絶対位置と誤差情報を算出する。この絶対位置で慣性航法のみによる速度と位置が求まる。また、位置推定部800は、絶対位置情報を取得するとカルマンフィルタなどのパラメトリックな手法を用いて位置を補正し誤差情報を算出する。この時、絶対位置補正フラグがTRUEとなる。絶対位置、誤差情報及び絶対位置補正フラグ(706)は補正パラメータ推定処理部700に通知される。
また、位置推定部800は、補正パラメータ推定処理部700から歩行速度推定パラメータ705を取得し、進行速度ベクトルを補正する。
デバイス座標系三軸加速度取得部120は、三軸の加速度センサからデバイス座標系の加速度を取得する。デバイス座標系の三軸加速度は座標系変換処理部106に通知される。
測位誤差・座標推定部330は送信機340などから絶対位置情報(位置を補正するための絶対位置情報なので絶対位置補正情報と称してもよい)を取得して絶対位置を抽出すると共に、測位誤差を推定する。測位誤差と絶対位置は位置推定部800に通知される。
<<姿勢情報推定部110>>
姿勢情報推定部110の機能について説明する。姿勢情報推定部110は、加速度取得部101と、角速度取得部102と、地磁気取得部103と、姿勢演算部104と、を含む。姿勢情報推定部110は、慣性装置の姿勢を示す回転行列DCM1(又は方向余弦行列、クォータニオン)を出力する。
なお、絶対座標系は、複数種類のセンサによって観測された座標値を統一的に扱うために使用される統一的な座標系であり、例えば、GPSで使用されるWGS84経緯度座標系や、UTM(ユニバーサル横メルカトル図法)など直交座標系がある。絶対座標系は、ワールド(世界)座標系とも呼ばれる。また、デバイス座標系は、ボディー座標系とも呼ばれ、慣性装置1上の1点を原点として定義し、この原点で互いに直交する三軸をそれぞれX´軸、Y´軸、Z´軸として定義した座標系のことを表す。
加速度取得部101は、加速度センサ14によって検出された三軸加速度の変化を取得する。角速度取得部102は、角速度センサ15によって検出された三軸角速度の変化を取得する。ここで取得される角速度は、加速度と同様に、デバイス座標系に固定されている。地磁気取得部103は、地磁気センサ16によって検出された磁北を表す三次元の地磁気ベクトルを取得し、慣性装置1の向きを表す情報を取得する。ここで取得される向きは、加速度と同様に、デバイス座標系に固定されている。
姿勢演算部104は、上記した、加速度取得部101、角速度取得部102及び地磁気取得部103が取得したセンサ情報を用いて、慣性装置1の現在の姿勢を演算し、その姿勢情報(回転行列)から、逆行列計算により逆回転行列を求める。この逆回転行列が後述する回転行列DCM1である。
姿勢演算部104は、一般的に知られている拡張カルマンフィルタ(非特許文献1、非特許文献2、図38−40)を利用して、慣性装置1の姿勢を表す行列を求め、その逆行列を出力する。以下、その処理について詳細に説明する。
((拡張カルマンフィルタの一般式))
図38は、一般的に知られている拡張カルマンフィルタの一般式を表す。カルマンフィルタでは、時間ステップをひとつ進めるために、時間発展と観測更新の二つの手続きが実行される。時間発展の手続きでは、前の時刻の推定状態から、現在時刻の推定状態が計算される。また、観測更新の手続では、現在時刻の観測を用い、推定値を補正して、より正確な状態が推定される。これらの手続を逐次繰り返すことによって、最適な状態変数が推定される。
図39は、拡張カルマンフィルタの時間発展の手続きにおいて使用される各変数を説明する図である。図38の時間発展の項目に含まれる(1)−(3)の式に対応して、各変数の説明が付されている。ここで、kは、離散的なステップ時間を表しており、k−1は、今の時間より1ステップ前の時間を表す。
図40は、拡張カルマンフィルタの観測更新の手続きにおいて使用される各変数を説明する図である。図38の観測更新の項目に含まれる(1)−(6)の式に対応して、各変数の説明が付されている。
((拡張カルマンフィルタの適用))
姿勢演算部104は、上記した拡張カルマンフィルタにおける「時間発展の手続き」を、角速度センサ15の出力による姿勢情報の更新に用いる(ロール角、ピッチ角、ヨー角)。また、拡張カルマンフィルタにおける「観測更新の手続き」を、加速度センサ14の出力による姿勢情報の更新に用いる(ロール角、ピッチ角)(以後、「第一の観測更新の手続き」と呼ぶ)。また、「観測更新の手続き」を、更に、地磁気センサ16の出力による姿勢情報の更新にも用いる(ヨー角)(以後、「第二の観測更新の手続き」と呼ぶ)。
このように、姿勢演算部104は、7状態の拡張カルマンフィルタを構築する。姿勢演算部104は、一つの時間発展の手続きと、二つの観測更新の手続きを、並列にループさせ、姿勢とジャイロゼロ点バイアス値の推定を行う。ここで、姿勢は、以下のようなクォータニオン(以下では、クォータニオンベクトルという)によって表される。
クォータニオンベクトルとは、4変数からなるベクトルであり、物体の姿勢を表す。ロール・ピッチ・ヨーの姿勢表現ではシンバルロックという特異点問題がある反面、クォータニオンでは、特異点が無く全ての姿勢を表すことができる。また、ジャイロゼロ点バイアス値は、三軸に対応する三変数(bx
k,by
k,bz
k)によって表される(bは定数)。
以下では、上記の三つの手続き(1)−(3)を、それぞれ説明する。
((時間発展の手続き))
まず、図5−図7を用いて、拡張カルマンフィルタにおける時間発展の手続きを説明する。姿勢演算部104は、当該手続きを実行し、後述する状態推定モデルにおいて、ジャイロ出力値を入力として、拡張カルマンフィルタの時間発展の手続きに従って時間積分を行う。そして、更新されたクォータニオンベクトルqと、誤差共分散行列Pを得る(ロール・ピッチ・ヨー角)。
図5は、本実施形態において、拡張カルマンフィルタの一般式における、(1)システムの状態推定モデルに含まれる項目を説明する図である。ここでは、現在時刻の状態推定値を、上記したクォータニオンベクトル及びジャイロゼロ点バイアス値により、図5の式(1)−1のように定義する。
また、入力量u
kを、角速度センサの出力値(ω
0xk,ω
0yk,ω
0zk)(rad/sec)を用いて、図5の式(1)−4のように定義する(bは定数)。
すなわち、(ω
xk,ω
yk,ω
zk)は、ゼロ点の値を減算し、オフセットがない角速度を表す(rad/sec)。そして、システムの状態推定モデルを、図5の式(1)−5のように表す。ここで、C
1、C
2、C
3は、任意の定数である。
図6は、本発明において、拡張カルマンフィルタの一般式における、(2)時間更新における偏微分行列(ヤコビアン)を説明する図である。図5を用いて説明したように、システムの状態推定モデルは、図5の式(1)−5によって表される。式(1)−5の右辺はfであるため、右辺の偏微分が、時間発展における偏微分行列となる。
図7は、本発明において、拡張カルマンフィルタの一般式における、(3)誤差共分散行列Pk|k-1を説明する図である。プロセスノイズQkは、システム同定によって求められる定数であり、予め算出される。
上記のプロセスノイズQ
kと、1ステップ前の時刻の誤差共分散行列と、時間発展における偏微分行列(ヤコビアン)F
k及びその転置行列F
k Tにより、現在時刻の誤差共分散行列P
k|k-1が求まる(図7の式(3)−5)。なお、現在時刻の誤差共分散行列P
k|k-1及びP
k-1|k-1は、7×7の行列となり、行列内の要素は実数となる。
姿勢演算部104は、上記のモデル及び変数の定義を用いて、拡張カルマンフィルタにおける時間発展の手続きを実行し、慣性装置1の、絶対座標系における姿勢を求め、その姿勢を表す行列の逆行列(逆回転行列)を求める。
((第一の観測更新の手続き))
次に、図8を用いて、拡張カルマンフィルタにおける第一の観測更新の手続きを説明する。姿勢演算部104は、当該手続きを実行し、加速度取得部101から得られる水平方向の角度情報と、現在のクォータニオンベクトルの水平角度情報とを比較し、その誤差分を補正する処理を行う(ロール・ピッチ角のみ)。
図8は、拡張カルマンフィルタの一般式における、
(1)観測残差
に含まれる項目を説明する図である。まず、1ステップ前の時間の観測値(ベクトル)hは、図8の式(1)−3のように表される。
上記の式に含まれる要素は、三次元の回転行列(4×4)に由来するものであり、予め決められているものとする。また、観測値(ベクトル)z
kは、図8の式(1)−2のように表される。
ここで、(a
x,a
y,a
z)は、加速度取得部101が取得した、加速度センサ14の出力値である。上記のh及びz
kにより、観測残差
また、拡張カルマンフィルタの一般式における、(2)観測更新における偏微分行列(ヤコビアン)Hkは、図8の式(1)−3に表される観測値hの偏微分を求めることによって求まる。
また、拡張カルマンフィルタの一般式における、(3)残差の共分散Skは、以下の観測ノイズ(行列)Rkと、観測更新における偏微分行列Hk、その転置行列Hk T及び現在時刻の誤差共分散行列Pk|k-1を用いて求められる。
ここで、(r
1,r
2,r
3)は、加速度センサ14のデバイス評価の結果、予め求められた分散値である。
また、拡張カルマンフィルタの一般式における、(4)カルマンゲインKkは、現在時刻の誤差共分散行列Pk|k-1と、観測更新における偏微分行列の転置行列Hk Tと、残差の共分散の逆行列Sk-1から求まる。なお、Kkは、7×3の行列であり、実数値を要素に持つ。
同様に、拡張カルマンフィルタの一般式における、(5)更新された状態推定値xk|k及び(6)更新された誤差共分散行列Pk|kは、これまでに求められた変数を用いて求めることができる。
姿勢演算部104は、上記のモデル及び変数を用いて、拡張カルマンフィルタにおける観測更新の手続きを実行し、水平方向の角度情報と、現在のクォータニオンベクトルの水平角度情報とを比較し、その誤差分を補正する(ロール・ピッチ角のみ)。
((第二の観測更新の手続き))
次に、図9を用いて、拡張カルマンフィルタにおける第二の観測更新の手続きを説明する。姿勢演算部104は、TRIADアルゴリズムに従って取得された姿勢情報から算出されるヨー角の情報を用いて、当該手続きを行い、クォータニオンベクトルのヨー角成分の補正を行う。TRIADアルゴリズムを用いた演算方法については後述する。
図9は、図8と同様に、拡張カルマンフィルタの一般式における、
(1)観測残差
に含まれる項目を説明する図である。図8の場合と同様に、1ステップ前の時間の観測値(ベクトル)hは、図9の式(1)−3のように表される。一方、図8の場合と異なり、観測値(ベクトル)z
kは、図9の式(1)−2のように表される。
上記のベクトルは、TRIADアルゴリズムによって求められた、ヨー角方向ベクトルである。
また、拡張カルマンフィルタの一般式における、(2)観測更新における偏微分行列(ヤコビアン)Hkは、第一の観測更新の手続きと同様に、1ステップ前の時間の観測値hの偏微分を求めることによって、求めることができる。
また、拡張カルマンフィルタの一般式における、(3)残差の共分散Skは、以下の観測ノイズ(行列)Rkと、観測更新における偏微分行列Hk、その転置行列Hk T及び現在時刻の誤差共分散行列Pk|k-1を用いて求められる。
ここで、(T
1,T
2,T
3)は、地磁気センサ16のデバイス評価の結果、予め求められ
た分散値である。
また、拡張カルマンフィルタの一般式における、(4)カルマンゲインKk、(5)更新された状態推定値xk|k及び(6)更新された誤差共分散行列Pk|kは、第一の観測更新の手続きと同様に、求めることができる。
<<TRIADアルゴリズムを用いた姿勢情報の計算>>
ここで、図10を用いて、姿勢演算部104が、既知のTRIADアルゴリズムにより、第二の姿勢情報を計算する例を説明する。
まず、ステップS10において、工場出荷時、あるいはユーザによる設定時に、初期化処理が実行され、加速度の鉛直下向き成分を表すベクトルと、地磁気ベクトルのリファレンスベクトルが、AccRef及びMagRefとして、それぞれ記憶される。ここで、加速度の鉛直下向き成分を表すベクトルは、以下のステップS20で説明されるように、クォータニオンから変換して求められる。また、地磁気ベクトルは、地磁気センサから入力された、磁北を表すベクトルである。当該ステップは、工場出荷時、あるいはユーザによる設定時にのみ実行される。すなわち、リファレンスベクトルは、初期化処理が実行されない限り、同じ値が保持される。
ステップS30において、AccFrameとMagFrameを用いて、3×3の行列であるMagFrameMを計算する。
ステップS32において、AccFrameとMagFrameの外積を求め、その結果を正規化することにより、AccCrossMagを求める。
ステップS34において、AccFrameとAccCrossMagの外積を求め、その結果を正規化することにより、AccCrossAcMを求める。
ステップS36において、AccFrameと、ステップS32で求めたAccCrossMagと、ステップS34で求めたAccCrossAcM(いずれも1×3の行列(ベクトル))とを組み合わせて、3x3の行列であるMagFrameMを生成する。
一方、ステップS40において、AccRefとMagRefを用いて、3×3の行列であるMagRefMを計算する。
ステップS42において、AccRefとMagRefの外積を求め、その結果を正規化することにより、MagCrossAccを求める。
ステップS44において、AccRefとMagCrossAccの外積を求め、その結果を正規化することにより、MagCrossを求める。
ステップS46において、AccRefと、ステップS42で求めたMagCrossAccと、ステップS44で求めたMagCross(いずれも1×3の行列(ベクトル))とを組み合わせて、3x3の行列であるMagRefMを生成する。
ここで、上記ステップS40(S42−S46)は、初期化処理が実行され、AccRef及びMagRefが変更されたときに実行され、以後、再度の初期化処理が実行されるまで、求めたMagRefMの値を保持して利用しても良い。
ステップS50において、MagFrameとMagRefMの内積を求める。求められた行列をmag_triad(3×3)と呼ぶ。行列mag_triadは、デバイス座標系から絶対座標系へ変換するための行列である。TRIADアルゴリズムにおいて、三つの列を、それぞれ、TRIAD1、TRIAD2、TRIAD3と呼ぶ。
ステップS60において、mag_triadの逆行列(絶対座標系からデバイス座標系への変換のための行列)を求め、これをクォータニオンに変換する。このようにして得られたクォータニオンが、第二の姿勢情報である。また、mag_triadの逆行列が回転行列DCM1である。
<姿勢情報変換処理部100>
図11を用いて図4の姿勢情報変換処理部100について説明する。姿勢情報変換処理部100は姿勢情報取得部105と座標系変換処理部106を有している。姿勢情報取得部105は、姿勢情報推定部110からDCM1を取得し、補正パラメータ推定処理部700からヨー角補正値を取得する。
そして、座標系変換処理部106は、デバイス座標系三軸加速度取得部120が取得したデバイス座標系加速度に、回転行列DCM1を乗算し、絶対座標系での三軸加速度(絶対加速度)を算出する。また、絶対座標系における三軸加速度ベクトルを回転行列DCM2で座標変換する。
図11に示すように、座標系変換処理部106は、回転行列DCM1を用いてデバイス座標系における三軸加速度ベクトルを、絶対座標系における加速度ベクトルへと座標変換する。
回転行列DCM2は方向余弦行列であり、回転行列DCM2の各要素が方向余弦と称される。回転行列DCM2は絶対座標系における三軸加速度ベクトルを、ヨー角補正値で進行方向の補正を加えた絶対座標系加速度ベクトル107へと座標変換する行列である。回転行列(DCM1とDCM2)は同じ3×3の行列である。
図12は、回転行列DCM2を説明する図の一例である。回転行列DCM2 はオイラー角を用い図12のように表わされる。図2に示したように、X軸を中心とする回転角をロール角φ、Y軸を中心とする回転角をピッチ角θ、Z軸を中心とする回転角をヨー角ψと称する。
ヨー角補正値によりヨー角ψを補正するので、座標系変換処理部106は、回転行列DCM2のロール角φとピッチ角θにゼロ(=0)を代入し、ヨー角ψに、補正パラメータ推定処理部700から通知され取得されたヨー角補正値(deg)を代入する。
これにより、回転行列DCM1とDCM2を用いた行列の積演算によって、絶対座標系において方位角に対して補正が加えられた((方位のオフセットがキャンセルされた)絶対座標系加速度ベクトル107を進行方向推定部200が取得できる。なお、ヨー角補正値の求め方については図35等において後述する。
<進行方向推定部>
図13は、進行方向推定部200の機能を説明する図の一例である。慣性装置1において進行方向推定機能を実行する進行方向推定部200は、バンドパスフィルタ201と、ピーク検出部204と、ピーク位置記憶部205と、移動加速度記憶部206と、水平成分移動速度特徴情報管理部207と、鉛直成分ピーク移動加速度取得部208と、水平成分移動速度特徴情報取得部209と、周期取得部210と、判定部211と、進行方向算出部212とを含む。進行方向推定部200は、姿勢情報変換処理部100によって得られた、絶対座標系加速度ベクトル107から、ユーザの一歩ごとの進行方向を求める。
バンドパスフィルタ201は、姿勢情報変換処理部100より出力された絶対座標系加速度ベクトル107から重力成分を除去する。通過帯域は、例えば、一般的な歩行周波数である1−3Hz程度とする。なお、通過帯域は、慣性装置1の測定対象の歩行又は移動に係る周波数に応じて、適宜変更され得る。ここでは、バンドパスフィルタ201が出力する、重力成分が除去された絶対加速度を、移動加速度202と呼ぶ。移動加速度202は、後述する、移動加速度記憶部206に記憶される。また、移動加速度の鉛直成分である、鉛直成分移動加速度203は、後述する、ピーク検出部204に渡される。
ピーク検出部204は、バンドパスフィルタ201から出力される移動加速度202のうち、鉛直成分移動加速度203の変化(時間変化)を観測し、その波形の谷ピークの位置(ピーク時刻又はピーク位置)を検出する。検出されたピーク位置は、後述する、ピーク位置記憶部205に記憶される。以下では、谷ピークの検出方法について説明する。
図14は、横軸を時間とした、鉛直成分移動加速度203(Z)と、移動加速度の水平成分(X、Y)の変化を表す波形である。このように、それぞれの波形は、移動周期(例えば、歩行周期)と整合する周期を持った波形となり、特に、鉛直成分移動加速度203は、±1m/s2程度の、水平成分に比べて大きな振幅波形が出力される。山ピークは、足が地面に着地した際に現れ、谷ピークは、片方の軸足をもう一方の足が追い越す瞬間に現れる。
図15は、歩行動作における鉛直方向の運動特性を説明する図である。歩行動作は、下肢の動作により、一般に、立脚相と遊脚相とに分類される。立脚相(Stance Phase)は、一側下肢の踵が接地して、同側足先が離れるまでの時期である。また、遊脚相(Swing Phase)は、一側下肢の足趾が地面から離れてから、同側足の踵が地面に接するまでの時期である。更に、歩行動作は、両脚支持期の有無によって特徴付けられる。一般に、歩行が緩徐になると、歩行の全体に占める両脚支持期の割合が増加し、歩行速度が速くなると少なくなることが知られている。また、両脚支持期は、走行状態になると、消失することも知られている。また、体の上下(鉛直方向)移動量及び左右(水平面方向)移動量は、立脚中期に最大となることが知られている。
立脚中期の前半は、蹴り出した片方の足が軸足追い抜く(蹴り出した片方の足が体幹直下を通り過ぎる)動作を含む。このとき、体は鉛直上方向に移動し、鉛直上向きの移動加速度が発生する。一方、立脚中期の後半は、蹴り出した片方の足の踵が接地するまでの動作を含む。このとき、体は鉛直下方向に移動し、鉛直下向きの移動加速度が発生する。
図16は、歩行動作における水平方向の運動特性を説明する図である。立脚中期の前半の移動加速度の水平方向成分は、目標となる地点への体幹の移動のために片足を蹴り出した際の加速度と、体の重心移動に伴う横揺れに起因した加速度の影響を含む。一方、立脚中期の後半の移動加速度の水平方向成分は、蹴り出した足を踵接地地点まで移動させることに伴う体幹の移動と、進行方向に対する体の横揺れで発生する加速度と、踵を接地した際の振動に伴う加速度の影響を含む。つまり、立脚中期の後半には、蹴り出しの際に必要となる純粋な移動加速度を包含していないということになる。
従って、本願では、体幹の移動のために片足を蹴り出した際の加速度が大きく反映される、立脚中期の前半の移動加速度を用いて、進行方向の推定を行う。
そこで、鉛直成分移動加速度203に着目し、シグナルの閾値判定により、谷ピークを検出することで、一歩を計測する。なお、谷ピークを用いて一歩を検出する理由は、足が着地する山ピーク位置での水平方向の移動加速度は、足が着地することに伴う振動やノイズを含む可能性があるためである。谷ピーク位置での水平方向の移動加速度は、足の着地に伴う影響を減少させ、歩行に伴う実際の加速度をより的確に反映する。
ピーク検出部204は、鉛直成分移動加速度203が、所定の閾値Thを下回り、再度、閾値Thを越える瞬間を捉えることにより、ピーク検出を行なう。ここで、所定の閾値Thを下回ったときの時刻taと、閾値Thを越えたときの時刻tbとの中間の時刻を、ピーク位置とすることができる。例えば、閾値Thは、実際に歩行した際の、鉛直方向の移動加速度の半分程度の値を設定するのが望ましい。なお、ピーク検出のために、上記以外の任意の方法を用いても良い。
また、過去のピーク位置を記録しておくことにより、現在のピーク位置と、過去のピーク位置との間の時間を表すピーク間隔を求めることができる。
ピーク位置記憶部205は、ピーク検出部204の検出したピーク位置を記憶する。ピーク位置記憶部205は、例えばリングバッファにより、ピーク位置(時間)を、最新のものから、過去の時間領域に相当するものまでを記憶する。ピーク位置記憶部205は、少なくとも、最新のピーク位置と一つ過去のピーク位置を記憶し、後に得られたピーク位置によって随時更新される。記憶されるピーク位置の数は、慣性装置1の記憶容量に応じて適宜変更されてもよい。
移動加速度記憶部206は、バンドパスフィルタ201から出力された移動加速度202を、観測された時刻情報を付加し時系列データとして記憶する。
水平成分移動速度特徴情報管理部207は、ピーク検出部204によるピーク位置の検出に応じて、そのピーク位置を中心とする所定の期間(τ)内の、水平成分の移動加速度を、各成分(X、Y)ごとに積算処理し、水平速度を算出する。この水平速度を、水平成分移動速度特徴情報と呼ぶ。水平成分移動速度特徴情報は、速さの方向と大きさの相対値を示すベクトルである。水平成分移動速度特徴情報管理部207は、水平成分速度特徴情報を、時刻情報tと共に記憶する。すなわち、水平成分移動速度特徴情報管理部207は、水平成分移動速度特徴情報を算出する算出部としての機能と、水平成分移動速度特徴情報を記憶する記憶部としての機能を有する。
図17は、図14に対応する、移動加速度の水平成分と、鉛直成分の、時間に対する変化を表す波形を示す。この例では、移動加速度の鉛直成分の波形から求められるピーク位置t1、t2、t3を中心とする、τの期間の、水平成分の移動加速度が、時間積分され、水平成分移動速度特徴ベクトルV1、V2、V3が算出される。
期間(τ)は、例えば、(tb−ta)時間以内の期間であることが望ましい。これは、全時間領域において積算を実行すると、進行方向に対する体の横揺れで発生する加速度と、踵を接地した際の振動に伴う加速度の影響を大きく含んだまま積算してしまい、進行方向の推定が正しく行なえない可能性があるためである。
水平成分移動速度特徴情報は、前述のピーク検出処理と、その後の一連の処理によって、一歩足を踏み出し、片足が軸足を追い越す瞬間に生成される。生成される特徴情報は、大きさと方向を特徴として持った水平成分移動速度特徴ベクトルVである。図18に例示されるように、水平成分移動速度特徴ベクトルは、片足が軸足を追い越す瞬間に歩行者の体が移動している方向(進行方向)と強度を表している。
鉛直成分ピーク移動加速度取得部208は、時刻tにおけるピーク位置(時間)に対応する、鉛直成分移動加速度203における移動加速度(鉛直成分ピーク移動加速度)を取得し、この加速度を、後述する、判定部211に渡す。
水平成分移動速度特徴情報取得部209は、水平成分移動速度特徴情報管理部207から、最新の水平成分移動速度特徴情報及び過去の水平成分移動速度特徴情報を取得し、取得した水平成分移動速度特徴情報を、後述する、判定部211に渡す。
周期取得部210は、ピーク位置記憶部205から複数のピーク位置を取得し、測定対象の移動周期(例えば、歩行周期)への換算処理を実行する。また、周期取得部210は、複数のピーク位置の差分を順次計算することによって、最新の移動周期及び過去の移動周期を取得することができる。周期取得部210は、取得した複数の移動周期を、後述する、判定部211に渡す。
判定部211は、図19に示される手順を実行することにより、これまでに検出された各種情報が歩行動作に由来して取得されたデータであるかどうかを判定する。歩行動作は、歩行や走行等を含んだ、測定対象による移動動作である。一方、歩行動作でない非歩行動作は、慣性装置を無作為又は意図的に振る動作や、外部環境から加速度を受けた場合(例えば、外部の移動物体により移動させられている場合)のように、測定対象による移動動作のみに起因しない動作である。以下、図19に示されるステップに沿って説明する。
まず、ステップS100において、鉛直成分ピーク移動加速度取得部208から取得した、鉛直成分ピーク移動加速度が、所定の範囲内である場合に、ステップS200に進む。そうでない場合には、ステップS600に進み、判定部211は、非歩行動作に起因すると判定する。なお、鉛直成分ピーク移動加速度についての所定の範囲は、慣性装置1の測位対象の特性(例えば、歩行者の歩行特性)に応じて、慣性装置1の製造者又はユーザにより、予め設定される。
次に、ステップS200において、水平成分移動速度特徴情報管理部207から取得した、水平成分移動速度特徴情報(ベクトル)の大きさが、所定の範囲内である場合に、ステップS300に進む。そうでない場合には、ステップS600に進み、判定部211は、非歩行動作に起因すると判定する。なお、水平成分移動速度特徴情報(ベクトル)の大きさについての所定の範囲は、慣性装置1の測位対象の特性(例えば、歩行者の歩行特性)に応じて、慣性装置1の製造者又はユーザにより、予め設定される。
次に、ステップS300において、周期取得部210から取得した、移動周期が、所定の範囲内である場合に、ステップS400に進む。そうでない場合には、ステップS600に進み、判定部211は、非歩行動作に起因すると判定する。なお、移動周期についての所定の範囲は、慣性装置1の測位対象の特性(例えば、歩行者の歩行特性)に応じて、慣性装置1の製造者又はユーザにより、予め設定される。
次に、ステップS400において、左右方向の移動速度のぶれ幅が、所定の範囲内である場合に、ステップS500に進み、判定部211は、歩行動作に起因すると判定する。そうでない場合には、ステップS600に進み、判定部211は、非歩行動作に起因すると判定する。
ここで、左右の移動速度のぶれ幅について、図20を用いて説明する。人の歩行では、右足を踏み出すと、右方向に移動速度ベクトルが発生し、左足を踏み出すと、左方向に移動速度ベクトルが発生する。本判定を行うために、水平成分移動速度特徴ベクトルがこの特性に適合するかどうか判定する。
判定部211は、まず、水平成分移動速度特徴ベクトルの始点と終点を、図20(b)のように結合する。次に、判定部211は、各Vnベクトルの中点を結ぶ直線と、各Vnベクトルの終点との距離dnを算出する。そして、判定部211は、dnが所定の範囲内であるか判定を行い、範囲内であれば、歩行動作に起因すると判定する。なお、dnについての所定の範囲は、慣性装置1の測位対象の特性(例えば、歩行者の歩行特性)に応じて、慣性装置1の製造者又はユーザにより、予め設定される。
ステップS500では、判定部211は、歩行動作に起因すると判定する。
ステップS600では、判定部211は、非歩行動作に起因すると判定する。
なお、上記したステップS100−S400における判定のうち、一部の判定が省略されてもよい。しかしながら、全ての判定がなされることで、より精度の高い進行方向の推定が可能となる。
進行方向算出部212は、上記した判定部211が、「歩行動作に起因する」であると判定した場合に、一歩ごとの進行方向213を求めるために、以下の処理を行う。
進行方向算出部212は、ユーザが、ゼロ歩目から一歩目に移るとき、水平成分移動速度特徴情報取得部209より、水平成分移動速度特徴ベクトルV0を取得する(図21(a))。また、進行方向算出部212は、ユーザが、一歩目から二歩目に移るとき、水平成分移動速度特徴情報取得部209より、水平成分移動速度特徴ベクトルV1を取得する(図21(a))。
ここで、進行方向算出部212は、ベクトルV0とベクトルV1とを規格化し、ベクトルV0'とベクトルV1'を得る(図21(b))。そして、進行方向算出部212は、それらの合成ベクトル(V0'+V1')を求め、合成ベクトルの向きを、一歩ごとの進行方向213と推定する((図21(c)))。上記処理は、ユーザが一歩踏み出すごとに、逐次実行される。
このように、慣性装置1は、一歩ごとの進行方向を推定するために、絶対座標系における加速度の、鉛直方向の加速度の谷ピーク位置を基点とし、その前後の一定の時間枠を用いて算出水平方向の速度ベクトルを用いる。これにより、歩行者の着地の際に生ずる振動等による影響を少なくすることができ、推定の精度が向上する。
更に、慣性装置1は、地磁気センサによるセンサ情報が信頼できるかどうかを評価し、センサ情報が信頼できる場合に、そのセンサ情報を、慣性装置1の姿勢を表すベクトルを補正(ヨー角成分)するために使用する。これにより、姿勢を表すベクトルを、地磁気センサ情報を用いて、より高い精度で補正することが可能となる。
<位置推定部>
図22は、位置推定部800が有する機能を説明する図の一例である。位置推定部800は現在位置を推定する。位置推定部800は、進行速度推定部300、絶対位置情報入力部400、絶対位置推定部500、速度補正部350及びマップマッチング部600を含む。以下、進行速度推定部300と、絶対位置情報入力部400と、絶対位置推定部500と、マップマッチング部600と、速度補正部350とが提供する機能について、それぞれ説明する。
<<進行速度推定部>>
図23は、進行速度推定部300の内部の詳細な機能ブロック図を表す。進行速度推定部300は、水平成分速度特徴情報取得部301と、鉛直成分ピーク加速度取得部302と、第2周期取得部303と、ぶれ幅取得部304と、進行方向取得部305と、変換部306と、進行速度波形生成部312と、進行速度波形合成部314と、進行速度波形記憶部315とを有する。
進行速度推定部300は、進行方向推定部200が算出する、水平成分速度特徴情報214、鉛直成分ピーク加速度215、周期216、ぶれ幅217及び進行方向213に基づき、測定対象の実際の進行速度を表す進行速度推定ベクトルを算出する。以下,その処理内容について説明する。
水平成分速度特徴情報取得部301は、水平成分移動速度特徴情報管理部207から水平成分速度特徴情報214(水平成分速度特徴ベクトル)を取得し、この情報を変換部306に入力する。
鉛直成分ピーク加速度取得部302は、移動加速度記憶部206から、鉛直成分ピーク加速度215を取得し、この情報を変換部306に入力する。
第2周期取得部303は、ピーク位置記憶部205に記憶された情報を用いて、測定対象の移動に係る周期216を取得し、この情報を変換部306に入力する。
ぶれ幅取得部304は、水平成分速度特徴情報214を用いて、左右の移動速度のぶれ幅(以後、ぶれ幅217とする)を取得し、この情報を変換部306に入力する。ぶれ幅217の算出方法は、図14等を用いて説明した通りである。
進行方向取得部305は、進行方向推定部200の進行方向算出部212が出力した進行方向213を取得し、この情報を変換部306に入力する。
変換部306は、入力された水平成分速度特徴情報214、鉛直成分ピーク加速度215、周期216、ぶれ幅217及び進行方向213を、それぞれ、速度パラメータPa307、強度パラメータPb308、周期パラメータPc309、ぶれ幅パラメータPd310及び進行方向パラメータPe311に変換する。当該変換は、例えば、符号214−217の各情報を、所定の規則に従って正規化する。変換部306は、入力された各種情報を、所定の範囲のパラメータに変換するために、任意の方法を用いることができる。
進行速度波形生成部312は、入力された各パラメータ307−311と、予め登録された測定対象の属性情報とをキーとして、パラメータDB313(例えば、図23)を参照し、対応する速度生成係数Ca−Ccを取得する。次に、進行速度波形生成部312は、入力された各パラメータと、取得した速度生成係数Ca−Ccを用いて、後述する式により、進行速度波形を生成する。
ここで、測定対象の属性情報は、例えば、ユーザの性別、年齢、身長等を含む。なお、測定対象の属性情報は、測定対象を特定するための任意の情報であってもよい。例えば、その他の属性情報の例として、測定対象の種別(人、動物、二足歩行ロボット等)、識別番号、型番、性質(高速、低速)等が挙げられるが、この限りではない。
図24に例示されるパラメータDB313は、測定対象の属性情報である、性別、年齢及び身長を、パラメータ307−310及び速度生成係数Ca−Ccと関連付けている。図24の例では、簡単のため、速度生成係数Ca−Ccを取得するために,進行方向パラメータPe311を用いていない(すなわち、速度生成係数は、進行方向に拠らず決定される)。しかしながら、進行方向パラメータPe311を用いて、パラメータDB313を構築しても良い。
パラメータDB313は、慣性装置1の提供者等が、属性情報によって特定される測定対象の集団に対して、パラメータ307−311で特定される動作を行わせることにより、取得したデータをもとに、予め構築するものである。なお、速度パラメータPaについては、そのノルムの値がテーブルに登録される。
進行速度波形生成部312は、入力された各種パラメータと一致するエントリが、パラメータDB313に存在しない場合には、パラメータの組み合わせと類似するエントリの速度生成係数Ca−Ccを取得する。例えば、進行速度波形生成部312は、類似するエントリを検出するために、各パラメータの二乗平均平方根(Root mean Square)の合計値が最も小さいものを選択することができる。類似するエントリを検出するために、任意の方法を用いることができる。進行速度波形生成部312は、取得した速度生成係数Ca−Ccを、進行速度波形生成部312に渡す。
進行速度波形生成部312は、速度生成係数Ca−Ccと、入力された各パラメータPa307−Pe311とを用いて、例えば、以下のような式により、時刻0≦t≦Pcにおける、測定対象の速度変化を表す進行速度波形を生成する。
なお、進行速度波形生成部312は、上記の式以外の任意の式を用いて、進行速度波形を生成しても良い。
生成された進行速度波形の例を図25に示す。ここで、θは、進行方向パラメータPeに対応する。数11によって示される式は、速度パラメータPa307による推定項(第1項)と、強度パラメータPb308での推定項(第2項)と、ぶれ幅パラメータPd310での推定項(第3項)を含む。各項は、それぞれの重みとして、速度生成係数Ca−Ccが乗算されており、複数のパラメータに基づく速度推定を行うことにより、精度の高い歩行速度波形を生成することが可能である。なお、速度生成係数Ca−Ccは、和が1となるように、予め正規化されて、パラメータDB313に格納されている。進行速度波形生成部312は、生成した進行速度波形を、進行速度波形合成部314に渡す。
進行速度波形合成部314は、進行速度波形生成部312から進行速度波形を取得すると、後述する、進行速度波形記憶部315に記憶された、過去に合成された進行速度波形を取得し、これらを合成する。このとき、進行速度波形合成部314は、対応する時間位置において、これらの進行速度波形を積算することにより合成する。進行速度波形合成部314は、進行速度波形を任意の方法で合成してもよい(例えば、複数の進行速度波形が重なった場合には各々を比較しその最大値を採用する、等)。
図26は、進行速度波形合成部314が、複数の進行速度波形を積算して合成した例を表している。進行速度波形合成部314は、合成した進行速度波形を、進行速度波形記憶部315に記憶する。
進行速度波形記憶部315は、進行速度波形合成部314により合成された進行速度波形を、時刻情報と共に記憶する。
このようにして、随時取得されるパラメータに基づいて合成された最新の進行速度波形が、進行速度波形記憶部315に記憶される。後述する絶対位置推定部500及びマップマッチング部600は、現時点における、合成された当該進行速度波形の値を参照することにより、最新の進行速度を表す進行速度情報を取得することができる。進行速度情報は、上記したように、水平方向の二成分によって表されるため、以下では、この進行速度情報を、進行速度推定ベクトル316と呼ぶ。
<<絶対位置情報入力部>>
図27は、絶対位置情報入力部400、絶対位置推定部500及びマップマッチング部600のそれぞれの内部の詳細な機能ブロック図を表す。
絶対位置情報入力部400は、第一の絶対位置取得部401と、第二の絶対位置取得部402と、測位時間計測部403と、誤差補正部404とを含む。絶対位置情報入力部400は、慣性装置1の絶対位置を表す情報と、その情報に含まれる誤差の大きさを表す誤差情報とを絶対位置推定部500に入力する。
第一の絶対位置取得部401は、例えば、Bluetooth(登録商標)通信により、外部に設けられた位置情報送信装置と通信し、慣性装置1の絶対位置を表す絶対位置情報を取得する。絶対位置情報は、例えば、緯度、経度、高度によって表される位置ベクトルX1、Y1、Z1であってもよいし、決まった地点を基点とする、位置を表す任意のベクトルであってもよい。また、第一の絶対位置取得部401は、位置情報送信装置から、上記情報の誤差を表す誤差情報σ1も取得する。誤差情報σ1は、上記の位置情報送信装置から取得される、絶対位置についての誤差共分散行列である。誤差情報σ1は、例えば、位置情報送信装置と慣性装置1との通信に係る電波強度に応じて決定される誤差値を含む。例えば、電波強度が低い値を示す場合には、より誤差の大きな値を有する誤差共分散行列が取得される。誤差情報は、予め慣性装置1に記憶されていてもよいし、位置情報送信装置から送信されてもよい。取得した位置ベクトル及び誤差情報は、後述する、絶対位置推定部500の第一の観測更新処理部502に渡される。
第二の絶対位置取得部402は、第一の絶対位置取得部401と異なる通信手段(例えば、GPS又はIMES(Indoor Messaging System))により、絶対位置情報X2、Y2、Z2と誤差情報σ2を取得する。なお、絶対位置取得部は、少なくとも一つ存在していればよく、その数は任意である。取得した位置ベクトル及び誤差情報は、後述する、絶対位置推定部500の第二の観測更新処理部503に渡される。
第一の絶対位置取得部401と第二の絶対位置取得部402が取得する絶対位置情報と誤差情報を図22では絶対位置情報405、誤差情報406で示した。
測位時間計測部403は、第一の絶対位置取得部401及び第二の絶対位置取得部402が、それぞれ絶対位置情報を取得する間隔を計測し、その間隔を、後述する誤差補正部404に渡す。
誤差補正部404は、測位時間計測部403から受け取った間隔が、所定の間隔であるかどうかを判定し、その間隔の大きさに応じて、各絶対位置取得部から出力される誤差共分散行列(σ1又はσ2)の各分散値が大きくなるよう補正をかける。このために、誤差補正部404は、予め用意した、間隔[秒]と補正量(分散値に乗算する値)とを関連付けるテーブルを利用することができる。あるいは、誤差補正部404は、間隔の大きさの閾値を設け、閾値を超えると、誤差共分散行列に対して所定の補正をかけてもよい。
特に、GPSやIMESの誤差情報は、マルチパスの影響を考慮して生成されないため、誤差情報に対する信頼度が、電波状況によって変化する可能性がある。信号量と雑音量の比率を表すS/N比と、測位間隔の大きさとの間には、図28に示すような相関関係がある。そこで、誤差補正部404は、測位間隔が所定の間隔を超える場合には、測位誤差がより大きくなるよう補正をかけることにより、誤差情報の信頼度の変化を吸収することができる。誤差情報が直接、送信されない場合については後述する。
なお、外部の位置情報送信装置は、上記の通信方式のほか、赤外線、無線LAN、可視光通信、カメラを用いた測位手段等の手段により、絶対位置情報及び誤差情報を送信してもよい。慣性装置1は、対応する受信手段により、絶対位置情報及び誤差情報を受信することができる。受信した絶対位置情報及び誤差情報は、上記した絶対位置取得部と同様に、絶対位置推定部500の観測更新処理部に入力される。絶対位置取得部及び観測更新処理部の数は、任意である。
<<絶対位置推定部>>
絶対位置推定部500は、時間発展処理部501と、第一の観測更新処理部502と、第二の観測更新処理部503と、第三の観測更新処理部504とを含む。
絶対位置推定部500は、進行速度推定ベクトル316と、進行速度推定ベクトル316と共に提供される進行速度測定誤差情報317とを用いて、現在位置と、その誤差を表す誤差情報(現在位置と誤差情報505)を求める。進行速度測定誤差情報317は、進行速度推定ベクトル316の誤差を表す誤差共分散行列σvであり、システム同定によって求めた固定のデータが使用される。あるいは、進行速度測定誤差情報317は、進行速度に応じて用意された、複数のデータが使用されてもよい。
また、絶対位置推定部500は、絶対位置情報入力部400が出力した絶対位置情報(位置ベクトル)及び誤差情報(誤差共分散行列)とを用いて、現在位置と誤差情報505を更新する。更に、絶対位置推定部500は、マップマッチング部600の出力する位置情報(位置ベクトル)と誤差情報(誤差共分散行列)とを用いて、現在位置と誤差情報505を更新する。
拡張カルマンフィルタを用いて、現在位置と誤差情報505を算出又は更新する。ここでは、一つの時間発展の手続き(時間発展処理部501)と、三つの観測更新の手続き(第一の観測更新処理部502、第二の観測更新処理部503、第三の観測更新処理部504)が、並列に実行される。各手続きにおいて用いられる変数及びモデルを以下で説明する。
時間発展処理部501は、図29−図31に示す変数及びモデルの定義に従って、拡張カルマンフィルタの時間発展の手続きを実行し、慣性装置1の現在位置と誤差情報505を算出又は更新する。ここで、拡張カルマンフィルタにおける変数とモデルを、図29に示すように定義する。すなわち、今の時刻の状態推定値は、図29の式(1)−1のように、三次元の位置ベクトルによって定義される。また、入力量は、図29の式(1)−4のように、進行速度推定部300によって出力された、進行速度推定ベクトル316によって定義される。また、システムの状態推定モデルは、図29の式(1)−4のように定義される。
また、図30に示す通り、時間発展における偏微分行列(ヤコビアン)は、システムの状態推定モデルにおける右辺の偏微分が、そのまま偏微分行列となる。
また、図31に示す通り、プロセスノイズQkは、システム同定によって予め求められる、進行速度測定誤差情報317(σv)であり、予め算出された定数である。なお、現在時刻の誤差共分散行列Pk|k-1及びPk-1|k-1は、3×3の行列となり、行列内の要素は実数となる。
第一の観測更新処理部502は、拡張カルマンフィルタの観測更新の手続きを実行し、慣性装置1の現在位置と誤差情報505を算出又は更新する。ここで、拡張カルマンフィルタにおける変数を、図32を用いて説明する。すなわち、1ステップ前の時間の観測値は、図32の式(1)−3のように、三次元の位置ベクトルとなる。また、観測値は、図32の式(1)−2のように、第一の絶対位置取得部401が出力する位置ベクトル(絶対位置情報)である。上記のh及びzkにより、観測残差
また、拡張カルマンフィルタの一般式における(2)観測更新における偏微分行列(ヤコビアン)Hkは、図32の式(1)−3に表される観測値hの偏微分を求めることによって求まる。
また、拡張カルマンフィルタの一般式における(3)残差の共分散Skは、第一の絶対位置取得部401の出力した誤差情報である観測ノイズ(行列)Rkと、観測更新における偏微分行列Hk、その転置行列HkT及び現在時刻の誤差共分散行列Pk|k-1を用いて求まる。
なお、r
1はX軸、r
2はY軸、r
3はZ軸の成分を表す。
また、拡張カルマンフィルタの一般式における、(4)カルマンゲインKkは、現在時刻の誤差共分散行列Pk|k-1と、観測更新における偏微分行列の転置行列HkTと、残差の共分散の逆行列Sk-1から求まる。
同様に、拡張カルマンフィルタの一般式における、(5)更新された状態推定値xk|k及び(6)更新された誤差共分散行列Pk|kは、これまでに求められた変数を用いて求めることができる。
第二の観測更新処理部503は、第一の観測更新処理部502と同様に、拡張カルマンフィルタの観測更新の手続きを実行し、慣性装置1の現在位置と誤差情報505を算出又は更新する。拡張カルマンフィルタにおける変数は、観測値zk及び観測ノイズRkが、第二の絶対位置取得部402によって出力される位置ベクトル及び誤差情報である点を除いて、第一の観測更新処理部502と同様である。
第三の観測更新処理部504は、第一の観測更新処理部502と同様に、拡張カルマンフィルタの観測更新の手続きを実行し、慣性装置1の現在位置と誤差情報505を算出又は更新する。拡張カルマンフィルタにおける変数は、観測値zk及び観測ノイズRkが、後述する、マップマッチング部600によって出力される位置ベクトル及び誤差情報である点を除いて、 第一の観測更新処理部502と同様である。
絶対位置推定部500は、以上の通り、拡張カルマンフィルタの枠組みにおいて、現在位置と誤差情報505を逐次更新することで、現在位置を高精度に推定することができる。慣性装置1に搭載された外部のアプリケーションは、現在位置と誤差情報505を参照することにより、現在位置の推定位置だけでなく、慣性装置の向き(ヘディング情報、ヨー角推定値)を取得することができる。
なお、第一の観測更新処理部502又は第二の観測更新処理部503が慣性装置1の現在位置と誤差情報505を算出又は更新すると、絶対位置補正フラグ319がTRUEとなる。TRUEとなった絶対位置補正フラグ319は補正パラメータ推定処理部700に通知され、これによりFALSEとなる(次ステップの絶対位置・誤差情報が通知される前にクリアされる)。
<<マップマッチング>>
マップマッチング部600は、マップマッチング処理部601を含み、最新の現在位置と誤差情報505と、進行速度推定ベクトル316を取得し、マップマッチング処理を実行する。
マップマッチング部600は、マップマッチング処理部601を含み、最新の現在位置と誤差情報505と、進行速度推定ベクトル316を取得し、マップマッチング処理を実行する。
マップマッチング処理部601は、最新の現在位置と誤差情報505と、進行速度推定ベクトル316を取得し、予め用意された地図データベース602を参照して、進行可能な領域を表す領域情報を取得する。次に、マップマッチング処理部601は、既知の技術であるパーティクルフィルタのアルゴリズムを適用し(例えば、非特許文献3)、マップマッチング処理を実行する。マップマッチング処理部601は、現在位置が、マップ上で進行が不可能な場所を示していた場合には、進行可能な領域に現在位置を修正する。また、マップマッチング処理部601は、内部処理により、修正された位置を表す位置ベクトルの各成分の誤差量を表す誤差共分散行列σmを算出する。マップマッチング処理部601は、修正された現在位置を表す位置ベクトル(Xm、Ym、Zm)及び誤差共分散行列σmを、上記したように、絶対位置推定部500の第三の観測更新処理部504に渡す。
<速度補正部>
速度補正部350は、補正パラメータ推定処理部700より通知される歩行速度推定パラメータに基づいて進行速度推定ベクトル316を補正する処理を実行する。補正の方法は、例えば、補正速度推定パラメータを進行速度推定ベクトル316に乗算することなどでよい。あるいは、ルックアップテーブルに補正速度推定パラメータと補正値を対応付けておき、補正速度推定パラメータに対応付けられた補正値をルックアップテーブルから取得して進行速度推定ベクトル316に乗算してもよい。この場合、非線形な補正が可能になる。なお、補正方法はこれらに限定はされない。
図33は、補正後の速度情報を模式的に示す。波形合成後の速度情報が補正速度推定パラメータにより補正され、補正後の速度情報316Aが得られる。なお、図33では速度情報が大きくなるように補正されているが、速度情報が小さく補正される場合もある。
<測位誤差・座標推定部>
測位誤差・座標推定部330は、絶対位置座標に対する誤差情報を推定する処理を実行する。例えばBluetooth(登録商標)通信モジュール20などの電波で絶対位置情報を慣性装置1が受信することができるが、絶対位置情報は絶対位置情報を送信する送信機340の位置を示しており、慣性装置1の位置を示していない。しかし、この場合、受信信号強度の大きさから誤差を推定することができる。
表1は図1で受信エリアAREA1と受信エリアAREA2を形成する送信機340が送信した絶対位置情報を示す。絶対位置情報IDに対応付けて緯度と経度が周囲の慣性装置1に送信されている。慣性装置1はこの絶対位置情報を受信した際の受信信号強度を測定している。測位誤差・座標推定部330はこの受信信号強度を用いることで誤差情報を推定できる。
図34は絶対位置情報の送信機340と慣性装置1との距離に応じた受信信号強度の変化の一例を示す図である。このような波形が既知であれば、測位誤差・座標推定部330は慣性装置1が絶対位置情報を受信した受信信号強度により送信機340までの距離を推定できる。具体的には、受信信号強度を距離に変換する関数やテーブルを設置者等が作っておき、測位誤差・座標推定部330はこれらから得られる距離を誤差情報として推定する。このようにして得られた誤差情報は図22の誤差情報406として使用される。
なお、送信機340が音波(音声や非可聴音)を用いて絶対位置情報を送信する場合、音声や非可聴音が到達するまでの時間や受信音量を元に誤差を推定することもできる。
また、送信機340は絶対位置情報そのものでなく絶対位置情報IDだけを送信する場合がある。この場合、慣性装置1は予め用意された絶対位置情報IDと絶対位置情報の対応テーブルを参照して絶対位置情報を取得する。
<補正パラメータ推定処理部700>
図4に示した補正パラメータ算出部701は、補正パラメータ算出部701、移動状態検出部702、及び、測位履歴記憶部703を有している。まず、測位履歴記憶部703には、位置推定部800より通知される絶対位置と誤差情報が時刻情報に対応付けて記憶される。なお、測位履歴記憶部703はRAM12やROM13に構築される。
測位履歴記憶部703には、パラメトリックな手法で第一の観測更新処理部502又は第二の観測更新処理部503により補正された絶対位置を記憶する。また、時間発展処理部501が慣性航法で取得した位置も記憶する。下記の移動状態の判定を行わない場合、補正された絶対位置だけを記憶してもよい。
移動状態検出部702は、慣性装置1によって検出される移動が人の歩行や走行によるものか、又は、絶対位置補正によるものかを検出する。具体的には、絶対位置補正フラグのTRUEを検出した時点から、次に絶対位置補正フラグがTRUEとなるまでの間の絶対位置を測位履歴記憶部703に記憶しておく。絶対位置補正フラグがTRUEとなってから次に絶対位置補正フラグがTRUEとなるまでに検出される絶対位置の変化は慣性航法のみにより得られた移動距離であるので、これが所定距離以上であれば、歩行又は走行による移動と判定することができる。この所定距離は適宜決定されるものとする。例えば、数メートルから数十メートルしてよいがこれに限られるものではない。
図1を用いて補正パラメータ算出部701の処理を説明する。最後に絶対位置補正フラグがTRUEとなった際に絶対位置情報を用いてパラメトリックな手法で補正された慣性装置1の絶対位置がP´1である。軌跡1012を構成する絶対位置、位置1001、及びP´1は測位履歴記憶部703に記憶されている。また、絶対位置情報を用いてパラメトリックな手法で慣性装置1の絶対位置が補正されるとP´2も記憶される。P´1、P´2には絶対位置補正フラグ(TRUE)が付与されている。
次に、絶対位置補正フラグがTRUEとなる前に、慣性装置1は慣性航法のみにより推定された位置1001にいる。一方、位置1001は受信エリアAREA2で受信される絶対位置情報で位置P´2に補正される。
従って、受信エリアAREA1から受信エリアAREA2に慣性装置1が移動する過程で距離と方位にずれが生じたことになる。補正パラメータ算出部701はこのずれを補正する補正パラメータを作成する。具体的には、P'1と位置1001の距離Bと、P´1とP´2の距離Aとを算出する。慣性航法が算出する位置は「距離A/距離B」だけ大きい傾向にあるので、「距離A/距離B」を歩行速度推定パラメータとする。
また、P'1とP'2を結ぶ方向と、P'1と位置1001を結ぶ方向の差が、慣性装置1の正しい方向に対する誤差なので、この差をヨー角補正値Δψとすることができる。P'1と、P'2又は位置1001を結ぶ方向でなくても、例えばP'1から歩行者が若干歩行した後の位置とP'2又は位置1001を結ぶ方向を用いてもヨー角補正値Δψを得ることができる。従って、方位の算出の起点となるのはP'1そのものでなくてもよく、測位履歴記憶部703に記憶されている過去の絶対位置でもよい場合がある。ただし、絶対位置補正フラグ(TRUE)が付与された後の絶対位置のうち絶対位置補正フラグ(TRUE)が付与された絶対位置に近いものほど好ましい。
上記した方向の差をそのままヨー角補正値Δψとするのでなく、補正パラメータ算出部701はパラメトリックな手法で位置が補正された際に算出される誤差情報の大きさに応じて方向の差を修正してもよい。例えば、誤差情報を予め定めた閾値で除算して、結果が1以下であればそのまま方向の差をそのままヨー角補正値Δψに決定する。除算した結果が1より大きい場合、除算した結果を方向の差に乗じた値をヨー角補正値Δψに決定する。これにより、誤差情報が大きい場合に進行方向を大きく補正することができる。
上記のように、姿勢情報変換処理部100はヨー角補正値で慣性装置1の方向(ヨー角)を補正するので慣性装置1の方向は絶対位置情報を取得する毎に正しい方向に補正される。また、位置推定部800は歩行速度推定パラメータで進行速度推定ベクトルを補正するので、慣性航法により検出される絶対位置を適切に補正できる。
<慣性装置の全体的な動作>
図35は、歩行速度推定パラメータとヨー角補正値の推定手順を示すフローチャート図の一例である。
位置推定部800は、進行方向推定部200が推定した進行方向213を用いて慣性航法による人の歩行速度と現在位置と推定する(S10)。
次に、慣性装置1が受信エリアに進入することで、絶対位置情報入力部400が送信機340から絶対位置情報を取得すると(S20のYes)、位置推定部800は現在位置をパラメトリックな手法で補正し、位置座標と該位置座標に対応する誤差情報を推定する(S30)。
次に、位置推定部800は、絶対位置情報を用いて補正された位置座標と誤差情報を、慣性航法による現在の位置に設定する(S40)。すなわち、慣性航法が歩行速度に基づいて位置座標を累積するための初期位置が、絶対位置情報を用いて補正された位置座標と誤差情報にリセットされる。なお、TRUEの絶対位置補正フラグが補正パラメータ推定処理部700に送信される。
次に、補正パラメータ算出部701は、補正条件が満たされるか否かを判定する(S50)。補正条件は、慣性航法のみによる位置座標の履歴が所定距離以上の任意の2点間を歩行したことを示したこと、かつ、絶対位置情報を用いたパラメトリックな手法で得られた位置座標の補正後の誤差情報がある閾値以下であることである。前者により絶対位置情報の受信により位置が補正されたものでないと判定でき(慣性装置1によって検出される移動が人の歩行や走行によるものであると判定でき)、後者により慣性航法による測位結果又は絶対位置情報になんらかのエラーがあると判定できる。なお、前者と後者の両方を判定することで好ましい状況で補正パラメータを作成できるが、前者と後者の一方のうちどちらのみを判定してもよいし、どちらも判定しなくてもよい。
ステップS50の判定がYesの場合、補正パラメータ算出部701は、今回、補正された最新の位置座標と前回の補正後の位置座標との2点間の座標差(距離)と方向を算出する(S60)。
次に、補正パラメータ算出部701は慣性航法のみで得られた座標差(距離)と2点間の方向を算出する(S70)。慣性航法のみで得られた2点の座標のうちの始点は上記のようにリセットされた値である。
そして、補正パラメータ算出部701は、ステップS60の座標差をS70の座標差で除算した結果を、歩行速度補正パラメータに決定し位置推定部800の速度補正部350に通知する(S80)。
また、補正パラメータ算出部701は、ステップS60の方向とS70の方向の差分をヨー角補正値として姿勢情報変換処理部100へフィードバックする(S90)。以降、ステップS10〜S90が繰り返し実行される。
<測位結果>
図36は、5箇所の受信エリアを移動した際に得られた測位結果と補正の様子を示す図の一例である。P1〜P5は送信機340が送信する絶対位置情報であり、P´1〜P´5は絶対位置情報を用いてパラメトリックな手法で補正された位置である。また、位置1002〜1005は慣性航法のみにより推定された位置である。AREA1〜5は送信機340からの絶対位置情報を受信可能な範囲である。E1〜E5はパラメトリックな手法で位置が補正された際に算出された誤差情報である。
受信エリアAREA2では慣性航法により測位された位置1002と位置P2とが大きくずれているが、絶対位置情報により位置P´2に補正されている。その後、受信エリアAREA3〜5では慣性航法により測位された位置1003〜1005と位置P´3〜5とが大きくずれていないが、各受信エリアで位置1002〜1005が補正されていることが分かる。
表2は、測定機が出力する絶対位置情報と、パラメトリックな手法で補正された位置とを数値で示す。
図36や表2から分かるように、パラメトリックな手法を適用することで、絶対位置情報を受信した際の誤差情報が加味されて絶対位置が更新される。このため、誤差が限りなく小さい場合にはP1〜P5の絶対位置情報が示す位置座標の値に補正される位置が近くなり、誤差が大きい場合には位置座標を中心として含まれる誤差分だけ外側にオフセットした位置に現在位置が更新される。
これはパラメトリックな手法を適用することで、誤差の大きさを元に推定値を収束させる信頼度(誤差共分散行列)が決定され、補正前後で位置座標の補間がなされるためである。よって、パラメトリックな手法を適用することによって、絶対位置情報を受信した場合に受信した位置座標に一致させる補正する手法に比べ、より正確な位置推定を行うことができる。
仮に、絶対位置情報を受信した場合に受信した位置座標そのものに一致させる補正手法を元に図35のステップS50〜S90が実行された場合、絶対位置情報の誤差が大きい場合でも該ステップが実行されてしまう。このため、歩行速度補正パラメータや、ヨー角補正値の推定精度が低下する。あるいは誤った推定を行ってしまう。
本実施形態では、パラメトリックな手法を適用した上でステップS50〜S90を実行することで、より精度の高い歩行速度補正パラメータや、ヨー角補正値の推定が可能である。
また、本実施形態では2点の座標から距離を求めるので、2点間を歩行者がどのように歩行したかによって歩行速度補正パラメータが影響を受けない。特許文献2には、ユーザの絶対位置を取得し所定距離移動したことをトリガとして、歩行テンポと、ユーザの歩幅又は移動速度を示すパラメータを対応させて学習する方法が開示されている。しかしながら、この方法では、絶対位置の2点間を蛇行して繋ぐように歩行した場合、歩数カウント数と歩行テンポを元にしていることから距離関係が崩れ、正しい学習結果を得ることはできない。
以上から、歩行速度補正パラメータやヨー角補正値を推定しフィードバックすることによる補正と、パラメトリックな手法を適用することによる歩行速度補正パラメータやヨー角補正値の推定精度の向上の双方によって、慣性航法における位置の推定精度を高めることができる。
図37は、従来技術と本実施形態の技術とで得られた歩行者の移動軌跡を比較する図の一例である。図37(a)は従来技術による移動軌跡を示し、図37(b)は本実施形態の技術による移動軌跡を示している。図37(a)では、方位角のずれと歩行速度の推定誤差が、慣性航法が推定する現在位置に重畳し位置の精度が低下している。図37(b)では、絶対位置情報により歩行速度パラメータとヨー角が補正されるので、位置の精度が向上している。
従来技術では、現在位置を絶対位置に一致させる補正が行われるが、ヨー角の誤差と歩行速度の推定誤差が慣性航法の際に重畳し位置精度が低下するため、測位精度を継続的に維持することが困難であった。本実施形態の技術では、絶対位置情報を観測したタイミングでヨー角と歩行速度推定パラメータが同時に更新されるため、測位精度を向上させることが可能になる。
<その他の適用例>
以上、本発明を実施するための最良の形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。
例えば、本実施形態ではパラメトリックな手法としてカルマンフィルタを用いて位置推定部800が慣性航法により位置を推定し、また、この位置を絶対位置情報で補正した。しかし、パラメトリックな手法として、パーティクルフィルタ、α−βフィルター等の推定手法を用いてもよい。また、これらが変形又は拡張されたフィルターが使用されてもよい。また、LMS(least-mean-square(LMS))や最小最急勾配法などで位置を推定してもよい。
また、図4などの構成例は、慣性装置1による処理の理解を容易にするために、主な機能に応じて分割したものである。処理単位の分割の仕方や名称によって本願発明が制限されることはない。慣性装置1の処理は、処理内容に応じて更に多くの処理単位に分割することもできる。また、1つの処理単位が更に多くの処理を含むように分割することもできる。
また、本実施形態では慣性装置1が全ての処理を行ったが、慣性装置1が加速度、角速度及び地磁気を検出してそれらをサーバに送信してもよい。この場合、サーバは慣性装置1の位置を推定して慣性装置1に送信するので、慣性装置1の処理負荷を低減できる。
なお、図1の位置座標(P´1)は所定位置の一例であり、第一の観測更新処理部502又は第二の観測更新処理部503は位置補正部の一例であり、補正パラメータ推定処理部700は補正パラメータ算出部の一例であり、歩行速度補正パラメータは移動速度補正パラメータの一例であり、座標系変換処理部106は進行方向補正部の一例である。測位履歴記憶部703は記憶装置の一例であり、移動状態検出部702は記憶部の一例である。距離Aは第1距離の一例であり、距離Bは第2距離の一例である。速度補正部350は移動速度補正部の一例であり、測位履歴記憶部703は位置記録部の一例である。絶対位置情報取得部は絶対位置情報入力部400の一例である。慣性装置1は自装置の一例であり、位置1001は推定位置の一例であり、P´2は補正位置の一例である。