JP2013200162A - Compact attitude sensor - Google Patents

Compact attitude sensor Download PDF

Info

Publication number
JP2013200162A
JP2013200162A JP2012067533A JP2012067533A JP2013200162A JP 2013200162 A JP2013200162 A JP 2013200162A JP 2012067533 A JP2012067533 A JP 2012067533A JP 2012067533 A JP2012067533 A JP 2012067533A JP 2013200162 A JP2013200162 A JP 2013200162A
Authority
JP
Japan
Prior art keywords
sensor
acceleration
error
angular velocity
quaternion
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.)
Granted
Application number
JP2012067533A
Other languages
Japanese (ja)
Other versions
JP5061264B1 (en
Inventor
Kenzo Nonami
健蔵 野波
Satoshi Suzuki
智 鈴木
Makoto Tawara
誠 田原
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chiba University NUC
Hirobo Ltd
Original Assignee
Chiba University NUC
Hirobo Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Chiba University NUC, Hirobo Ltd filed Critical Chiba University NUC
Priority to JP2012067533A priority Critical patent/JP5061264B1/en
Application granted granted Critical
Publication of JP5061264B1 publication Critical patent/JP5061264B1/en
Publication of JP2013200162A publication Critical patent/JP2013200162A/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Landscapes

  • Gyroscopes (AREA)
  • Navigation (AREA)

Abstract

PROBLEM TO BE SOLVED: To provide a compact attitude sensor capable of estimating attitude and orientation of a mobile object under acceleration with high precision using an algorithm that takes the gyro bias error into account.SOLUTION: A compact attitude sensor includes an angular velocity sensor 101 for measuring angular velocity of a mobile object on three independent axes, an acceleration sensor 102, a magnetic sensor 103, and an arithmetic processing unit 104 which computes a quaternion estimate based on the values measured by the angular velocity sensor 101, the acceleration sensor 102, and the magnetic sensor 103. The arithmetic processing unit 104 computes a bias error, a scale factor error, a misalignment error, an angular velocity sensor error including a noise factor, and a quaternion level error for angular velocity data, and obtains a quaternion estimate representing the current attitude of the mobile object by using an extended Kalman filter obtained by using the error values. The compact attitude sensor outputs the quaternion estimate.

Description

本発明は姿勢制御を要する移動物体に搭載され、ジャイロバイアス誤差を導入したアルゴリズムを用いて、加速度環境下における高精度の姿勢・方位推定を実現する小型姿勢センサに関する。   The present invention relates to a small posture sensor that is mounted on a moving object that requires posture control and realizes high-precision posture / orientation estimation under an acceleration environment using an algorithm that introduces a gyro bias error.

小型無人ヘリコプタなどの姿勢制御を要する移動物体には、その三次元姿勢角を計測するため姿勢センサを備えている。一般的に姿勢センサは移動物体の独立な3軸上の角速度を計測する角速度センサ、及び角速度センサにより得られる角速度の計測値を時間的に積分することにより、移動物体の三次元姿勢角を取得する処理を行う演算処理部を備えている。三次元姿勢角はロール・ピッチ・ヨー角等を用いてオイラー角により表現されるものがあり、姿勢センサにより計測される三次元姿勢角を基に移動物体の姿勢制御が行われる。このような移動物体の姿勢制御技術の従来例としては、例えば特許文献1に記載されたものがある。また、本発明の基礎となる技術として特許文献2に記載された発明がある。   A moving object that requires posture control, such as a small unmanned helicopter, is equipped with a posture sensor to measure its three-dimensional posture angle. In general, the posture sensor acquires the three-dimensional posture angle of a moving object by temporally integrating the angular velocity sensor that measures the angular velocity on three independent axes of the moving object and the angular velocity measurement value obtained by the angular velocity sensor. An arithmetic processing unit for performing the processing is provided. Some three-dimensional posture angles are expressed by Euler angles using roll, pitch, yaw angle, etc., and the posture control of the moving object is performed based on the three-dimensional posture angle measured by the posture sensor. As a conventional example of such a moving object posture control technique, for example, there is one described in Patent Document 1. Moreover, there exists invention described in patent document 2 as a technique used as the foundation of this invention.

特開平9−5104号公報Japanese Patent Laid-Open No. 9-5104 特開2007−183138号公報JP 2007-183138 A

以下、従来の移動物体の姿勢制御技術について説明する。本発明の前提となるセンサとして特許文献2に示されたセンサを作製し、これをMARGセンサという。このMARGセンサの仕様を表1に示す。
Hereinafter, a conventional attitude control technique for a moving object will be described. A sensor disclosed in Patent Document 2 is manufactured as a sensor that is a premise of the present invention, and this is referred to as a MARG sensor. Table 1 shows the specifications of the MARG sensor.

表1に示されるような比較的低い精度ではあるが、非常に軽量で小型無人ヘリコプタにも充分搭載可能なハードウェアに対して提案アルゴリズムを実装し、検証を行う。このような小型軽量のセンサが実現されれば、重量の上では小型無人ヘリコプタに限らず一般的な航空機や移動ロボット等への応用も期待できる。以下、この説明で用いる座標系及びクォータニオンについて説明を行う。   Although the accuracy is relatively low as shown in Table 1, the proposed algorithm is implemented and verified for hardware that is very lightweight and can be fully installed in a small unmanned helicopter. If such a small and lightweight sensor is realized, it can be expected to be applied not only to a small unmanned helicopter but also to a general aircraft or a mobile robot in terms of weight. Hereinafter, the coordinate system and quaternion used in this description will be described.

1)座標系及びクォータニオン
先ず、この説明で用いる座標系および各座標系上で表されるベクトルの表記方法についての定義を行う。ここで用いる座標系を図19に示す。ここで、図中のr-frame はReference frame (参照フレーム)といい、地上の任意の点を原点として、磁北をXr軸、重力方向をZr軸、XrZr平面の垂直方向をYr軸にとった座標系である。また、b-frame はBody frame (機体フレーム)といい、機体の重心を原点として、機体前方をXb軸、機体下方をZb軸、XbZb平面の垂直方向をYb軸にとった座標系である。ここで、三次元空間中の任意の幾何ベクトルをrとしたとき、rを各座標系上の代数ベクトルとして表したものをそれぞれrr,rbと定義する。また、r-frameに対する b-frameの姿勢を機体姿勢として定義する。
1) Coordinate system and quaternion First, the coordinate system used in this description and the notation method of vectors represented on each coordinate system are defined. The coordinate system used here is shown in FIG. Here, the r-frame in the figure is referred to as a reference frame, with an arbitrary point on the ground as the origin, magnetic north as the Xr axis, gravity direction as the Zr axis, and XrZr plane vertical direction as the Yr axis. Coordinate system. The b-frame is called a body frame, and is a coordinate system in which the center of gravity of the aircraft is the origin, the Xb axis is the front of the aircraft, the Zb axis is below the aircraft, and the Yb axis is the vertical direction of the XbZb plane. Here, assuming that an arbitrary geometric vector in the three-dimensional space is r, those expressing r as an algebraic vector on each coordinate system are defined as r r and r b , respectively. Also, the attitude of the b-frame relative to the r-frame is defined as the body attitude.

次に、この説明で用いる機体姿勢の表現方法として用いられるクォータニオンの概要を説明する。クォータニオンとは、複素数(通常、a+biのように虚数部分が一次元で表現される)を三次元に拡張した数であり、次式で示すようにクォータニオンqは1つの実数部と3つの虚数部から構成されている。
q=qo +q1i+q2j+q3k ・・・(1)
ただしここで、式中のi,j,kは虚数単位であり、それぞれ次式のような関係を満たしている。
2=j2=k2=−1 ・・・(2)
ij=k, jk=i, ki=j ・・・(3)
Next, an outline of a quaternion used as a method for expressing the body posture used in this description will be described. A quaternion is a number obtained by expanding a complex number (usually the imaginary part is represented in one dimension as in a + bi) to three dimensions, and the quaternion q has one real part and three imaginary parts as shown in the following equation. It is composed of
q = q o + q 1 i + q 2 j + q 3 k (1)
However, here, i, j, and k in the formula are imaginary units, and satisfy the relationship as shown in the following formula.
i 2 = j 2 = k 2 = −1 (2)
ij = k, jk = i, ki = j (3)

クォータニオンは四元数、又はオイラーパラメータと呼ばれることもあり(例えば文献、田島洋、マルチボディダイナミクス−3次元運動方程式の立て方、東京電機大学出版局、2006年 を参照)、工学的には物体の回転姿勢の表現方法として用いられることが多い。回転姿勢を表現するクォータニオンは上記文献のSimple Rotationを用いて定義することができる。いま、オイラーの定理によれば、図10中のb-frame がr-frameに対してどのような回転姿勢をとっていた場合でも、回転軸を表す単位ベクトルとその軸周りの回転角度を用いてその回転姿勢を表現することが可能である。上記ベクトルと回転角度の組を用いた回転姿勢の表現法をSimple Rotationと呼ぶ。ただしここで、ベクトルは式4のように成分で表すことができ、そのノルムは(5)式に示すように1となる。
V=[l m n]T
・・・(4)
ここで、V:ベクトル
(l2+m2+n21/2=1
・・・(5)
Quaternions are sometimes referred to as quaternions or Euler parameters (see, for example, Hiroshi Tajima, Multibody Dynamics-3D Equations of Motion, Tokyo Denki University Press, 2006). It is often used as a method for expressing the rotational posture of The quaternion that expresses the rotation posture can be defined using Simple Rotation in the above document. Now, according to Euler's theorem, the unit vector representing the rotation axis and the rotation angle around that axis are used regardless of the rotation posture of the b-frame in FIG. 10 relative to the r-frame. It is possible to express the rotation posture. A method of expressing a rotation posture using a set of the vector and the rotation angle is called Simple Rotation. However, the vector can be expressed by a component as shown in Equation 4, and its norm is 1 as shown in Equation (5).
V = [l m n] T
... (4)
Where V: vector (l 2 + m 2 + n 2 ) 1/2 = 1
... (5)

以上のようなSimple Rotationを用いてr-frameに対する b-frameの回転姿勢を表現するクォータニオンは次式のように定義される。
ここで、クォータニオンのノルムは式(5)より常に1となる。
The quaternion that expresses the rotation orientation of the b-frame with respect to the r-frame using Simple Rotation as described above is defined as follows.
Here, the norm of the quaternion is always 1 from the equation (5).

クォータニオンはq=[qo q123Tというように虚数単位を省いたベクトル表記で表されることも多く、この表記法によってクォータニオンp,qの和、差、積は次式のように定義される。
q+p=[qo+po q1+p12+p23+p3T
・・・(7)
q−p=[qo−po q1−p12−p23−p3T
・・・(8)
・・・(9)
Quaternions are often expressed in vector notation such that q = [qo q 1 q 2 q 3 ] T , and the sum, difference, and product of quaternions p and q are expressed as Is defined as
q + p = [qo + po q 1 + p 1 q 2 + p 2 q 3 + p 3 ] T
... (7)
q-p = [qo-po q 1 -p 1 q 2 -p 2 q 3 -p 3] T
... (8)
... (9)

また、qに共役なクォータニオンq*は次式のように定義される。
*=[qo −q1 −q2 −q3T
・・・(10)
いま、r-frame上の代数ベクトルrrからb-frame上の代数ベクトルrbへの座標変換は次式で表すことができる。
・・・(11)
ただしここで、式中の(rbqおよび(rrqは三次元ベクトルをクォータニオン表記したものであり、それぞれ次式のように定義される。
・・・(12)
A quaternion q * conjugated to q is defined as follows.
q * = [qo−q 1 −q 2 −q 3 ] T
... (10)
Now, the coordinate transformation from the algebraic vector r r on r-frame to the algebraic vector r b on b-frame can be expressed by the following equation.
(11)
However, here, (r b ) q and (r r ) q in the equation are quaternion notation of a three-dimensional vector, and are respectively defined as the following equations.
(12)

[MARGセンサ]
ここでは、先行技術において開発されたMARGセンサのハードウェアの諸元と従来の機体姿勢推定アルゴリズムの概要およびその問題点について説明する。MARGセンサの諸元を表2に示す。
[MARG sensor]
Here, the hardware specifications of the MARG sensor developed in the prior art, the outline of the conventional aircraft attitude estimation algorithm, and its problems will be described. Table 2 shows the specifications of the MARG sensor.

このセンサは、3軸の加速度センサ、ジャイロセンサ、磁気センサの出力を用いて現在の姿勢を表すクォータニオンを推定するアルゴリズムを実装した姿勢センサの一種である。   This sensor is a kind of attitude sensor that implements an algorithm for estimating a quaternion that represents the current attitude by using outputs of a three-axis acceleration sensor, a gyro sensor, and a magnetic sensor.

従来のMARGセンサに実装されていた姿勢推定アルゴリズムの概要について述べると、先ず、r-frame上での重力ベクトルをgr =[0 0 g]T,地磁気ベクトルをmr =[mn 0 mdT とする。ただし、gは重力加速度、mnは水平地磁気、mdは垂直地磁気をそれぞれ表している。ここで、現在の姿勢を表すクォータニオンの推定値からb-frame上での重力ベクトルおよび地磁気ベクトル推定値を求めると次式となる。
The outline of the posture estimation algorithm implemented in the conventional MARG sensor will be described. First, the gravity vector on the r-frame is g r = [0 0 g] T , and the geomagnetic vector is m r = [m n 0 m d ] T. However, g is the gravitational acceleration, m n denotes the horizontal geomagnetic, m d is the vertical geomagnetic respectively. Here, when the gravity vector and the geomagnetic vector estimation value on the b-frame are obtained from the quaternion estimation value representing the current posture, the following equations are obtained.

続いて、機体に搭載された加速度センサおよび磁気センサの出力をそれぞれameasure=[axyzT、mmeasure=[mxyzTとすると、これらと上述したb-frame上での重力ベクトルおよび地磁気ベクトル推定値との誤差ベクトルは次式となる。
Subsequently, assuming that the outputs of the acceleration sensor and the magnetic sensor mounted on the body are a measure = [a x a y a z ] T and m measure = [m x m y m z ] T , respectively, and b described above The error vector between the gravity vector and the geomagnetic vector estimate on the -frame is as follows.

いま、各センサの誤差がないと仮定すれば、最も確からしいクォータニオン推定値は式(15)において、
とした同次連立方程式を解くことで得られることが分かる。しかし、式(15)は4個の未知数に対して式が6個ある過剰条件の連立方程式であるため、一般には解くことができない。よって、最小二乗法を用いて誤差の二乗和を最小とするような近似解を求める。Gauss−Newton反復法を用いれば、あるステップにおいて誤差の二乗和を最小に近づけるような修正ベクトルは次式となる。
ただしここでXは以下に示すようなヤコビアンである。
Assuming that there is no error in each sensor, the most probable quaternion estimate is
It can be seen that it can be obtained by solving the simultaneous equations. However, since equation (15) is an over-condition simultaneous equation with six equations for four unknowns, it cannot generally be solved. Therefore, an approximate solution that minimizes the sum of squared errors is obtained using the least square method. If the Gauss-Newton iteration method is used, a correction vector that brings the sum of squares of errors close to the minimum in a certain step is represented by the following equation.
Where X is a Jacobian as shown below.

一方、現在の姿勢を示す真のクォータニオンの時間微分と機体角速度ωbとは以下の関係式であることが知られている。
これにより、クォータニオンの推定値の時間微分に関しても同様に次式のような関係が成り立つことが分かる。
On the other hand, it is known that the true quaternion time derivative indicating the current posture and the airframe angular velocity ω b are the following relational expressions.
Accordingly, it can be understood that the following relationship is similarly established for the time differentiation of the estimated value of the quaternion.

式(16)と式(20)を次のように組み合わせることによって、クォータニオン推定値を得るというのが従来のアルゴリズムである。ただしここで、kはスカラゲインである。
The conventional algorithm is to obtain the quaternion estimated value by combining the equations (16) and (20) as follows. Here, k is a scalar gain.

以上のアルゴリズムを信号線で表すと図20のようになる。このアルゴリズムは単純なフィルタアルゴリズムであり、角速度から計算されたクォータニオンの時間微分を加速度、地磁気のデータを用いて補正し、それを積分した後に正規化して単位クォータニオンの推定値を求めるというものである。基本的に、低周波数の姿勢変動は加速度、地磁気のデータから推定する一方、高周波数の姿勢変動はジャイロを用いて推定するという構成となっている。しかし、このアルゴリズムには以下のような問題点がある。   The above algorithm is represented by signal lines as shown in FIG. This algorithm is a simple filter algorithm, which corrects the quaternion time derivative calculated from the angular velocity using acceleration and geomagnetic data, integrates it, and normalizes it to obtain an estimated value of unit quaternion. . Basically, low-frequency posture fluctuations are estimated from acceleration and geomagnetism data, while high-frequency posture fluctuations are estimated using a gyro. However, this algorithm has the following problems.

式(15)において、amには重力をb−frame上に座標変換した成分だけではなく、センサが動くことで発生する動的加速度の成分も含まれている。このため、動的加速度が印加されることで姿勢データが大きな誤差を持ってしまうのである。実際にセンサを水平に保った状態で大きな加速度を与えた際の姿勢データを図21に示す。ただし、図21に示されているのは推定されたクォータニオンを変換して得られたオイラー角である。図21より、加速度を加えることで最大で20[deg]程度の姿勢誤差が生じていることが分かる。このような動的加速度による誤差は移動体に搭載するためのセンサとしては致命的であり、なんとしても誤差を低減化する必要がある。   In Expression (15), am includes not only a component obtained by coordinate conversion of gravity on b-frame but also a component of dynamic acceleration generated by the movement of the sensor. For this reason, the posture data has a large error by applying the dynamic acceleration. FIG. 21 shows posture data when a large acceleration is applied while the sensor is actually kept horizontal. However, what is shown in FIG. 21 is the Euler angle obtained by converting the estimated quaternion. From FIG. 21, it can be seen that a posture error of about 20 [deg] at maximum is generated by applying acceleration. Such an error due to dynamic acceleration is fatal for a sensor to be mounted on a moving body, and it is necessary to reduce the error.

本発明はこのような従来技術の不具合に鑑みてなされたもので、その目的は、拡張カルマンフィルタのアルゴリズムを用いて、加速度環境下、とりわけ低周波加速度外乱の下における高精度の姿勢・方位推定を実現する小型姿勢センサを提供することである。   The present invention has been made in view of the above-described problems of the prior art. The purpose of the present invention is to perform highly accurate posture / orientation estimation under an acceleration environment, particularly under a low-frequency acceleration disturbance, using an extended Kalman filter algorithm. It is to provide a small attitude sensor to be realized.

本発明は上記目的を達成するために、移動物体の独立な3軸上の角速度を計測する角速度センサ、加速度センサ及び磁気センサと、前記角速度センサ、前記加速度センサ、及び前記磁気センサの計測値を基にクォータニオンの推定値を算出する演算処理部と、を備え、前記演算処理部は、予め推定されたクォータニオン推定値を用いて座標変換行列を算出し、当該座標変換行列を用いて地上固定座標系成分の地磁気ベクトル及び重力加速度ベクトルを変換し、現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルを推定する座標変換部と、前記磁気センサ及び前記加速度センサから現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルを取得する磁気・加速度センサデータ取得部と、前記座標変換部で推定した現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルと、前記磁気・加速度センサデータ取得部で取得した現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルとの誤差を算出するベクトル誤差算出部と、前記角速度センサから現ステップのセンサ座標系成分の角速度ベクトルを取得する角速度センサデータ取得部と、前記取得した角速度データについて、バイアス誤差、スケールファクタ誤差、ミスアラインメント誤差、ノイズの要素を含む角速度センサ誤差を演算により求める角速度データ誤差算出部と、前ステップのクォータニオンと現ステップのクォータニオンの時間変化を求めるクォータニオン算出部と前記座標変換部、前記磁気・加速度センサデータ取得部、前記ベクトル誤差算出部、前記角速度データ誤差算出部、及び前記クォータニオン算出部を基に現在の姿勢を表すクォータニオンおよび角速度誤差を推定する拡張カルマンフィルタを構成し、前記拡張カルマンフィルタによるクォータニオン推定値を現在の姿勢として出力する小型姿勢センサを提供する。 In order to achieve the above object, the present invention provides an angular velocity sensor, an acceleration sensor, and a magnetic sensor that measure angular velocities on three independent axes of a moving object, and measured values of the angular velocity sensor, the acceleration sensor, and the magnetic sensor. An arithmetic processing unit that calculates a quaternion estimated value based on the quaternion estimated value, the coordinate processing matrix is calculated using the estimated quaternion estimated value, and the ground fixed coordinates are calculated using the coordinate transformation matrix. A coordinate conversion unit for converting a geomagnetic vector and a gravitational acceleration vector of a system component to estimate a geomagnetic vector and a gravitational acceleration vector of a sensor coordinate system component of the current step; and a sensor coordinate system component of the current step from the magnetic sensor and the acceleration sensor A magnetic / acceleration sensor data acquisition unit that acquires a geomagnetic vector and a gravitational acceleration vector of A vector for calculating an error between the geomagnetic vector and gravitational acceleration vector of the sensor coordinate system component of the determined current step and the geomagnetic vector and gravitational acceleration vector of the sensor coordinate system component of the current step acquired by the magnetic / acceleration sensor data acquisition unit. An error calculation unit, an angular velocity sensor data acquisition unit that acquires an angular velocity vector of a sensor coordinate system component of the current step from the angular velocity sensor, and elements of bias error, scale factor error, misalignment error, and noise for the acquired angular velocity data An angular velocity data error calculating unit that calculates an angular velocity sensor error including calculation, a quaternion calculating unit that calculates a time change between the quaternion of the previous step and the quaternion of the current step, the coordinate conversion unit, the magnetic / acceleration sensor data acquisition unit, and the vector Error calculation A quaternion representing the current posture and an extended Kalman filter for estimating the angular velocity error based on the angular velocity data error calculating unit and the quaternion calculating unit, and outputting a quaternion estimated value by the extended Kalman filter as the current posture An attitude sensor is provided.

本発明によれば、角速度データについて、バイアス誤差、スケールファクタ誤差、ミスアラインメント誤差、ノイズの要素を含む角速度センサ誤差を演算により求めてこの誤差を補正することにより、加速度環境下における高精度の姿勢・方位推定を実現する小型姿勢センサを実現できる。   According to the present invention, with respect to angular velocity data, a bias error, a scale factor error, a misalignment error, and an angular velocity sensor error including a noise element are obtained by calculation, and this error is corrected, whereby a highly accurate posture in an acceleration environment is obtained.・ A small-size attitude sensor that realizes direction estimation can be realized.

また、本発明では、移動体に一定値の低周波数の加速度外乱が加えられたときでも、正常に低周波数の加速度による外乱が推定され、低周波加速度環境下における高精度の姿勢・方位推定を実現する小型姿勢センサを実現できる。   In the present invention, even when a constant low frequency acceleration disturbance is applied to the moving body, the disturbance due to the low frequency acceleration is normally estimated, and the posture / orientation estimation with high accuracy in the low frequency acceleration environment is performed. A small attitude sensor can be realized.

本発明の第1の実施の形態で用いる小型姿勢センサに備えられる各構成をブロックで示した図The figure which showed each structure with which the small attitude | position sensor used in the 1st Embodiment of this invention was equipped with the block 上記第1の実施の形態で用いる小型姿勢センサの演算処理部で行う処理をブロックで示した図The figure which showed the process performed by the arithmetic processing part of the small attitude | position sensor used in the said 1st Embodiment with the block 上記第1の実施の形態で用いる小型姿勢センサの演算処理部で行う処理をフローで示した図The figure which showed the process performed by the arithmetic processing part of the small posture sensor used in the said 1st Embodiment with the flow 上記第1の実施の形態において、小型姿勢センサを移動物体に搭載し、動的加速度による誤差が少なくなるように観測ノイズの共分散行列の値を決定したときのシミュレーション結果を示す図The figure which shows the simulation result when the value of the covariance matrix of observation noise is determined so that the small posture sensor may be mounted on a moving object and the error due to dynamic acceleration may be reduced in the first embodiment. 上記第1の実施の形態において、小型姿勢センサを移動物体に搭載し、姿勢の追従を充分に速くなるように観測ノイズの共分散行列の値を決定したときのシミュレーション結果を示す図The figure which shows the simulation result when the value of the covariance matrix of observation noise is determined so that a small posture sensor may be mounted on a moving object and the posture tracking may be made sufficiently fast in the first embodiment. 上記第1の実施の形態において、小型姿勢センサを移動物体に搭載し、動的加速度による誤差を少なく、且つ姿勢の追従を充分に速くできるように観測ノイズの共分散行列の値を決定したときのシミュレーション結果を示す図In the first embodiment, when a small attitude sensor is mounted on a moving object, the value of the covariance matrix of the observation noise is determined so that the error due to dynamic acceleration is small and the attitude can be tracked sufficiently quickly Of simulation results 上記第1の実施の形態において、移動物体に搭載した小型姿勢センサに、線形誤差推定カルマンフィルタ(KF)を用いた場合と、本発明の拡張カルマンフィルタ(EKF)用いた場合におけるロール角についてのシミュレーション結果を比較して示す図In the first embodiment, a simulation result of the roll angle when the linear error estimation Kalman filter (KF) is used for the small attitude sensor mounted on the moving object and when the extended Kalman filter (EKF) of the present invention is used. Figure comparing and showing 上記第1の実施の形態において、移動物体に搭載した小型姿勢センサに、線形誤差推定カルマンフィルタ(KF)を用いた場合と、本発明の拡張カルマンフィルタ(EKF)用いた場合におけるピッチ角についてのシミュレーション結果を比較して示す図In the first embodiment, the simulation results for the pitch angle when using the linear error estimation Kalman filter (KF) and the extended Kalman filter (EKF) of the present invention for the small attitude sensor mounted on the moving object. Figure comparing and showing 上記第1の実施の形態において、移動物体に搭載した小型姿勢センサに、線形誤差推定カルマンフィルタ(KF)を用いた場合と、本発明の拡張カルマンフィルタ(EKF)用いた場合におけるヨー角についてのシミュレーション結果を比較して示す図In the first embodiment, the simulation result of the yaw angle when the linear error estimation Kalman filter (KF) is used for the small posture sensor mounted on the moving object and when the extended Kalman filter (EKF) of the present invention is used. Figure comparing and showing 上記第1の実施の形態において、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るような角速度環境下で、ロール角についてリアルタイムで姿勢・方位推定を行った結果を示す図In the first embodiment, the roll angle is obtained in an angular velocity environment in which the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer and greatly exceed the angular velocity generated when a small unmanned helicopter moves. Showing the result of posture / orientation estimation in real time 上記第1の実施の形態において、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るような角速度環境下で、ピッチ角についてリアルタイムで姿勢・方位推定を行った結果を示す図In the first embodiment, the pitch angle is obtained in an angular velocity environment in which the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer and greatly exceed the angular velocity generated when a small unmanned helicopter moves. Showing the result of posture / orientation estimation in real time 上記第1の実施の形態において、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るような角速度環境下で、ヨー角についてリアルタイムで姿勢・方位推定を行った結果を示す図In the first embodiment, the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer, and the yaw angle is under an angular velocity environment that greatly exceeds the angular velocity generated during the movement of a small unmanned helicopter. Showing the result of posture / orientation estimation in real time 上記第1の実施の形態において、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの飛行中に加わる加速度と比較して充分大きな加速度を印加し、ピッチ角についてリアルタイムで姿勢・方位推定を行った結果を示す図In the first embodiment, the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer, and a sufficiently large acceleration is applied as compared with the acceleration applied during the flight of a small unmanned helicopter. Figure showing the results of posture / orientation estimation in real time 本発明の第2の実施の形態において、移動体に一定値の低周波数の加速度外乱を加えたときの姿勢推定シミュレーションに用いた各軸の加速度データを示す図The figure which shows the acceleration data of each axis | shaft used for the attitude | position estimation simulation when applying the low-frequency acceleration disturbance of a fixed value to a moving body in the 2nd Embodiment of this invention. 上記第2の実施の形態の姿勢推定シミュレーションにおける各アルゴリズムの姿勢推定結果を示す図The figure which shows the attitude | position estimation result of each algorithm in the attitude | position estimation simulation of the said 2nd Embodiment 上記第2の実施の形態の姿勢推定シミュレーションにおける加速度外乱推定値と加速度外乱の真値との比較を示す図The figure which shows the comparison with the acceleration disturbance estimated value in the attitude | position estimation simulation of the said 2nd Embodiment, and the true value of an acceleration disturbance 上記第2の実施の形態において、センサを車に搭載して一定値の低周波数の加速度外乱を加えた実験での各アルゴリズムの姿勢計測結果を示す図The figure which shows the attitude | position measurement result of each algorithm in the experiment which mounted the sensor in the said 2nd Embodiment and added the acceleration disturbance of the low frequency of a fixed value in the said 2nd Embodiment. 上記第2の実施の形態の実験において、本発明のアルゴリズムにおける加速度外乱推定値と加速度外乱の真値との比較を示す図The figure which shows the comparison of the acceleration disturbance estimated value and the true value of an acceleration disturbance in the algorithm of this invention in the experiment of the said 2nd Embodiment. 本発明の前提となるクォータニオンの概要を説明するための三次元座標を示した図The figure which showed the three-dimensional coordinate for demonstrating the outline | summary of the quaternion used as the premise of this invention 従来技術において用いられた小型姿勢センサの演算処理部で処理を行う際に用いる数式やデータの流れを示した図A diagram showing the flow of mathematical formulas and data used when processing is performed by the arithmetic processing unit of the small attitude sensor used in the prior art 従来技術における小型姿勢センサを移動物体に搭載し、クォータニオンにより表現する三次元姿勢角を計測した実験データを示した図Figure showing experimental data obtained by measuring a three-dimensional attitude angle expressed by a quaternion with a conventional attitude sensor mounted on a moving object

以下、本発明を実施するための形態について説明する。本発明においても、上記特許文献2に開示された装置と同様の構成を有するセンサ装置回路が用いられる。ただし、本発明のセンサ装置回路と従来のセンサ装置回路においては演算処理の内容が異なる。   Hereinafter, modes for carrying out the present invention will be described. Also in the present invention, a sensor device circuit having the same configuration as the device disclosed in Patent Document 2 is used. However, the content of arithmetic processing differs between the sensor device circuit of the present invention and the conventional sensor device circuit.

(第1の実施の形態)
図1は本発明の第1の実施の形態に係る小型姿勢センサに備えられる各ハードウェア構成を示すブロック図である。小型姿勢センサ100は、移動物体の独立な3軸上の角速度を計測する角速度センサ101、移動物体の独立な3軸上の加速度を計測する加速度センサ102、移動物体の独立な3軸上の地磁気を計測する磁気センサ103、及び角速度センサ101、加速度センサ102、及び磁気センサ103の計測値を基にクォータニオンを算出する演算処理部104を備えている。この第1の実施の形態の小型姿勢センサ100は1入力1出力の線形システムに適用し、定常カルマンフィルタを構成するとともに、演算項目の中にジャイロセンサ(角速度センサ)の誤差を導入することにより、「カルマンゲイン」や「クォータニオン」の算出をより精密にすることができ、現実的な実装を可能にし、さらに、小型化及び実時間処理を可能とする。
(First embodiment)
FIG. 1 is a block diagram showing each hardware configuration provided in the small posture sensor according to the first embodiment of the present invention. The small posture sensor 100 includes an angular velocity sensor 101 that measures angular velocities on three independent axes of a moving object, an acceleration sensor 102 that measures acceleration on three independent axes of the moving object, and a geomagnetism on three independent axes of the moving object. And an arithmetic processing unit 104 that calculates a quaternion based on measurement values of the angular velocity sensor 101, the acceleration sensor 102, and the magnetic sensor 103. The small posture sensor 100 according to the first embodiment is applied to a linear system with one input and one output, constitutes a stationary Kalman filter, and introduces an error of a gyro sensor (angular velocity sensor) into a calculation item. Calculation of “Kalman gain” and “quaternion” can be made more precise, enabling realistic implementation, and further enabling miniaturization and real-time processing.

図2は演算処理部104で行う処理を機能ブロックで示したものである。演算処理部104は、拡張カルマンフィルタ200、磁気・加速度センサデータ取得部202、角速度センサデータ取得部205、及び角速度データ誤差算出部208を備えている。また、拡張カルマンフィルタ200は、座標変換部201、観測更新演算部203、時間更新演算部207を備えている。また、図3は演算処理部104で行う処理をフローで示したものである。   FIG. 2 shows processing performed by the arithmetic processing unit 104 as functional blocks. The arithmetic processing unit 104 includes an extended Kalman filter 200, a magnetic / acceleration sensor data acquisition unit 202, an angular velocity sensor data acquisition unit 205, and an angular velocity data error calculation unit 208. The extended Kalman filter 200 includes a coordinate conversion unit 201, an observation update calculation unit 203, and a time update calculation unit 207. FIG. 3 is a flowchart showing processing performed by the arithmetic processing unit 104.

本発明の基本的な考えとして、動的加速度が姿勢センサシステムの観測ノイズに相当しているということに着目すると、ノイズを含む観測信号よりシステムの真の状態を推定することに適したカルマンフィルタを用いて姿勢推定アルゴリズムを構成することで誤差を低減化できると考えられる。   Focusing on the fact that the dynamic acceleration corresponds to the observation noise of the attitude sensor system as a basic idea of the present invention, a Kalman filter suitable for estimating the true state of the system from the observation signal including noise is provided. It is considered that the error can be reduced by using the attitude estimation algorithm.

[姿勢・方位推定アルゴリズム]
プロセスモデルの構築
カルマンフィルタを構成するためには、システムのプロセスモデルが必要となる。以下では、本発明で対象とするシステムの離散時間プロセスモデルを導出する。先ず状態方程式を求める。式(19)の右辺において、機体角速度ωbはジャイロセンサから得ることができるが、ジャイロセンサや加速度センサといった慣性センサには様々な要因による誤差があることが知られている。より精度の高い推定を行うためにはこのような誤差も推定することが望ましいため、ジャイロセンサ誤差を状態量として導入する。いま、ジャイロセンサより得られた角速度をωmeasure =[ωx ωy ωzT、ジャイロセンサの誤差をΔωb=[δωx δωy δωzTとすると、これらとωbとは以下のような関係となる。
ωmeasure =ωb+Δωb
・・・(22)
この式を式(19)に代入することで次式が得られる。
[Attitude and orientation estimation algorithm]
Construction of a process model To construct a Kalman filter, a process model of the system is required. In the following, a discrete time process model of the system targeted by the present invention is derived. First, the equation of state is obtained. On the right side of Equation (19), the airframe angular velocity ω b can be obtained from a gyro sensor, but it is known that inertial sensors such as a gyro sensor and an acceleration sensor have errors due to various factors. Since it is desirable to estimate such an error in order to perform more accurate estimation, a gyro sensor error is introduced as a state quantity. Assuming that the angular velocity obtained from the gyro sensor is ω measure = [ω x ω y ω z ] T and the error of the gyro sensor is Δω b = [δω x δω y δω z ] T , these and ω b are as follows: It becomes a relationship like this.
ω measure = ω b + Δω b
(22)
By substituting this equation into equation (19), the following equation is obtained.

一般的に、ジャイロセンサ誤差の要因としてはバイアス誤差、スケールファクタ誤差、ミスアラインメント誤差、ノイズ等が考えられ、これを式で表すと次式のようになる。
In general, as a cause of a gyro sensor error, a bias error, a scale factor error, a misalignment error, noise, and the like can be considered, and this is expressed by the following equation.

上式において、右辺第1項目がバイアス誤差、第2項目がスケールファクタ誤差、第3項目がミスアラインメント誤差、第4項目がノイズをそれぞれ表している。この中で、スケールファクタ誤差およびミスアラインメント誤差に関しては予め計測試験を行ったデータを基に補正することが可能であるため、リアルタイムに推定する必要があるのはバイアス誤差とノイズのみとなる。また、ノイズをリアルタイムに推定するのは不可能であるため、ここでは
として、バイアス誤差の推定のみを行う。
ジャイロセンサのバイアス誤差は経験的に次のようなダイナミクスを持つことが知られている。
In the above equation, the first item on the right side represents a bias error, the second item represents a scale factor error, the third item represents a misalignment error, and the fourth item represents noise. Among them, the scale factor error and the misalignment error can be corrected based on data that has been subjected to a measurement test in advance, so that only the bias error and noise need to be estimated in real time. Also, since it is impossible to estimate noise in real time, here
Only the bias error is estimated.
The bias error of the gyro sensor is empirically known to have the following dynamics.

ただしここでw = [wx wy wzTはホワイトノイズである。また、βx,βyzはパラメータであり、これらの値およびホワイトノイズwの分散は、静置試験のデータおよび式(25)の出力データのアラン分散[13]が一致するように調整し、決定した。各パラメータの値を示すと表3のようになる。
Here, w = [w x w y w z ] T is white noise. Also, β x , β y , β z are parameters, and these values and the variance of the white noise w are such that the static test data and the Allan variance [13] of the output data of Equation (25) match. Adjusted and decided. Table 3 shows the values of each parameter.

以上より、
を状態量としたときの状態方程式は次式となる。
ただし、ここで各行列は次式のようになっている。
From the above,
The state equation when is a state quantity is as follows.
However, here, each matrix is as follows.

式(26)を離散化すると式(30)となる。ただし、xtはtステップ目の状態量であり、Δtはサンプリングタイムである。ここではサンプリング周波数は50[Hz]であるため、Δt=0.02[s]となっている。
ここで、改めて
とおくと、離散時間の状態方程式は次式となる。
When equation (26) is discretized, equation (30) is obtained. However, x t is the state quantity at the t step, and Δt is the sampling time. Here, since the sampling frequency is 50 [Hz], Δt = 0.02 [s].
Here again
Then, the state equation in discrete time is as follows.

続いて、観測方程式を求める。機体に印加される動加速度を含む加速度誤差Δa=[δax δay δazTおよび地磁気誤差Δm=[δmx δmy δmzTを考えると、加速度センサ出力および磁気センサ出力は次式のように表すことができる。
ここで、
とすると、ytを出力とする観測方程式は次式となる。
ただし、ここでht(xt)は次式のようになっている。
Subsequently, the observation equation is obtained. Acceleration error Δa including dynamic acceleration applied to the body = [δa x δa y δa z ] Given the T and geomagnetic error Δm = [δm x δm y δm z] T, the acceleration sensor output and the magnetic sensor output is expressed by the following equation It can be expressed as
here,
When, observation equation to output a y t is expressed by the following equation.
Here, h t (x t ) is expressed by the following equation.

[拡張カルマンフィルタアルゴリズム]
先に示した式(31)、式(34)のような離散時間のプロセスモデルを得ることができた。しかし、両式はともに非線形方程式となっているため、線形カルマンフィルタアルゴリズムをそのまま適用することができない。このような場合、式(31)、式(34)を真値に対する推定誤差に関して線形化し、推定誤差を改めて状態量に選び、線形カルマンフィルタを適用して誤差推定カルマンフィルタを構成する方法と、非線形システムに線形カルマンフィルタを適用するための近似手法である拡張カルマンフィルタを構成するという2つの方法があるが、ここでは後者を選択した。以下では、拡張カルマンフィルタアルゴリズムについて説明を行う。いま、式(31)、式(34)のシステムに関して、tステップ目におけるxtの濾波推定値及び予備推定値を次のように表記したとき、
行列FtおよびHtを次式のように定義する。
[Extended Kalman filter algorithm]
It was possible to obtain a discrete-time process model such as Equation (31) and Equation (34) shown above. However, since both equations are nonlinear equations, the linear Kalman filter algorithm cannot be applied as it is. In such a case, a method of linearizing Equations (31) and (34) with respect to an estimation error with respect to a true value, selecting an estimation error as a state quantity, and applying a linear Kalman filter to configure an error estimation Kalman filter, and a nonlinear system There are two methods of constructing an extended Kalman filter which is an approximation method for applying a linear Kalman filter to the above, but the latter was selected here. Hereinafter, the extended Kalman filter algorithm will be described. Now, regarding the systems of Equation (31) and Equation (34), when the filtered estimated value and preliminary estimated value of x t at the t-th step are expressed as follows,
The matrices F t and H t are defined as follows:

以上のように定義された行列を用いた拡張カルマンフィルタアルゴリズムは以下に示す式で与えられる。
The extended Kalman filter algorithm using the matrix defined as above is given by the following equation.

ただしここで、式中のKtはカルマンゲイン、Pt/tおよびPt/t-1は推定誤差の共分散行列、Qtはシステムノイズの共分散行列、Rtは観測ノイズの共分散行列をそれぞれ表している。上記の式(37)〜(41)のアルゴリズムは、得られた観測値を用いて推定値を濾波する部分と、次ステップの推定値を予測する部分の2つから構成されており、それぞれ式(37)〜(39)が前者に当たり、式(40)〜(41)が後者に当たる。また、前者を観測更新、後者を時間更新と呼ぶこともある。以上のように与えられた式を順に計算することで、tステップ目におけるxtの最も確からしい濾波推定値
を得ることができる。
Where K t is Kalman gain, P t / t and P t / t-1 are covariance matrices of estimation errors, Q t is a covariance matrix of system noise, and R t is a covariance of observation noise. Each matrix is represented. The algorithms of the above equations (37) to (41) are composed of two parts, a part for filtering the estimated value using the obtained observation value and a part for predicting the estimated value of the next step. (37) to (39) correspond to the former, and expressions (40) to (41) correspond to the latter. In addition, the former may be called observation update and the latter may be called time update. By calculating the given formulas in order as described above, the most probable filtered estimated value of x t at the t-th step.
Can be obtained.

以下、図2及び図3を用いて、演算処理部104で行う処理について説明する。ステップ301において、座標変換部201は初期姿勢における磁気、加速度について、クォータニオンで表された式(11)を用いて上述のr-frameから b-frameへの座標変換を行う。クォータニオンの各要素は、Simple Rotationと呼ばれている姿勢表現方法を用いて式(6)で定義されている。式(6)は3つの虚数と1つの実数からなる4次のベクトル表現を用いた予め設定されるクォータニオンの初期値である式(1)及び式(10)を用いて表される。初期姿勢における磁気は、地磁気ベクトルの地上固定座標系(r-frame)成分であり、一方、初期姿勢における加速度は、重力加速度ベクトルの地上固定座標系(r-frame)成分である。   Hereinafter, processing performed by the arithmetic processing unit 104 will be described with reference to FIGS. 2 and 3. In step 301, the coordinate conversion unit 201 performs coordinate conversion from r-frame to b-frame described above using the equation (11) expressed in quaternions for the magnetism and acceleration in the initial posture. Each element of the quaternion is defined by Expression (6) using a posture expression method called “Simple Rotation”. Expression (6) is expressed using Expression (1) and Expression (10), which are preset quaternion initial values using a fourth-order vector expression composed of three imaginary numbers and one real number. The magnetism in the initial posture is the ground fixed coordinate system (r-frame) component of the geomagnetic vector, while the acceleration in the initial posture is the ground fixed coordinate system (r-frame) component of the gravitational acceleration vector.

次に、ステップ302では、予め設定されるクォータニオンの初期値を基に現ステップのセンサ座標系の磁気、加速度を推定することができ、これらを観測更新演算部203に入力する。   Next, in step 302, the magnetism and acceleration of the sensor coordinate system in the current step can be estimated based on the preset initial value of the quaternion, and these are input to the observation update calculation unit 203.

ステップ303では、磁気・加速度センサデータ取得部202は、磁気センサ103及び加速度センサ102から現ステップの磁気センサ出力(ベクトル)ameasure=[axyzT、加速度センサ出力(ベクトル)mmeasure=[mxyzTを取得する。取得する現ステップの磁気センサ出力、加速度センサ出力は、センサ座標系成分のデータであり、これらを観測更新演算部203に入力する。 In Step 303, the magnetic / acceleration sensor data acquisition unit 202 outputs the current sensor's magnetic sensor output (vector) a measure = [a x a y a z ] T , acceleration sensor output (vector) from the magnetic sensor 103 and the acceleration sensor 102. m measure = [m x m y m z ] T is acquired. The acquired magnetic sensor output and acceleration sensor output at the current step are sensor coordinate system component data, which are input to the observation update calculation unit 203.

ステップ304では、観測更新演算部203は、ステップ302で座標変換部201が算出した推定される現ステップの地磁気ベクトル推定値及び重力加速度ベクトル推定値と、ステップ303で磁気・加速度センサデータ取得部202から取得した現ステップの地磁気ベクトル及び重力加速度ベクトルとの誤差を求め、さらに式(34)を使って式(37)〜式(39)の演算を行って推定値を濾波する。   In step 304, the observation update calculation unit 203 calculates the estimated geomagnetic vector estimated value and gravitational acceleration vector estimated value of the current step calculated by the coordinate conversion unit 201 in step 302, and in step 303 the magnetic / acceleration sensor data acquisition unit 202. An error between the geomagnetic vector and the gravitational acceleration vector of the current step acquired from the above is obtained, and further, the estimated values are filtered by performing the calculations of Expressions (37) to (39) using Expression (34).

上述のステップ301〜ステップ304の処理とは別に、或いは並列的に、角速度の検出及びその誤差の算出が行われる。ステップ305において、角速度センサデータ取得部205は、角速度センサ101から現ステップの角速度を取得する。取得する現ステップの角速度は、センサ座標系成分のデータであり、これを角速度データ誤差算出部208に入力して角速度データの誤差を算出する。   The angular velocity is detected and its error is calculated separately or in parallel with the processing in steps 301 to 304 described above. In step 305, the angular velocity sensor data acquisition unit 205 acquires the angular velocity of the current step from the angular velocity sensor 101. The angular velocity of the current step to be acquired is sensor coordinate system component data, which is input to the angular velocity data error calculation unit 208 to calculate the angular velocity data error.

ステップ306では、角速度データ誤差算出部208は角速度センサ101の誤差であるバイアス誤差、スケールファクタ誤差、ミスアライメント誤差を考慮して角速度データの誤差求める。この角速度データの誤差算出に当たっては、上述の式(22)〜式(31)の演算を行い、式(26)の状態方程式、およびこれをさらに変形させた式(31)の離散時間の状態方程式を求める。   In step 306, the angular velocity data error calculation unit 208 calculates an angular velocity data error in consideration of a bias error, a scale factor error, and a misalignment error that are errors of the angular velocity sensor 101. In calculating the error of the angular velocity data, the above-described equations (22) to (31) are calculated, and the state equation of equation (26) and the state equation of discrete time of equation (31) obtained by further modifying the equation (26). Ask for.

先のステップ304における演算結果とステップ306における角速度データの誤差算出結果とは加算されて時間更新演算部207へ入力される。時間更新演算部207では、式(26)を使って式(40)及び式(41)の演算を行い推定値を濾波する。以上の演算処理の結果クォータニオン補正もなされ、次のステップにおけるxtの最も確からしい濾波推定値を得ることができる。 The calculation result in the previous step 304 and the error calculation result of the angular velocity data in step 306 are added and input to the time update calculation unit 207. The time update calculation unit 207 performs calculations of Expression (40) and Expression (41) using Expression (26), and filters the estimated value. As a result of the above arithmetic processing, quaternion correction is also performed, and the most probable filtered estimated value of x t in the next step can be obtained.

最後に、ステップ308において、上記のようにして得られた算出された次ステップの濾波推定値は、座標変換部201に入力される。すなわち、座標変換部201では、算出された次ステップの濾波推定値を用いて次ステップの座標変換行列を算出し、その結果は観測更新演算部203に入力されてクォータニオン推定値として小型姿勢センサの最終出力となる一方、さらに次のステップにおける観測更新演算処理、さらには時間更新演算処理に用いる。上記ステップ301乃至ステップ308の処理が、小型姿勢センサの演算処理部104で行う一ステップの処理であり、これを繰り返し行う。   Finally, in step 308, the calculated filtered estimated value of the next step obtained as described above is input to the coordinate conversion unit 201. That is, the coordinate conversion unit 201 calculates the next step coordinate conversion matrix using the calculated next-step filtered estimated value, and the result is input to the observation update calculation unit 203 as a quaternion estimated value. While it is the final output, it is used for observation update calculation processing in the next step, and further for time update calculation processing. The processing from step 301 to step 308 is one-step processing performed by the arithmetic processing unit 104 of the small posture sensor, which is repeated.

[シミュレーション1]
本発明の主目的である動的加速度環境下における姿勢誤差の低減化を実現するためには、観測ノイズの共分散行列Rtの決め方が重要になる。ここでは、Rtの値をシミュレーションにより決定する。先ず動的加速度による誤差が少なくなるようにRtの値を決定した。このときの値をR1とする。図4にR1によるシミュレーション結果を示す。図4(a)が姿勢を変動させたとき、図4(b)がセンサを水平に保った状態で、図19のb-frame上のXb軸方向、すなわちピッチ角度に大きな誤差が生じやすい方向に正弦波状の動的加速度を印加したときのシミュレーション結果である。また、それぞれ破線が従来のアルゴリズムによる結果、実線が拡張カルマンフィルタによる結果、点線がその際の真の姿勢値である。図4より、R1を用いると動的加速度による誤差は低減化するが、姿勢の追従が極端に遅くなってしまうことが分かる。一方、姿勢の追従を充分に速くなるようなRtの値をR2として、R2による結果を図5に示す。図5よりR2を用いると姿勢の追従はよくなるが、動的加速度誤差が低減化できないことが分かる。そこで、上述のR1,R2を用いてRtを次式のような関数とした。ただし、KRは適当なゲインである。また、Tに関してはシミュレーションを繰り返すことによって適当な値を決定した。この事例では、Tの値は11[s]となっている。
[Simulation 1]
In order to realize the reduction of the attitude error under the dynamic acceleration environment, which is the main object of the present invention, it is important to determine the covariance matrix R t of the observation noise. Here, determined by simulation the value of R t. First, the value of R t was determined so that the error due to dynamic acceleration was reduced. The value at this time is R 1 . FIG. 4 shows a simulation result by R 1 . When FIG. 4A changes the posture, FIG. 4B keeps the sensor horizontal, and the Xb axis direction on the b-frame in FIG. 19, that is, the direction in which a large error is likely to occur in the pitch angle. It is a simulation result when sinusoidal dynamic acceleration is applied to. Also, the broken lines are the results of the conventional algorithm, the solid lines are the results of the extended Kalman filter, and the dotted lines are the true posture values at that time. From FIG. 4, it can be seen that when R 1 is used, errors due to dynamic acceleration are reduced, but posture tracking becomes extremely slow. On the other hand, the value of such R t is sufficiently fast tracking attitude as R 2, shows the results of R 2 in FIG. From FIG. 5, it can be seen that using R 2 improves the tracking of the posture, but the dynamic acceleration error cannot be reduced. So, the R t and function as the following equation using the R 1, R 2 described above. However, K R is an appropriate gain. For T, an appropriate value was determined by repeating the simulation. In this example, the value of T is 11 [s].

このときのシミュレーション結果を図6に示す。以上のようなRtを用いることで、姿勢の追従が速く、動的加速度による誤差は低減化できるような結果を得ることができた。続いて図7〜図9に線形誤差推定カルマンフィルタ(KF)と本発明において用いた拡張カルマンフィルタ(EKF)の比較を示す。これらのうち、図7は、移動物体に搭載した小型姿勢センサに、線形誤差推定カルマンフィルタ(KF)を用いた場合と、本発明の拡張カルマンフィルタ(EKF)用いた場合におけるロール角についてのシミュレーション結果を比較して示す図である。図7中において、点線で示された線はKFによるシミュレーション結果であり、実線で示された線はEKFによるシミュレーション結果である。両者を比較すると、KFの結果はEKFと比べて振動的になってしまっていることが分かる。これは、式(34)の観測方程式に含まれる座標変換の非線形性をKFでは充分に考慮できていないためであると考えられる。 The simulation result at this time is shown in FIG. By using the Rt as described above, it was possible to obtain a result that the posture can be followed quickly and errors due to dynamic acceleration can be reduced. 7 to 9 show a comparison between the linear error estimation Kalman filter (KF) and the extended Kalman filter (EKF) used in the present invention. Among these, FIG. 7 shows the simulation results of the roll angle when the linear error estimation Kalman filter (KF) is used for the small posture sensor mounted on the moving object and when the extended Kalman filter (EKF) of the present invention is used. It is a figure shown in comparison. In FIG. 7, a line indicated by a dotted line is a simulation result by KF, and a line indicated by a solid line is a simulation result by EKF. Comparing both, it can be seen that the result of KF has become more vibrant than that of EKF. This is presumably because the nonlinearity of the coordinate transformation included in the observation equation of Equation (34) cannot be sufficiently considered by KF.

図8は、移動物体に搭載した小型姿勢センサに、線形誤差推定カルマンフィルタ(KF)を用いた場合と、本発明の拡張カルマンフィルタ(EKF)用いた場合におけるピッチ角についてのシミュレーション結果を比較して示す図である。図8中において、点線で示された線はKFによるシミュレーション結果であり、実線で示された線はEKFによるシミュレーション結果である。両者を比較すると、図7の場合と同様に、KFの結果はEKFと比べて振動的になってしまっていることが分かる。   FIG. 8 shows a comparison of simulation results for pitch angles when a linear error estimation Kalman filter (KF) is used for a small attitude sensor mounted on a moving object and when an extended Kalman filter (EKF) of the present invention is used. FIG. In FIG. 8, the line indicated by the dotted line is the simulation result by KF, and the line indicated by the solid line is the simulation result by EKF. Comparing the two, it can be seen that the result of KF is more oscillating than that of EKF, as in FIG.

図9は、移動物体に搭載した小型姿勢センサに、線形誤差推定カルマンフィルタ(KF)を用いた場合と、本発明の拡張カルマンフィルタ(EKF)用いた場合におけるヨー角についてのシミュレーション結果を比較して示す図である。図9中において、点線で示された線はKFによるシミュレーション結果であり、実線で示された線はEKFによるシミュレーション結果である。両者を比較すると、図7及び図8の場合と同様に、KFの結果はEKFと比べて振動的になってしまっていることが分かる。   FIG. 9 shows a comparison of simulation results for yaw angles when a linear error estimation Kalman filter (KF) is used for a small posture sensor mounted on a moving object and when an extended Kalman filter (EKF) of the present invention is used. FIG. In FIG. 9, the line indicated by a dotted line is a simulation result by KF, and the line indicated by a solid line is a simulation result by EKF. When both are compared, it can be seen that the result of KF is more oscillating than that of EKF, as in FIGS.

[実験1]
従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、リアルタイムで姿勢・方位推定を行った結果を図10〜図13に示す。比較対象として、Crossbow社製の高精度姿勢センサAHRS400を選択した。このAHRS400は移動体に搭載することを目的とした姿勢センサであり、動的加速度を印加しても姿勢誤差が生じにくいという特徴を持っている。しかし、その重量は700(g)であり、MARGセンサの10倍以上となっている。本来であればAHRS400とMARGセンサをともに小型無人ヘリコプタ上のマウントに搭載し、飛行しつつ姿勢データの比較を行うことが望ましいが、AHRS400の重量の関係上それは難しいため、この実験では両方のセンサを同一マウント上に搭載した状態で、そのマウントを地上で動かすという試験方法をとった。ちなみに、比較を行いやすくするために、得られたクォータニオン推定値をオイラー角に変換した後にプロットしている。図10は、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るような角速度環境下で、ロール角についてリアルタイムで姿勢・方位推定を行った結果を示す図である。この図において、実線は拡張カルマンフィルタアルゴリズムによる方位推定結果を示し、細かい点線は従来のアルゴリズムによる方位推定結果を示す。また、粗い目の鎖線は上記高精度姿勢センサAHRS400による方位推定結果を示す。これらのグラフ図から、本発明の拡張カルマンフィルタアルゴリズムは高精度の姿勢・方位推定結果が得られることが分かる。
[Experiment 1]
The conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer and the results of posture / orientation estimation in real time are shown in FIGS. As a comparison target, a high-precision posture sensor AHRS400 manufactured by Crossbow was selected. The AHRS 400 is a posture sensor intended to be mounted on a moving body, and has a feature that a posture error hardly occurs even when dynamic acceleration is applied. However, the weight is 700 (g), which is more than 10 times that of the MARG sensor. Originally, it would be desirable to have both AHRS400 and MARG sensors mounted on a small unmanned helicopter mount and compare attitude data while flying, but because of the weight of AHRS400, it is difficult to do so. The test method of moving the mount on the ground while mounted on the same mount was used. Incidentally, in order to facilitate comparison, the obtained quaternion estimated value is plotted after being converted to Euler angle. FIG. 10 shows a conventional algorithm and the extended Kalman filter algorithm of the present invention mounted on a microcomputer, and in real-time posture and roll angle in an angular velocity environment that greatly exceeds the angular velocity generated when a small unmanned helicopter moves. It is a figure which shows the result of having performed azimuth | direction estimation. In this figure, the solid line shows the direction estimation result by the extended Kalman filter algorithm, and the fine dotted line shows the direction estimation result by the conventional algorithm. A coarse chain line indicates a direction estimation result obtained by the high-accuracy posture sensor AHRS400. From these graphs, it can be seen that the extended Kalman filter algorithm of the present invention can obtain a highly accurate posture / orientation estimation result.

図11は、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るような角速度環境下で、ピッチ角についてリアルタイムで姿勢・方位推定を行った結果を示す図である。グラフ図の線の内訳については図10と同じである。また図11からも、本発明の拡張カルマンフィルタアルゴリズムは高精度の姿勢・方位推定結果が得られることが分かる。図12は、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るような角速度環境下で、ヨー角についてリアルタイムで姿勢・方位推定を行った結果を示す図である。グラフ図の線の内訳については図10、図11と同じである。また図12からも、本発明の拡張カルマンフィルタアルゴリズムは高精度の姿勢・方位推定結果が得られることが分かる。図13は、従来のアルゴリズムおよび本発明の拡張カルマンフィルタアルゴリズムをマイクロコンピュータ上に実装し、センサを水平に保った状態で、シミュレーションと同様に、小型無人ヘリコプタの飛行中に加わる加速度と比較して充分大きく、図19のb-frame上Xb軸方向に対して正弦波状の動的加速度を印加し、ピッチ角についてリアルタイムで姿勢・方位推定を行った結果を示す図である。グラフ図の線の内訳については図10−図12と同じである。また図13からも、本発明の拡張カルマンフィルタアルゴリズムは高精度の姿勢・方位推定結果が得られることが分かる。   FIG. 11 shows the pitch angle in real time in an angular velocity environment in which the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer, and greatly exceed the angular velocity generated when a small unmanned helicopter moves. It is a figure which shows the result of having performed azimuth | direction estimation. The breakdown of the lines in the graph is the same as in FIG. FIG. 11 also shows that the extended Kalman filter algorithm of the present invention can obtain a highly accurate posture / orientation estimation result. FIG. 12 shows the yaw angle in real time in an angular velocity environment in which the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer and greatly exceed the angular velocity generated when a small unmanned helicopter moves. It is a figure which shows the result of having performed azimuth | direction estimation. The breakdown of the lines in the graph is the same as in FIGS. Also from FIG. 12, it can be seen that the extended Kalman filter algorithm of the present invention can obtain a highly accurate posture / orientation estimation result. FIG. 13 shows that the conventional algorithm and the extended Kalman filter algorithm of the present invention are mounted on a microcomputer and the sensor is kept horizontal, as in the simulation, compared with the acceleration applied during the flight of a small unmanned helicopter. FIG. 20 is a diagram showing a result of performing posture / orientation estimation in real time for a pitch angle by applying a sinusoidal dynamic acceleration in the Xb axis direction on the b-frame in FIG. The breakdown of the lines in the graph is the same as in FIGS. Also from FIG. 13, it can be seen that the extended Kalman filter algorithm of the present invention can obtain a highly accurate posture / orientation estimation result.

このときの実験では、図10〜図12に関しては小型無人ヘリコプタの運動の際に発生する角速度を大きく上回るように最大250[deg/s]の角速度環境下において試験を行っており、図13に関しては、小型無人ヘリコプタの飛行中に加わる加速度と比較して充分大きな加速度を印加しているが、これは、小型姿勢センサを小型無人ヘリコプタに搭載することを想定したからである。表4に小型無人ヘリコプタの各飛行状態およびこの試験における加速度のRMS値を示している。
In the experiment at this time, with respect to FIGS. 10 to 12, tests were performed in an angular velocity environment of a maximum of 250 [deg / s] so as to greatly exceed the angular velocity generated during the movement of the small unmanned helicopter. Applies a sufficiently large acceleration compared to the acceleration applied during the flight of the small unmanned helicopter, because it is assumed that the small attitude sensor is mounted on the small unmanned helicopter. Table 4 shows each flight state of the small unmanned helicopter and the acceleration RMS value in this test.

それぞれ、「Hover」がホバリング飛行中、「Circular flight」が比較的大きな加速度が印加される旋回飛行中、「Experiment」が図13に示す試験中の値である。図10〜図12より、概ねAHRS400とMARGセンサのデータは一致しており、正確な推定が行われているということが分かる。しかし、図12の6〜16秒付近においてAHRS400と比較して提案手法および従来手法のデータが誤差を持っていることが分かる。これは磁気外乱の影響であると考えられ、小型無人ヘリコプタの飛行する屋外では特に大きな問題ではないが、今後、このような磁気外乱の影響も低減化するように改良する必要があると考えられる。また、AHRS400のデータが一定値となっている箇所がいくつかあるが、これは、AHRS400は角速度が200[deg/s]以下でしか姿勢の推定を行えず、200[deg/s]以上の角速度を与えると姿勢の推定が停止してしまうためである。これに対して、MARGセンサでは200[deg/s]以上の角速度を与えた際にも姿勢の推定を行うことができている。一方、動的加速度を印加しても姿勢誤差を充分に低減化できており、仕様を充分に満たしていることが分かる。   “Experiment” is a value during the test shown in FIG. 13, while “Hover” is during hovering flight, “Circular flight” is during a turning flight to which a relatively large acceleration is applied. 10 to 12, it can be seen that the data of the AHRS 400 and the MARG sensor are almost the same, and accurate estimation is performed. However, it can be seen that the data of the proposed method and the conventional method have an error in comparison with the AHRS 400 in the vicinity of 6 to 16 seconds in FIG. This is considered to be the effect of magnetic disturbance, and is not a major problem outdoors where small unmanned helicopters fly, but it is thought that it will be necessary to improve in the future to reduce the influence of such magnetic disturbance. . In addition, there are some places where the data of AHRS 400 has a constant value. This is because AHRS 400 can estimate the posture only when the angular velocity is 200 [deg / s] or less, and is 200 [deg / s] or more. This is because posture estimation stops when angular velocity is given. On the other hand, the MARG sensor can estimate the posture even when an angular velocity of 200 [deg / s] or more is given. On the other hand, it can be seen that the posture error can be sufficiently reduced even when dynamic acceleration is applied, and the specification is sufficiently satisfied.

以上説明したように、本発明によれば、先行技術において開発されたMARGセンサの問題点である動的加速度環境下での姿勢誤差を低減化するために、慣性センサのバイアス推定を含む拡張カルマンフィルタを姿勢推定アルゴリズムとして適用し、動的加速度環境下での精度の高い姿勢・方位推定を実現した。   As described above, according to the present invention, the extended Kalman filter including the bias estimation of the inertial sensor is used to reduce the posture error under the dynamic acceleration environment, which is a problem of the MARG sensor developed in the prior art. Was applied as a posture estimation algorithm, and the posture and orientation estimation with high accuracy under the dynamic acceleration environment was realized.

(実施の形態2)
本発明の第2の実施の形態は、加速度環境下における高精度の姿勢・方位推定を実現するに当って、低周波数の加速度環境下での姿勢推定誤差の低減化を図るものである。
(Embodiment 2)
The second embodiment of the present invention is intended to reduce posture estimation errors under a low-frequency acceleration environment in order to realize highly accurate posture / orientation estimation under an acceleration environment.

[加速度信号のモデル化と拡張カルマンフィルタへの組込]
第1の実施の形態の方法においては、加速度環境下における高精度の姿勢・方位推定を或る程度は実現している。しかしながら、一定値に近い低周波数の加速度外乱が印加されると姿勢推定値がドリフトしてしまう。そこで、低周波数の加速度外乱をシステムの状態の一部として推定することで姿勢誤差を低減化する方法を以下に提案する。加速度外乱を推定するためのもう一つの手法として、カルマンフィルタの状態量に加速度外乱を含めることで、姿勢と加速度外乱を同時に推定する手法を提案する。移動や回転運動に伴って発生する加速度外乱は加速度センサのサンプリング周期に比べて十分に低周波であることから、サンプリング間における加速度外乱のダイナミクスを下記の式で示す。
[Modeling of acceleration signal and incorporation into extended Kalman filter]
In the method of the first embodiment, high-precision posture / orientation estimation under an acceleration environment is realized to some extent. However, when a low-frequency acceleration disturbance close to a certain value is applied, the estimated posture value drifts. Therefore, a method for reducing the attitude error by estimating low-frequency acceleration disturbance as a part of the system state is proposed below. As another method for estimating acceleration disturbance, we propose a method to estimate posture and acceleration disturbance at the same time by including acceleration disturbance in the state quantity of Kalman filter. Since the acceleration disturbance generated due to the movement and rotation motion is sufficiently low in frequency compared to the sampling period of the acceleration sensor, the dynamics of the acceleration disturbance during sampling is expressed by the following equation.

まず低周波数の加速度外乱を
dr=[abxbybzT
とし、外乱のダイナミクスを次式のように表す。
式(44)より、
を状態量としたときの状態方程式は次の式(45)のように置き換えられる。
First, the low-frequency acceleration disturbance is expressed as a dr = [a bx a by a bz ] T
And the disturbance dynamics is expressed as
From equation (44)
The state equation where is a state quantity is replaced as the following equation (45).

ただし、ここで各行列は次式のようになっている。
However, here, each matrix is as follows.

式(45)を離散化すると下記の式(49)のようになる。ただし、xt はt ステップ目の状態量であり、Δt はサンプリングタイムである。
ここで、改めて
とおくと、離散時間の状態方程式は次式となる。
When the expression (45) is discretized, the following expression (49) is obtained. Here, x t is the state quantity at the t step, and Δt is the sampling time.
Here again
Then, the state equation in discrete time is as follows.

続いて、観測方程式を求める。式(45)で推定された低周波数の加速度外乱を観測方程式に含めると次式のように書くことができる。
ここで、
とすると、yt を出力とする観測方程式は次式となる。
ただし、ここでht(xt)は次式のようになっている。
Subsequently, the observation equation is obtained. When the low-frequency acceleration disturbance estimated by the equation (45) is included in the observation equation, it can be written as the following equation.
here,
When, observation equation to output a y t is expressed by the following equation.
Here, h t (x t ) is expressed by the following equation.

第1の実施の形態の方法と同様に行列Ft及びHtを次式のように定義する。
以上のように定義された行列と式(50)、式(53)を用いた拡張カルマンフィルタアルゴリズムは以下に示す式で与えられる。
Similar to the method of the first embodiment, the matrices F t and H t are defined as follows:
The extended Kalman filter algorithm using the matrix defined above and the equations (50) and (53) is given by the following equation.

[シミュレーション2]
以下においては、本発明の拡張フィルタアルゴリズムによって正常に低周波数の加速度の外乱が推定されていることを確認するために移動体に一定値の低周波数の加速度外乱を加え、それが正しく推定されていることをシミュレーションにより確認する。本シミュレーションでは、機体が停止状態から最大加速度で加速し、一定時間の移動の後、最大加速度で原則すると仮定し、印加される外乱の大きさは0.2Gとした。また、各アルゴリズムの性能を明確にするため、加減速による加速度外乱の印加時間は実フライトよりも十分に長い20秒間に設定した。図14は本シミュレーションに用いた各軸の加速度データを示す。図14中において、実線はX軸方向の加速度データを示し、一点鎖線はY軸方向の加速度データを示し、点線はZ軸方向の加速度データを示す。図14により、開始から約15秒後にx軸加速度に一定値の外乱が負方向及び正方向に1回ずつ印加されていることが分かる。y軸には加速度がなく、z軸には重力加速度が印加されている。図15は各アルゴリズムの姿勢推定結果を示し、図15(a)はロール方向における姿勢推定結果を示し、図15(b)はピッチ方向における姿勢推定結果を示し、図15(c)はヨー方向における姿勢推定結果を示す。また、図15において、Trueの表示は姿勢推定結果の真値を実線で示し、Prevの表示は従来のアルゴリズムによる姿勢推定結果を点線(細線)示し、Newの表示は本発明において状態量に加速度外乱を含めた拡張カルマンフィルタのアルゴリズムによる姿勢推定結果を点線(太線)示す。図16は本発明のアルゴリズムにおける加速度外乱推定値と加速度外乱の真値との比較を示す図であり、図16(a)はX方向における加速度外乱の比較結果を示し、図16(b)はY方向における加速度外乱の比較結果を示し、図16(c)はZ方向における加速度外乱の比較結果を示す。図16において、Trueの表示は加速度外乱の真値を実線で示し、Newの表示は本発明のアルゴリズムによる加速度外乱推定結果を点線(太線)で示す。
[Simulation 2]
In the following, in order to confirm that the low-frequency acceleration disturbance is normally estimated by the extended filter algorithm of the present invention, a constant low-frequency acceleration disturbance is applied to the moving object, and the low-frequency acceleration disturbance is correctly estimated. Confirm by simulation. In this simulation, it is assumed that the aircraft accelerates at a maximum acceleration from a stationary state, moves for a certain period of time, and then, as a rule, uses the maximum acceleration, and the magnitude of the applied disturbance is 0.2G. Further, in order to clarify the performance of each algorithm, the application time of acceleration disturbance due to acceleration / deceleration was set to 20 seconds, which is sufficiently longer than the actual flight. FIG. 14 shows acceleration data of each axis used in this simulation. In FIG. 14, the solid line indicates acceleration data in the X-axis direction, the alternate long and short dash line indicates acceleration data in the Y-axis direction, and the dotted line indicates acceleration data in the Z-axis direction. FIG. 14 shows that a constant disturbance is applied to the x-axis acceleration once in the negative direction and the positive direction about 15 seconds after the start. There is no acceleration on the y-axis, and gravitational acceleration is applied on the z-axis. FIG. 15 shows the posture estimation result of each algorithm, FIG. 15A shows the posture estimation result in the roll direction, FIG. 15B shows the posture estimation result in the pitch direction, and FIG. 15C shows the yaw direction. The posture estimation result in is shown. In FIG. 15, True indicates the true value of the posture estimation result with a solid line, Prev indicates the posture estimation result according to the conventional algorithm with a dotted line (thin line), and New displays the state quantity in the present invention. The posture estimation result by the extended Kalman filter algorithm including disturbance is shown by a dotted line (thick line). FIG. 16 is a diagram showing a comparison between the acceleration disturbance estimated value and the true value of the acceleration disturbance in the algorithm of the present invention. FIG. 16A shows a comparison result of the acceleration disturbance in the X direction, and FIG. The comparison result of the acceleration disturbance in the Y direction is shown, and FIG. 16C shows the comparison result of the acceleration disturbance in the Z direction. In FIG. 16, True displays the true value of the acceleration disturbance with a solid line, and New displays the acceleration disturbance estimation result according to the algorithm of the present invention with a dotted line (thick line).

上述のシミュレーションにおいて、計測される加速度は重力と移動による加速度外乱であり、図15の真値から姿勢には変動がないことが分かる。図14の加速度と、重力の各軸成分が、(gx gy gz)=[0,0,1][G]であることから、時刻15[sec]から35[sec]にかけてX軸には加速による−0.2[G]の外乱が印加され、また、時刻50[sec]から70[sec]にかけてX軸には減速による0.2[G]の外乱が印加されているが、図16の本発明のアルゴリズムによる加速度外乱推定値は真値によく一致している。図15の各アルゴリズムによる姿勢推定値を比較すると、従来アルゴリズムによる姿勢推定ではピッチ角に最大12[deg]の誤差を生じているのに対して、提案した2種類のアルゴリズムでは共に誤差が0.1[deg]以内に収まっており、その有効性が確認できる。従来アルゴリズムにおける姿勢推定誤差は、印加された加速度外乱をピッチ角が傾いたときにX軸へ印加される重力成分として扱ったために生じていると考えられる。 In the above simulation, the measured acceleration is an acceleration disturbance due to gravity and movement, and it can be seen from the true values in FIG. Since each axis component of acceleration and gravity in FIG. 14 is (g x g y g z ) = [0, 0, 1] [G], the X axis extends from time 15 [sec] to 35 [sec]. A disturbance of -0.2 [G] due to acceleration is applied to the X axis, and a disturbance of 0.2 [G] due to deceleration is applied to the X axis from time 50 [sec] to 70 [sec]. The acceleration disturbance estimated value according to the algorithm of the present invention in FIG. 16 agrees well with the true value. Comparing the estimated pose values of the algorithms shown in FIG. 15, the pose estimation using the conventional algorithm has a maximum error of 12 [deg] in the pitch angle, whereas the two proposed algorithms have an error of 0. It is within 1 [deg], and its effectiveness can be confirmed. The posture estimation error in the conventional algorithm is considered to be caused by treating the applied acceleration disturbance as a gravity component applied to the X axis when the pitch angle is tilted.

[実験2]
実際の加速度外乱環境下における測定結果を基に、従来の姿勢推定アルゴリズムと本発明の姿勢推定アルゴリズムの性能比較を行う。アルゴリズムを実装するハードウェアとして本発明の小型姿勢センサを用いるが、この小型姿勢センサに搭載されているマイクロプロセッサ(MPU)は小型のため、複数種類のアルゴリズムを同時にリアルタイムで演算することはできない。そこで、カルマンフィルタを用いた加速度外乱推定を行うアルゴリズムを実装し、測定によって得られた加速度、角速度、磁気データを基にオフラインで推定演算を行った。演算に用いる行列の次数からカルマンフィルタを用いた加速度外乱推定アルゴリズムの演算量が多いのは明らかであるため、本手法がMPUに実装可能であれば、他の手法においても十分に実現可能である。姿勢推定アルゴリズムの比較検証を行うための真値は、Crossbow社製の高精度姿勢センサAHRS440から取得する。
[Experiment 2]
Based on the measurement results in an actual acceleration disturbance environment, the performance comparison between the conventional posture estimation algorithm and the posture estimation algorithm of the present invention is performed. The small attitude sensor of the present invention is used as hardware for implementing the algorithm. However, since the microprocessor (MPU) mounted on the small attitude sensor is small, a plurality of types of algorithms cannot be calculated simultaneously in real time. Therefore, we implemented an algorithm for estimating the acceleration disturbance using the Kalman filter, and performed an off-line estimation calculation based on the acceleration, angular velocity, and magnetic data obtained by the measurement. Since it is clear that the calculation amount of the acceleration disturbance estimation algorithm using the Kalman filter is large from the order of the matrix used for the calculation, if this method can be implemented in the MPU, it can be sufficiently realized by other methods. A true value for performing comparison verification of the posture estimation algorithm is acquired from a high-precision posture sensor AHRS440 manufactured by Crossbow.

この実験では、車載による姿勢計測実験を行うこととし、一般車両に各センサ及びデータ収集用ノートパソコン、電源を搭載し、走行実験を行った結果を図17及び図18に示す。本来ならば、上記センサ等は小型無人ヘリコプターに搭載して測定することが望ましいが、姿勢センサAHRS440は小型無人ヘリコプターへの搭載に適切な重量ではない。そこで、移動体へ適用した際の有効性を確かめる手段として、車両への搭載を選択した。試験は平坦な道路上を直線的に走行し、40km/hへ到達した段階で、急減速により停止する。このとき、センサの姿勢変動はほとんどなく、前後方向に大きな加速度外乱が印加される。   In this experiment, an on-vehicle attitude measurement experiment is performed, and the results of a running experiment in which each sensor, a data collection notebook personal computer, and a power source are mounted on a general vehicle are shown in FIGS. 17 and 18. Originally, it is desirable that the above-described sensor or the like is mounted on a small unmanned helicopter for measurement, but the attitude sensor AHRS 440 is not of an appropriate weight for mounting on a small unmanned helicopter. Therefore, as a means of confirming the effectiveness when applied to a moving body, it was selected to be mounted on a vehicle. The test runs linearly on a flat road and stops due to sudden deceleration when it reaches 40 km / h. At this time, there is almost no change in the attitude of the sensor, and a large acceleration disturbance is applied in the front-rear direction.

図17及び図18は上記実験におけるセンサによる姿勢計測結果を表す図である。図17は各アルゴリズムの姿勢計測結果を示し、図17(a)はロール方向における姿勢計測結果を示し、図17(b)はピッチ方向における姿勢計測結果を示し、図17(c)はヨー方向における姿勢計測結果を示す。また、図17において、Refの表示はCrossbow社製の高精度姿勢センサAHRS440による姿勢計測結果(真値)を実線で示し、Prevの表示は従来のアルゴリズムによる姿勢計測結果を点線(細線)で示し、Newの表示は本発明の拡張カルマンフィルタのアルゴリズムによる姿勢計測結果を点線(太線)で示す。図18は本発明のアルゴリズムにおける加速度外乱推定値と加速度外乱の真値との比較を示す図であり、図18(a)はX方向における加速度外乱の比較結果を示し、図18(b)はY方向における加速度外乱の比較結果を示し、図18(c)はZ方向における加速度外乱の比較結果を示す。図18において、Trueの表示は加速度外乱の真値を実線で示し、Newの表示は本発明のアルゴリズムによる加速度外乱推定結果を点線(太線)で示す。   FIGS. 17 and 18 are diagrams showing the result of posture measurement by the sensor in the above experiment. FIG. 17 shows the attitude measurement results of each algorithm, FIG. 17A shows the attitude measurement results in the roll direction, FIG. 17B shows the attitude measurement results in the pitch direction, and FIG. 17C shows the yaw direction. The attitude measurement result in is shown. In FIG. 17, the Ref display indicates the posture measurement result (true value) by the high accuracy posture sensor AHRS440 manufactured by Crossbow, and the Prev display indicates the posture measurement result by the conventional algorithm by a dotted line (thin line). , New indicates the posture measurement result by the extended Kalman filter algorithm of the present invention with a dotted line (thick line). FIG. 18 is a diagram showing a comparison between the acceleration disturbance estimated value and the true value of the acceleration disturbance in the algorithm of the present invention. FIG. 18A shows a comparison result of the acceleration disturbance in the X direction, and FIG. The comparison result of the acceleration disturbance in the Y direction is shown, and FIG. 18C shows the comparison result of the acceleration disturbance in the Z direction. In FIG. 18, True indicates the true value of the acceleration disturbance with a solid line, and New indicates the acceleration disturbance estimation result according to the algorithm of the present invention with a dotted line (thick line).

図17から、従来アルゴリズムでは加速度外乱の印加に応じてピッチ角に最大18[deg]の誤差を生じており、さらに8秒付近からヨー角の揺動も見られる。他方、本発明の拡張カルマンフィルタのアルゴリズムによる加速度外乱推定では図18から明らかなように、Crossbow社製の高精度姿勢センサAHRS440と同等の外乱抑制効果が得られており、その結果図17においても、ロール角(図17(a))、ピッチ角(図17(b))、ヨー角(図17(c))のいずれにおいても誤差が0.1[deg]以内に収まっている。これにより、本発明の拡張カルマンフィルタのアルゴリズムによれば、正しく姿勢推定が行われていることが分かり、その有効性が確認できる。   From FIG. 17, in the conventional algorithm, an error of a maximum of 18 [deg] is generated in the pitch angle in accordance with the application of the acceleration disturbance, and the yaw angle fluctuates from around 8 seconds. On the other hand, in the acceleration disturbance estimation by the extended Kalman filter algorithm of the present invention, as is apparent from FIG. 18, the disturbance suppression effect equivalent to the high accuracy attitude sensor AHRS440 manufactured by Crossbow is obtained. In any of the roll angle (FIG. 17A), pitch angle (FIG. 17B), and yaw angle (FIG. 17C), the error is within 0.1 [deg]. Thus, according to the extended Kalman filter algorithm of the present invention, it can be seen that the posture estimation is correctly performed, and the effectiveness thereof can be confirmed.

Claims (4)

移動物体の独立な3軸上の角速度を計測する角速度センサと、
移動物体の独立な3軸上の加速度を計測する加速度センサと、
移動物体の独立な3軸上の地磁気を計測する磁気センサと、
前記角速度センサ、前記加速度センサ、及び前記磁気センサの計測値を基にクォータニオンの推定値を算出する演算処理部とを備え、
前記演算処理部は、予め推定されたクォータニオン推定値を用いて座標変換行列を算出し、当該座標変換行列を用いて地上固定座標系成分の地磁気ベクトル及び重力加速度ベクトルを変換し、現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルを推定する座標変換部と、
前記磁気センサ及び前記加速度センサから現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルを取得する磁気・加速度センサデータ取得部と、
前記座標変換部で推定した現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルと、前記磁気・加速度センサデータ取得部で取得した現ステップのセンサ座標系成分の地磁気ベクトル及び重力加速度ベクトルとの誤差を算出するベクトル誤差算出部と、
前記角速度センサから現ステップのセンサ座標系成分の角速度ベクトルを取得する角速度センサデータ取得部と、
前記取得した角速度データについて、バイアス誤差、スケールファクタ誤差、ミスアラインメント誤差、ノイズの要素を含む角速度センサ誤差を演算により求める角速度データ誤差算出部と、
前ステップのクォータニオンと現ステップのクォータニオンの時間変化を求めるクォータニオン算出部と
前記座標変換部、前記磁気・加速度センサデータ取得部、前記ベクトル誤差算出部、前記角速度データ誤差算出部、及び前記クォータニオン算出部を基に現在の姿勢を表すクォータニオンおよび角速度誤差を推定する拡張カルマンフィルタを構成し、
前記拡張カルマンフィルタによるクォータニオン推定値を現在の姿勢として出力することを特徴とする小型姿勢センサ。
An angular velocity sensor that measures angular velocities on three independent axes of the moving object;
An acceleration sensor that measures acceleration on three independent axes of the moving object;
A magnetic sensor for measuring geomagnetism on three independent axes of a moving object;
An arithmetic processing unit that calculates an estimated value of a quaternion based on measurement values of the angular velocity sensor, the acceleration sensor, and the magnetic sensor;
The arithmetic processing unit calculates a coordinate transformation matrix using a pre-estimated quaternion estimated value, converts a geomagnetic vector and a gravitational acceleration vector of a ground fixed coordinate system component using the coordinate transformation matrix, A coordinate conversion unit for estimating the geomagnetic vector and the gravitational acceleration vector of the coordinate system component;
A magnetic / acceleration sensor data acquisition unit for acquiring a geomagnetic vector and a gravitational acceleration vector of a sensor coordinate system component of the current step from the magnetic sensor and the acceleration sensor;
The geomagnetic vector and gravitational acceleration vector of the current step sensor coordinate system component estimated by the coordinate conversion unit, and the geomagnetic vector and gravitational acceleration vector of the current step sensor coordinate system component acquired by the magnetic / acceleration sensor data acquisition unit. A vector error calculation unit for calculating an error;
An angular velocity sensor data acquisition unit that acquires an angular velocity vector of a sensor coordinate system component of the current step from the angular velocity sensor;
About the acquired angular velocity data, an angular velocity data error calculating unit for calculating an angular velocity sensor error including a bias error, a scale factor error, a misalignment error, and a noise element;
A quaternion calculation unit and a coordinate conversion unit, a magnetic / acceleration sensor data acquisition unit, a vector error calculation unit, an angular velocity data error calculation unit, and a quaternion calculation unit that obtain time changes between the quaternion of the previous step and the quaternion of the current step Construct a quaternion representing the current posture based on, and an extended Kalman filter that estimates angular velocity error,
A small attitude sensor that outputs a quaternion estimated value by the extended Kalman filter as a current attitude.
前記拡張カルマンフィルタは、
の式および、
低周波数の加速度外乱
dr=[abxbybzT
に対して
のダイナミクスを定義し、式(61)およびこのダイナミクスについて、
を状態量としたときの状態方程式
ただし、前記状態方程式において各行列は次式のようになっている、
を演算することを特徴とする請求項1に記載の小型姿勢センサ。
The extended Kalman filter is
And the formula
Low-frequency acceleration disturbance a dr = [a bx a by a bz ] T
Against
Define the dynamics of Eq. (61) and for this dynamics:
Equation of state where is a state quantity
However, each matrix in the state equation is as follows:
The small posture sensor according to claim 1, wherein:
前記拡張カルマンフィルタは、請求項2の演算に追加して、式(62)の状態方程式を離散化して、
ただし、xt はt ステップ目の状態量であり、Δt はサンプリングタイムである、
の式を求め、ここで、改めて
とおいて、離散時間の状態方程式を、
の形に求め、続いて、式(62)の状態方程式で推定された低周波数の加速度外乱を観測方程式に含めて、
なる観測方程式を求め、
ここで、
として、
ただし、ここでht(xt)は、
なる形の、yt を出力とする観測方程式を求め、この観測方程式により推定値を演算することを特徴とする請求項1記載の小型姿勢センサ。
In addition to the operation of claim 2, the extended Kalman filter discretizes the state equation of equation (62),
However, x t is the state quantity at the t step, Δt is the sampling time,
Where, again,
Where the discrete-time equation of state is
Next, the low-frequency acceleration disturbance estimated by the state equation of Equation (62) is included in the observation equation,
Find the observation equation
here,
As
Where h t (x t ) is
The small attitude sensor according to claim 1, wherein an observation equation having yt as an output is obtained and an estimated value is calculated by the observation equation.
前記拡張カルマンフィルタは、請求項3の演算に追加して、
行列Ft及びHtを、
のように定義し、
以上のように定義された行列と式(64)、式(67)を用いた拡張カルマンフィルタアルゴリズムを以下に示す式で与え、
これらのアルゴリズムにより、演算処理することを特徴とする請求項1に記載の小型姿勢センサ。
The extended Kalman filter is added to the operation of claim 3,
The matrices F t and H t are
Defined as
The extended Kalman filter algorithm using the matrix defined above and Equations (64) and (67) is given by the following equation:
The small-sized attitude sensor according to claim 1, wherein arithmetic processing is performed by these algorithms.
JP2012067533A 2012-03-23 2012-03-23 Small attitude sensor Active JP5061264B1 (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2012067533A JP5061264B1 (en) 2012-03-23 2012-03-23 Small attitude sensor

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2012067533A JP5061264B1 (en) 2012-03-23 2012-03-23 Small attitude sensor

Publications (2)

Publication Number Publication Date
JP5061264B1 JP5061264B1 (en) 2012-10-31
JP2013200162A true JP2013200162A (en) 2013-10-03

Family

ID=47189593

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2012067533A Active JP5061264B1 (en) 2012-03-23 2012-03-23 Small attitude sensor

Country Status (1)

Country Link
JP (1) JP5061264B1 (en)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104569914A (en) * 2014-12-29 2015-04-29 南京航空航天大学 Motion attitude sensing method based on electromagnetic wave polarization three-dimensional electric/magnetic signal stationary characteristics
CN107084743A (en) * 2016-02-12 2017-08-22 通用汽车环球科技运作有限责任公司 Skew and misalignment compensation using the six degree of freedom Inertial Measurement Unit of GNSS/INS data
JP2017526906A (en) * 2014-07-03 2017-09-14 日本テキサス・インスツルメンツ株式会社 Pedestrian navigation device and method
JP2018060326A (en) * 2016-10-04 2018-04-12 株式会社豊田中央研究所 Tracking device and program thereof
CN108534772A (en) * 2018-06-24 2018-09-14 西宁泰里霍利智能科技有限公司 Attitude angle acquisition methods and device
JP2019120587A (en) * 2018-01-05 2019-07-22 ローム株式会社 Positioning system and positioning method
US11361465B2 (en) 2019-04-24 2022-06-14 Canon Kabushiki Kaisha Image capturing apparatus and control method thereof, and orientation angle calculation apparatus for estimating orientation of image capturing apparatus
JP7120361B2 (en) 2021-03-03 2022-08-17 カシオ計算機株式会社 Correction device, attitude calculation device, correction method, and program

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103363993B (en) * 2013-07-06 2016-04-20 西北工业大学 A kind of aircraft angle rate signal reconstructing method based on Unscented kalman filtering
CN104880201B (en) * 2015-03-26 2016-01-13 武汉大学 MEMS gyro automatic calibration method
CN104898681B (en) * 2015-05-04 2017-07-28 浙江工业大学 A kind of quadrotor attitude acquisition method for approximately finishing card quaternary number using three ranks
CN106526644A (en) * 2016-10-10 2017-03-22 上海联适导航技术有限公司 Method of calculating relative angle variation between carrier parts
FR3069634B1 (en) * 2017-07-28 2021-06-11 Sysnav METHOD AND DEVICE FOR CHARACTERIZING A DETERMINED HEADING FROM MEASUREMENT OF THE MAGNETIC FIELD
CN108363298A (en) * 2018-01-17 2018-08-03 合肥工业大学 A kind of quadrotor drone Fast Convergent control method based on quaternion representation
JP6870635B2 (en) * 2018-03-08 2021-05-12 セイコーエプソン株式会社 Inertial measurement units, mobiles, portable electronic devices, and electronic devices
CN108427427B (en) * 2018-03-16 2021-03-26 北京控制工程研究所 Method for calculating attitude angle of spacecraft to earth surface orientation target
CN110132287B (en) * 2019-05-05 2023-05-05 西安电子科技大学 Satellite high-precision joint attitude determination method based on extreme learning machine network compensation
JP7421923B2 (en) * 2019-12-23 2024-01-25 フォルシアクラリオン・エレクトロニクス株式会社 Position estimation device and position estimation method
CN112304330B (en) * 2020-10-29 2024-05-24 腾讯科技(深圳)有限公司 Method for displaying running state of vehicle and electronic equipment
CN113514789B (en) * 2021-04-23 2022-06-07 北京大学 Magnetic sensor array calibration method
CN113268070A (en) * 2021-05-18 2021-08-17 哈尔滨理工大学 Unmanned aerial vehicle attitude angle tracking system based on multisensor
CN114176576B (en) * 2021-12-11 2024-05-24 江苏智恒文化科技有限公司 Method for identifying human motion state based on acceleration

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4830559B2 (en) * 2006-03-16 2011-12-07 セイコーエプソン株式会社 Positioning device and positioning method
JP4726134B2 (en) * 2006-04-03 2011-07-20 国立大学法人 東京大学 MOBILE BODY CONTROL DEVICE AND MOBILE BODY CONTROL METHOD
JP5569681B2 (en) * 2010-04-23 2014-08-13 国立大学法人 東京大学 POSITION ESTIMATION DEVICE AND POSITION ESTIMATION METHOD OF MOBILE BODY USING INDUCTION SENSOR, MAGNETIC SENSOR AND SPEED METER

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017526906A (en) * 2014-07-03 2017-09-14 日本テキサス・インスツルメンツ株式会社 Pedestrian navigation device and method
CN104569914A (en) * 2014-12-29 2015-04-29 南京航空航天大学 Motion attitude sensing method based on electromagnetic wave polarization three-dimensional electric/magnetic signal stationary characteristics
CN107084743A (en) * 2016-02-12 2017-08-22 通用汽车环球科技运作有限责任公司 Skew and misalignment compensation using the six degree of freedom Inertial Measurement Unit of GNSS/INS data
CN107084743B (en) * 2016-02-12 2020-07-07 通用汽车环球科技运作有限责任公司 Offset and misalignment compensation for six degree of freedom inertial measurement units using GNSS/INS data
JP2018060326A (en) * 2016-10-04 2018-04-12 株式会社豊田中央研究所 Tracking device and program thereof
JP2019120587A (en) * 2018-01-05 2019-07-22 ローム株式会社 Positioning system and positioning method
JP7025215B2 (en) 2018-01-05 2022-02-24 ローム株式会社 Positioning system and positioning method
CN108534772A (en) * 2018-06-24 2018-09-14 西宁泰里霍利智能科技有限公司 Attitude angle acquisition methods and device
CN108534772B (en) * 2018-06-24 2021-07-02 西宁泰里霍利智能科技有限公司 Attitude angle acquisition method and device
US11361465B2 (en) 2019-04-24 2022-06-14 Canon Kabushiki Kaisha Image capturing apparatus and control method thereof, and orientation angle calculation apparatus for estimating orientation of image capturing apparatus
JP7120361B2 (en) 2021-03-03 2022-08-17 カシオ計算機株式会社 Correction device, attitude calculation device, correction method, and program

Also Published As

Publication number Publication date
JP5061264B1 (en) 2012-10-31

Similar Documents

Publication Publication Date Title
JP5061264B1 (en) Small attitude sensor
CN106959110B (en) Cloud deck attitude detection method and device
Aydemir et al. Characterization and calibration of MEMS inertial sensors for state and parameter estimation applications
JP4876204B2 (en) Small attitude sensor
Cheviron et al. Robust nonlinear fusion of inertial and visual data for position, velocity and attitude estimation of UAV
JP2012173190A (en) Positioning system and positioning method
JP2014089113A (en) Posture estimation device and program
JP2011227017A (en) Device and method for attitude estimation of moving body using inertial sensor, magnetic sensor, and speed meter
CN106813679B (en) Method and device for estimating attitude of moving object
JPH095104A (en) Method and apparatus for measurement of three-dimensional attitude angle of moving body
JP2015179002A (en) Attitude estimation method, attitude estimation device and program
Allibert et al. Velocity aided attitude estimation for aerial robotic vehicles using latent rotation scaling
Al-Fahoum et al. Design of a modified Madgwick filter for quaternion-based orientation estimation using AHRS
RU2564379C1 (en) Platformless inertial attitude-and-heading reference
Bazhenov et al. Software for the mobile robot spatial orientation system
JP2013122384A (en) Kalman filter and state estimation device
WO2019127092A1 (en) State estimatation
KR20170092359A (en) System for detecting 3-axis position information using 3-dimention rotation motion sensor
JP2013061309A (en) Kalman filter, state estimation device, method for controlling kalman filter, and control program of kalman filter
Haotian et al. Accurate attitude estimation of HB2 standard model based on QNCF in hypersonic wind tunnel test
Li et al. Real-time indoor navigation of uav based on visual delay compensation
Lima et al. Performance evaluation of attitude estimation algorithms in the design of an ahrs for fixed wing uavs
Zhang et al. An improved Kalman filter for attitude determination of multi-rotor UAVs based on low-cost MEMS sensors
JP6516332B2 (en) Mobile terminal, program and method for estimating forward direction of user using angular velocity sensor
McLean et al. A robotic joint sensor

Legal Events

Date Code Title Description
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: 20120710

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20120806

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20150810

Year of fee payment: 3

R150 Certificate of patent or registration of utility model

Ref document number: 5061264

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313117

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

S111 Request for change of ownership or part of ownership

Free format text: JAPANESE INTERMEDIATE CODE: R313113

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350