(57)【要約】【課題】カルマンフィルタを単純に用いて検出値から積分量を推定する場合に比較して、カルマンフィルタを用いた検出値からの積分量の推定精度を向上させることを目的とする。【解決手段】加速度を検出する加速度センサ１２と、カルマンフィルタを用いて観測値から積分量としての変位量を推定する積分量推定プログラムを記憶し、当該プログラムを実行することで、カルマンフィルタの推移行列に予め定めた係数を付加した行列を用いて、加速度センサ１２の検出結果から変位量を推定するＰＣ１８と、を備える。【選択図】図１
2016061742.tif 000002
【特許請求の範囲】
【請求項１】
特定の物理量を検出する検出部と、
カルマンフィルタの推移行列に予め定めた係数を付加した行列を用いて、前記検出部の検出結果から積分量を推定する推定部と、
を備えた積分量推定装置。
【請求項２】
前記係数は、前記推移行列の１行目に少なくとも１つ以上付加されている請求項１に記載の積分量推定装置。
【請求項３】
前記推移行列に直流的な変動を補正する項を更に設けた請求項１又は請求項２に記載の積分量推定装置。
【請求項４】
前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２）
Ｈ＝（１，−２，１）／Δｔ２
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ〜Ｆｃの何れかを用いて２階積分した積分量を推定する請求項１に記載の積分量推定装置。
【数１】
2016061742.tif 000022
【請求項５】
前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２，Ｄｎ）、
Ｈ＝（１／Δｔ２，−２／Δｔ２，１／Δｔ２，１）
Ｄｎ＝（時刻ｎでの予め定めた超低周波数ノイズ）
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ〜Ｆｃの何れかを用いて２階積分した積分量を推定する請求項１に記載の積分量推定装置。
【数２】
2016061742.tif 000023
【請求項６】
前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１）
Ｈ＝（１，−１）／Δｔ
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ又は行列Ｆｃを用いて積分量を推定する請求項１に記載の積分量推定装置。
【数３】
2016061742.tif 000024
【請求項７】
前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，Ｄｎ）、
Ｈ＝（１，−１，Δｔ）／Δｔ
Ｄｎ＝（時刻ｎでの予め定めた超低周波数ノイズ）
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ又は行列Ｆｃを用いて積分量を推定する請求項１に記載の積分量推定装置。
【数４】
2016061742.tif 000025
【請求項８】
コンピュータを、請求項１〜７の何れか１項に記載の積分量推定装置の前記推定部として機能させるための積分量推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、積分量推定装置及び積分量推定プログラムに関する。
【背景技術】
【０００２】
積分量推定装置は、観測値から観測値の積分量を推定する技術が知られている。例えば、例えば、加速度を観測値として検出し、検出した加速度を１回積分することで速度を推定し、２回積分することで変位量を推定することができる。
【０００３】
加速度センサ等によって検出された加速度信号等の観測値には、重力や、重心を中心とする回転に伴う遠心力等、変位には直接関係ない加速度成分、検出器のゼロ点のドリフトなどが超低周波数ノイズとして含まれる。これらの観測値に含まれるノイズは、積分結果に多大な影響を与えるので、ノイズを除去する技術が提案されている（例えば、特許文献１〜３）。
【０００４】
特許文献１では、加速度を２回積分することによって変位相当パラメータを演算する。また、変位相当パラメータに含まれるドリフト成分が単調に増加あるいは減少することに基づき変位相当パラメータの演算結果から変位相当パラメータに含まれるドリフト成分を推定する。そして、変位相当パラメータから推定したドリフト成分を減じて変位信号を演算することが提案されている。
【０００５】
特許文献２では、加速度センサの検出信号に対してフィルタ処理を施して、低周波成分が除去された加速度データを生成する。加速度データを積分して得られた仮の速度データに対する一次回帰直線を求め、該一次回帰直線を用いて仮の速度データを補正することにより、補正速度データを生成する。そして、補正速度データを積分して得られた仮の位置データに対する二次回帰直線を求め、該二次回帰曲線を用いて仮の位置データを補正することにより、補正位置データを生成することが提案されている。
【０００６】
特許文献３には、加速度計測値を積分する際にカルマンフィルタによって算出されたトレンド成分を補正値として用い除去することで１回積分を実行し、加速度記録を速度記録に変換する。そして、算出された速度記録に対して、再度カルマンフィルタによるトレンド成分を算出し、これを除去して再度積分計算を行うことで速度記録を変位記録にすることが提案されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特許２００３−１３０６２９号公報
【特許文献２】特開２００６−１７７７４９号公報
【特許文献３】特許２００６−１９４８２２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明は、カルマンフィルタを単純に用いて検出値から積分量を推定する場合に比較して、カルマンフィルタを用いた検出値からの積分量の推定精度を向上させることを目的とする。
【課題を解決するための手段】
【０００９】
請求項１に記載積分量推定装置は、特定の物理量を検出する検出部と、カルマンフィルタの推移行列に予め定めた係数を付加した行列を用いて、前記検出部の検出結果から積分量を推定する推定部と、を備えている。
【００１０】
請求項２に記載の発明は、請求項１に記載の発明において、前記係数は、前記推移行列の１行目に少なくとも１つ以上付加されている。
【００１１】
請求項３に記載の発明は、請求項１又は請求項２に記載の発明において、前記推移行列に直流的な変動を補正する項を更に設けた。
【００１２】
請求項４に記載の発明は、請求項１に記載の発明において、前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２）
Ｈ＝（１，−２，１）／Δｔ２
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ〜Ｆｃの何れかを用いて２回積分した積分量を推定する。
【００１３】
【数１】
2016061742.tif 000003
【００１４】
請求項５に記載の発明は、請求項１に記載の発明において、前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２，Ｄｎ）、
Ｈ＝（１／Δｔ２，−２／Δｔ２，１／Δｔ２，１）
Ｄｎ＝（時刻ｎでの予め定めた超低周波数ノイズ）
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ〜Ｆｃの何れかを用いて２回積分した積分量を推定する。
【００１５】
【数２】
2016061742.tif 000004
【００１６】
請求項６に記載の発明は、請求項１に記載の発明において、前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１）
Ｈ＝（１，−１）／Δｔ
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ又は行列Ｆｃを用いて積分量を推定する。
【００１７】
【数３】
2016061742.tif 000005
【００１８】
請求項７に記載の発明は、請求項１に記載の発明において、前記推定部は、前記物理量を測定した時間間隔をΔｔ、カルマンフィルタの状態ｘｎ、観測行列Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，Ｄｎ）、
Ｈ＝（１，−１，Δｔ）／Δｔ
Ｄｎ＝（時刻ｎでの予め定めた超低周波数ノイズ）
とし、前記推移行列に前記係数を付加した行列として以下の行列Ｆａ又は行列Ｆｃを用いて積分量を推定する。
【００１９】
【数４】
2016061742.tif 000006
【００２０】
請求項８に記載の積分量推定プログラムは、コンピュータを、請求項１〜７の何れか１項に記載の積分量推定装置の前記推定部として機能させる。
【発明の効果】
【００２１】
請求項１に記載の発明によれば、カルマンフィルタの推移行列そのものを用いて検出値から積分量を推定する場合に比較して、カルマンフィルタを用いた検出値からの積分量の推定精度を向上させた積分量推定装置を提供することができる。
【００２２】
請求項２に記載の発明によれば、カルマンフィルタの推移行列そのものを用いて検出値から積分量を推定する場合に比較して、安定して計算を終了することが可能となる。
【００２３】
請求項３に記載の発明によれば、カルマンフィルタの推移行列そのものを用いて検出値から積分量を推定する場合に比較して、計算の発散を抑制することができる。
【００２４】
請求項４に記載の発明によれば、カルマンフィルタの推移行列そのものを用いる場合に比べて、２回積分の積分量の推定精度を向上することができる。
【００２５】
請求項５に記載の発明によれば、カルマンフィルタの推移行列そのものを用いる場合に比べて、計算の発散を抑制して２回積分の積分量の推定精度を向上することができる。
【００２６】
請求項６に記載の発明によれば、カルマンフィルタの推移行列そのものを用いる場合に比べて、積分量の推定精度を向上することができる。
【００２７】
請求項７に記載の発明によれば、カルマンフィルタの推移行列そのものを用いる場合に比べて、計算の発散を抑制して積分量の推定精度を向上することができる。
【００２８】
請求項８に記載の発明によれば、カルマンフィルタの推移行列そのものを用いて検出値から積分量を推定する場合に比較して、カルマンフィルタを用いた検出値からの積分量の推定精度を向上させた積分量推定プログラムを提供することができる。
【図面の簡単な説明】
【００２９】
【図１】本実施形態に係る積分量推定装置の概略構成を示すブロック図である。
【図２】（ａ）は画像形成装置の感光体のクリーニングブレードに加速度センサを設けて、トナー処理量を極端に減らしたストレス条件下で得られた加速度の測定結果をに示す図であり、（ｂ）は（ａ）の時間軸方向の拡大図であり、（ｃ）はフーリエ変換の結果を示す図であり、（ｄ）は平滑化処理後の信号を示す図である。
【図３】（ａ）は検出した加速度を時間に対して二重積分した結果を示す図であり、（ｂ）はベースラインを時間の二次関数で近似して求めて（ａ）から減じた結果を示す図であり、（ｃ）は図２（ｂ）を平滑化処理してベースライン処理した結果を示す図であり、（ｄ）は（ｃ）の拡大図である。
【図４】（ａ）は人工データの高周波数振動成分を示す図であり、（ｂ）は人工データの低周波数振動成分を示す図であり、（ｃ）は人工データのドリフト成分を示す図であり、（ｄ）は人工データの全加速度の拡大図であり、（ｅ）は人工データの二回積分後を示す図である。
【図５】（ａ）は人工データの全変位量を示す図であり、（ｂ）は人工データの高周波数成分を示す図であり、（ｃ）は人工データの低周波数成分を示す図であり、（ｄ）はカルマンフィルタによる全変位量の推定値を示す図であり、（ｅ）はカルマンフィルタによる高周波数成分の推定値を示す図であり、（ｆ）はカルマンフィルタによる低周波数成分の推定値を示す図である。
【図６】（ａ）はＦ１ａで推定計算を行って得られた条件毎の全変位量のｒｍｓ値の比を示す図であり、（ｂ）はＦ１ａで推定計算を行って得られた条件毎の高周波数成分のｒｍｓ値の比を示す図であり、（ｃ）Ｆ１ａで推定計算を行って得られた条件毎の低周波数成分のｒｍｓ値の比を示す図であり、（ｄ）はＦ１ｂで推定計算を行って得られた条件毎の全変位量のｒｍｓ値の比を示す図であり、（ｅ）はＦ１ｂで推定計算を行って得られた条件毎の高周波数成分のｒｍｓ値の比を示す図であり、（ｆ）Ｆ１ｂで推定計算を行って得られた条件毎の低周波数成分のｒｍｓ値の比を示す図である。
【図７】（ａ）はドリフト項有無の比較結果を示す図であり、（ｂ）はドリフト項導入時の各成分の再現性を示す図である。
【図８】ドリフト成分を含む場合の再現性を示す図である。
【図９】（ａ）は人工データの全変位量の再現性を示す図であり、（ｂ）は人工データの高周波数成分の再現性を示す図であり、（ｃ）は人工データの低周波数成分の再現性を示す図であり、（ｄ）は人工データのドリフト成分の再現性を示す図である。
【図１０】本実施形態の推定計算を実測データに適用した結果を示す図であり、（ａ）は高周波数成分（ストレス時）を示す図であり、（ｂ）は高周波数成分（ストレス時）の拡大図であり、（Ｃ）は低周波数成分（ストレス時）を示す図であり、（ｄ）はドリフト成分（ストレス時）を示す図であり、（ｅ）は高周波数成分（正常時）を示す図であり、（ｆ）は高周波数成分（正常時）の拡大図であり、（ｇ）は低周波数成分（正常時）を示す図であり、（ｈ）はドリフト成分（正常時）を示す図である。
【発明を実施するための形態】
【００３０】
以下、本実施の形態について図面を参照しつつ説明する。なお、以下では、積分量推定装置として、加速度から変位量を推定する積分量推定装置を一例として説明する。図１は、本実施形態に係る積分量推定装置の概略構成を示すブロック図である。
【００３１】
積分量推定装置１０は、加速度センサ１２、センサアンプ１４、データロガー１６、及びパーソナルコンピュータ（ＰＣ）１８を備えている。
【００３２】
加速度センサ１２によって物理量として加速度を検出して、加速度検出信号をセンサアンプ１４へ出力する。
【００３３】
センサアンプ１４は、加速度センサ１２によって検出された加速度検出信号を増幅してデータロガー１６へ出力する。
【００３４】
データロガー１６は、加速度センサ１２によって検出されてセンサアンプ１４によって増幅された加速度検出信号を順次記憶する。
【００３５】
ＰＣ１８は、データロガー１６によって記憶された加速度の検出値を観測値として、観測値を２回積分することにより変位量を推定する。具体的には、ＰＣ１８は、カルマンフィルタを用いて観測値から積分量としての変位量を推定する積分量推定プログラムを記憶し、ＰＣ１８が当該プログラムを実行することで、加速度から変位量を推定する。このとき、本実施形態では、カルマンフィルタの推移行列として、予め定めた係数を付加した行列を用いて変位量を推定するようになっている。なお、推移行列に付加する係数は、例えば、カルマンフィルタの推移行列の１行目に少なくとも１つ以上設けた行列を使用する。
【００３６】
さらに、本実施形態では、直流的な変動を補正する項を推移行列に付加することにより、ドリフト成分を考慮して変位量を推定するようになっている。
【００３７】
以下では、ＰＣ１８が加速度センサ１２によって検出された加速度から変位量を推定する際に用いる推移行列に係数を付加したカルマンフィルタの導出について詳細に説明する。
【００３８】
まず、カルマンフィルタの概要を説明する。観測可能な時系列的な現象について、現象の状態を特徴づける量をｘｎ（ｎ＝１，２，・・・、Ｎ）、その特徴量を測定した結果をｙｎとすれば、経時による事象の変化は、以下の（１）、（２）となる。
【００３９】
ｘｎ＝Ｆｘｎ−１＋Ｇｖ  ・・・（１）
ｙｎ＝Ｈｘｎ−１＋ω    ・・・（２）
【００４０】
ここで、太字はベクトルまたは行列を表し、（１）式は、状態ｘｎ−１がｘｎへ変化する過程を、（２）式は現象を特徴づける量が観測されるメカニズム（センサアンプのゲインなど）を表す物である。それらのメカニズムは各々二つの行列、Ｆ（推移行列）、Ｈ（観測行列）で数学的に表現可能でなければならない。各式の第二項（ｖ、ω）は現象そのものや観測の過程で含まれるであろうノイズ（雑音）であり、平均が０で分散がＱ、Ｒの白色ノイズが含まれるものとする。Ｇは計算式を行列・ベクトル表記するために導入された意味合いが強く、システムノイズが適切な状態を表す要素に加わるように決められる。
【００４１】
ここで、状態ｘｎを観測値ｙｊから推定する場合を考える。測定結果にはノイズが含まれ、解析結果には大きなばらつきがあると予測される場合には、推定される状態ｘｎにもばらつきが含まれるはずである。そこで、時間１からｎのうちの特定の時刻ｊに得られた観測ベクトルｙｊからの推定値を、以下のように定義する。
【００４２】
ｘｎ｜ｊ＝Ｅ［ｘｎ｜ｙｊ］ （E[f]はfの平均値を表す）・・・（３）
Ｖｎ｜ｊ＝Ｅ［（ｘｎ−ｘｎ｜ｊ）２｜ｙｊ］            ・・・（４）
【００４３】
すなわち、状態ｘｎ｜ｊを推定値の平均値と分散で代表される「分布を有する量」と表す。各時刻における推定結果一つ一つを、平均値と分散で代表される「曖昧さを許したもの」とする点が、カルマンフィルタによる推定方法の特徴の一つである。
【００４４】
上記式より、観測値ｙｎが得られる前に推定される状態ベクトルｘｎ｜ｎ−１は、一期前の推定値から以下のように予測される。
【００４５】
［一期先予測］
【００４６】
ｘｎ｜ｎ−１＝Ｆｘｎ−１｜ｎ−１      ・・・（５）          
Ｖｎ｜ｎ−１＝ＦＶｎ−１｜ｎ−１Ｆｔ＋ＧＱＧｔ  ・・・（６）
（式中のtは転置演算子）                                
【００４７】
状態ベクトルが平均と分散を有する確率変数を要素とするため、平均、分散の両方が計算される。
【００４８】
この予測された平均値に、分散Ｖｎ｜ｎ−１と新たに観測されたｙｎを用いて修正を加えることで、より確からしい推定値ｘｎ｜ｎを得るのが次のステップである。
【００４９】
［フィルタ］
【００５０】
Ｋｎ＝Ｖｎ｜ｎ−１Ｈｔ（ＨＶｎ｜ｎ−１Ｈｔ＋Ｒ）−１    ・・・（７）
ｘｎ｜ｎ＝ｘｎ｜ｎ−１＋Ｋｎ（ｙｎ−Ｈｘｎ｜ｎ−１）  ・・・（８）
Ｖｎ｜ｎ＝（Ｉ−ＫｎＨ）Ｖｎ｜ｎ−１              ・・・（９）
【００５１】
（７）式で、事前に予測される不確かさ（Ｖｎ｜ｎ−１、Ｒ）から補正率（Ｋｎ：カルマンゲイン）を求め、観測値を基に量を決めて補正を加えているのが（８）式、（９）式である。補正量を毎回計算しなおすため、多数回の繰り返し計算を行っても誤差が蓄積し難くなっているのも、この方法の特徴といえる。
【００５２】
さらに、全ての観測値と推定値が得られた時点（時刻Ｎ）であれば、以下のようにしてより高精度に推定される。
【００５３】
［固定区間平滑化］
【００５４】
Ａｎ＝Ｖｎ｜ｎＦｔＶ−１ｎ＋１｜ｎ     ・・・（１０）                 
ｘｎ｜Ｎ＝ｘｎ｜ｎ＋Ａｎ（ｘｎ＋１｜Ｎ−ｘｎ＋１｜ｎ）  ・・・（１１）     
Ｖｎ｜Ｎ＝Ｖｎ｜ｎ＋Ａｎ（Ｖｎ＋１｜Ｎ−Ｖｎ＋１｜ｎ）Ａｔｎ  ・・・（１２）
【００５５】
最後に推定された値から逆順に遡って計算が進むため、特に、初期値が不明な場合に有効である。
【００５６】
次に、加速度を観測値としてカルマンフィルタを用いた変位量の推定への応用を説明する。
【００５７】
加速度は速度の微分量であり、また、位置の二次微分量であるから、時間に対する離散値の差分方程式として以下のように記述される。
【００５８】
ａｎ＝（ｖ＋ｎ−ｖ−ｎ）／Δｔ
＝｛（ｐｎ＋１−ｐｎ）／Δｔ−（ｐｎ−ｐｎ−１）／Δｔ｝Δｔ
＝（ｐｎ＋１−２Ｐｎ＋Ｐｎ−１）／Δｔ２＝観測量（ｙｎ）  ・・・（１３）
【００５９】
ここで、ａｎは時刻ｎにおける加速度、ｐｎは位置、Δｔはサンプリングの時間間隔である。ｖは速度であり、ｖ＋ｎは位置ｐｎからｐｎ＋１へ移動する間の平均速度、ｖ−ｎは位置ｐｎ−１からｐｎへ移動する間の平均速度である。よって、観測量である加速度は、時刻（ｎ＋１、ｎ、ｎ−１）における３点の位置情報から求められる。しかし、このままでは未来の位置で現在を推定することになってしまうため、サンプリング間隔が振動周波数に対して十分早ければ、以下の式が成り立つと仮定する。
【００６０】
観測量（ｙｎ）＝ａｎ≒（ｐｎ−２ｐｎ−１＋ｐｎ−２）／Δｔ２  ・・・（１４）
よって、推定すべき状態ベクトルｘｎと観測量ベクトルｙｎは、
ｙｎ＝（ａｎ），ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２）  ・・・（１５）
と表され、（２）式より、
【００６１】
【数５】
2016061742.tif 000007
【００６２】
が得られる。
【００６３】
Ｆについては、（１）式の統合が厳密に成り立つ必要がないため様々な求め方で導出される。今回は、サンプリング間隔が振動周波数に対して十分に速いとして求められる（１７）式と、そのような仮定が必要ない（１８）式について推定結果を比較することにした。
【００６４】
【数６】
2016061742.tif 000008
【００６５】
ただし、（７）式中の逆行列の計算に失敗してＫを求めることができないことが多々あったことから、安定した計算が可能なように、以下のように、推移行列に係数δを付加している。
【００６６】
【数７】
2016061742.tif 000009
【００６７】
なお、δは、結果に擾乱を与えない小さな値であり、経験的に決定したものである。係数δについては詳細を後述する。
【００６８】
ところで、定義式中に含まれるノイズ項は「白色ノイズｖ、ω」だけであるため、直流的なドリフト成分は別に導入する必要がある。変化の挙動が大きくなる複数成分を同時に取り扱う方法としては、景気変動の時系列的解析で用いられる「多項式トレンドモデル」が知られている。「多項式トレンドモデル」では、多数年にわたるトレンド変化の成分と季節変動による周期的な変動を以下の式として取り扱う。
【００６９】
【数８】
2016061742.tif 000010
【００７０】
そこで、以下のようにして各小行列を求め、その適用可否を検証した。
【００７１】
まず、観測される加速度は振動ａｎ、ドリフトＤｎの二つの成分の和であるので、以下の式で表される。
【００７２】
ｙｎ＝（ａｎ＋Ｄｎ）  ・・・（２１）    
【００７３】
よって、（１６）式を求めたのと同様にして、以下の式が得られる。
【００７４】
【数９】
2016061742.tif 000011
【００７５】
また、状態ベクトルｘｎは確率変数的な分布を有するものであり、ノイズや分布の範囲を含めた上で（１）式、（２）式が成り立てば十分であるから、以下の式が得られる。
【００７６】
Ｄｎ＝Ｆドリフト＋Ｄｎー１＋ｖ、Ｄｎ＝（Ｄｎ）、Ｆドリフト＝（１）・・・（２４）
【００７７】
よって、（２４）式、（１９）式、及び（２０）式から以下の式が得られる。
【００７８】
【数１０】
2016061742.tif 000012
【００７９】
なお、Ｆ１ｃは、ｘｎとｘｎ−１との差が誤差程度の差しかない場合から求めたものである。
【００８０】
最後に、行列Ｒは観測値が１種類のままであるから（１６）式と同様の要素が１個の行列となり、Ｇ、Ｑについては状態ベクトルが（２３）式の形になったことから、以下のように導かれる。
【００８１】
【数１１】
2016061742.tif 000013
【００８２】
なお、実際の計算では、初期条件には無情報事前分布
ｘ０｜０＝０，Ｖ０｜０＝１×１０８    ・・・（２６）
を適用し、システムノイズ、観測ノイズには計算の柔軟性を損なわないよう大きめの値、
Ｒ＝１，Ｑ振動＝１，Ｑドリフト＝０．１      ・・・（２７）
を、全ての計算で統一して使用した。
【００８３】
ここで、比較例として積分法を用いて加速度信号から変位量を求める場合について考察する。
【００８４】
画像形成装置の感光体のクリーニングブレードに加速度センサ１２を設けて、トナー処理量を極端に減らしたストレス条件下で得られた加速度の測定結果を図２（ａ）に示す。測定は、２０秒間（Ａ３用紙２０枚相当）にわたって行ったが、生データから伺える振動の状態は単調なものである。時間軸方向へ拡大すると（図２（ｂ））、数ミリ秒周期で振動していることが分かり、フーリエ変換の結果（図２（ｃ））からその周波数は２００Ｈｚと求められ、金属支持部材の固有振動数と一致する。さらに、平滑化処理（ローパスフィルタ処理）によってフーリエ変換では検出が難しい、低繰り返し性の信号を抽出したところ、用紙周期に相当する時間間隔で単発的な上下動を示すと思われる信号（図２（ｄ））が抽出された。金属支持部材の固有振動数である２００Ｈｚの振動からゴムブレードの状態変化や運動に関する情報を得ることは困難と思われるため、この図２（ｄ）の信号の発生源（動き）について、さらに検討する。
【００８５】
支持部材の動きを把握することを目的として、時間に対して二重積分した結果を図３（ａ）に示すが、わずか２０秒間で６ｍほど沈み込む結果となってしまった。そこで、ベースラインを時間の二次関数（-0.002236+0.002582t-0.015304t2）に近似して求め、図３（ａ）から減じて図３（ｂ）を得た。重心の移動がない振動運動の加速度データの平均は０となるはずなので、積分定数を適切に処理できれば二重積分後も０のままであるはずだが、
（１）直流的なドリフト（原点からのずれ＝-0.015304m-2）
が加わったため、二次関数的なベースラインの曲がりが生じたものと考えられる。
ところで、ベースライン処理を施した図３（ｂ）からは、
（２）５秒を越える時間領域での±１〜２ｍｍ程度のゆっくりした大きな変位
（３）約１秒おきに発生しているスパイク状（下向き数百μｍ）の変位
の存在が推測される。また、拡大図（図２（ｂ））を平滑化処理してベースライン処理する（ハイパスフィルタ処理）することで図３（ｃ）が得られる。その結果からは、
（４）±３μｍ程度のミリ秒周期での微小変位変動（２００Ｈｚ振動に相当）
（５）±５〜１０μｍ程度の約１秒周期での単発的変位変動
の存在が推測される。ただし、（５）に相当する信号は１００ｍｓｅｃ周期とゆっくりしたものであり、本来は１秒周期で単発的な（３）と同様の現象を示すものとして検討を加えるのが妥当である。平滑化で求めたベースラインが柔軟性に欠けたため、信号の変化に追随できなかったところに問題が有ると考えられる。なお、（５）を除いた領域での２００Ｈｚの振動の大きさは、二乗平均平方根値（以下、ｒｍｓ値と略記する）は２．８μｍである。
【００８６】
次に、比較例の積分法を考慮して本実施形態における上述の係数を付加した推移行列を用いたカルマンフィルタによる状態推定を検討する。
【００８７】
変位量（周波数、振幅量）が既知の加速度データを人工的に作成し、モデル計算の結果から推移行列Ｆの具体的な形を決定した後に、実測データへの適用を試みた。
【００８８】
上記の積分法を用いた考察から、以下の成分を足し合わせることにした。
【００８９】
（Ａ）「振動」等の周期的な変位の二階微分相当量（（３）、（４）、（５））
【００９０】
【数１２】
2016061742.tif 000014
【００９１】
（Ｂ）ＤＣ〜低周波数の「ベースラインのうねり」（（１）、（２））
【００９２】
Ｄ（ｔ）＝Ｄ０（１＋ｓｉｎ（２π（０．０５）ｔ）    ・・・（２９）
Ｄ０：ベースライン変動量の平均値［ｍｓ−２］                        
【００９３】
ｒ１，ｆ１は図３（ｄ）に示した高周波数成分に対応するものである。ｒ２、ｆ２は図３（ｃ）で見られる１秒周期で発生している下向きのスパイク状の変位に対応するもので、正弦関数を３乗して波形を先鋭化した。ドリフト成分としては、図３（ｂ）から１０秒を越える周期でのゆっくりした変位変動から発生する加速度成分と、図２（ａ）から予測される直流的な原点シフトの重ね合わせを仮定した。図３から得られる各パラメータの値（ｒ１＝３μｍ、ｆ１＝２００Ｈｚ、ｒ２＝１ｍｍ、ｆ２＝１Ｈｚ、Ｄ０＝−０．０１ｍｓ−１）を用いて計算した信号の形を、図４に示す。（２８）式からえられた値（図（ａ）＋（ｂ））を数値的に２階微分した後、（２９）式で得られた値（図４（ｃ））を加えて人工の加速度データ（図４（ｄ））とした。このデータを２重積分した結果（図４（ｅ））は、図３（ａ）と同様、ベースラインが大きく曲がる。
【００９４】
続いて、部分行列Ｆ１を選定する。
まず、ドリフト成分を加えない人工データを用いてＦ１を決定する。（２４）式で示したように、決定のためには基本となる時系列方程式の形と共に、計算の異常終了を防ぐための経験的なパラメータ（δ）の決定が必要である。
【００９５】
そこで、図４に示した各パラメータが含まれるように大きく値を振った各二水準の直積を条件とした人工データを作成し、かつ、δを０．１、０．０１、０．００１と変えながら推定計算を行って再現性を比較した。設定した具体的な値を条件番号と共に表１に、計算例として、推移行列にＦ１ａ、δ＝０．０１とした時の計算結果例を図５に例示した。推定値は分布の平均値ｘｎと分散で表されているが、表示を簡略化するためにｘｎだけを示した。
【００９６】
【表１】
2016061742.tif 000015
図５（ａ）は、高周波数成分（図５（ｂ））と低周波数成分（図５（ｃ））を足し合わせて作成された全変位量を示す波形であり、正弦関数の３乗である先鋭化された低周波数変動に高周波数成分が重畳して線幅が広くなったような波形が得られている。この波形を２階微分したデータにカルマンフィルタを適用して得られた推定値ｘｎを、図５（ｄ）に示す。図５（ｄ）はフィルタで得られた結果と、固定区間平滑化で得られた結果を示したものである。図５（ｄ）中のハッチングで示す部分が固定区間平滑化で得られた結果に対してずれているフィルタで得られた結果の部分を示す。フィルタのｔ＝０〜２秒の領域で大きなずれが見られたものの、それ以外の部分では図５（ａ）を再現することができた。また、固定区間平滑化後のデータにローパス処理、ハイパス処理を加えて得られた結果（図５（ｅ），（ｆ））も、人工データの信号波形を良好に再現することができている。高周波数成分の推定値をフーリエ変換して求めた周波数は１００Ｈｚと人工データと一致し、各変位量のｒｍｓ値も、表２に示すように一致を示している。
【００９７】
【表２】
2016061742.tif 000016
このようにして全ての条件（表３）について推定計算を行い、得られたｒｍｓ値の比を図６に示す。なお、高周波数成分の周波数は全て１００Ｈｚであった。また、図６では、「▲」はδ＝０．００１、「●」はδ＝０．０１、「＋」はδ＝０．１を示す。
【００９８】
【表３】
2016061742.tif 000017
【００９９】
図６によれば、δ＝０．０１とした場合は、Ｆ１ａ、Ｆ１ｂともに全条件で良好な結果を得ることができた。特に、高周波数成分では全ての比がほぼ１となり、図中の点が重なって表示されているほどである。しかし、δ＝０．１では低周波数成分の周波数が０．５Ｈｚの水準での再現性が悪く、Ｆ１ｂとδ＝０．００１の組み合わせでは、高周波数成分、低周波数成分ともに大きく外れる場合が見られている。このような挙動を示す原因を詳細に把握できていないが、δ＝０．００１では行列計算の発散を抑えきれず、δ＝０．１では理想的な推移方程式からの逸脱が無視できなくなるためと考えられる。
【０１００】
以上より、Ｆ１ａとδ＝０．０１の組み合わせが安定して良い再現性を得られることが分かった。
【０１０１】
続いて、ドリフト項の導入効果について説明する。
表３で示した条件のうち、信号強度が最も小さくなる条件である１番のパラメータを選択し、これに想定される値の１／１０〜１０倍ドリフト成分を加えた場合に得られる全変位推定量の再現性を図７（ａ）に示す。なお、推移行列はＦ１ａ，δ＝０．０１を使用している。図７（ａ）の「▲」で示す点線はドリフト項を考慮していない場合に得られた結果であるが、δ＝０．００１で既に計算が発散してしまう。これに対して、ドリフト項を考慮した手法で計算した結果（図７（ｂ）の「○」）は、想定した値の２倍（Ｄ０＝−００２）付近まで、ほぼ人工データを再現することができている。また、このときの他の成分の再現性を図７（ｂ）に示したが、低周波数成分の再現性に困難が有るものの、高周波数成分およびドリフト成分については、全領域でほぼ満足される結果が得られた。
【０１０２】
次に、Ｄ０を想定される−０．０１に固定し、他のパラメータを表３と同様に変えて推定計算による再現性を検討した結果を図８に示す。高周波数成分とドリフト成分はほぼ１と良好な再現性を示しており、低周波数成分の再現が困難な場合に全変位量の再現性も悪化していることが分かる。また、１から大きく外れているのは低周波数成分の周波数が０．５Ｈｚの水準（条件番号１Ｄ、２Ｄ、５Ｄ、６Ｄ、９Ｄ、１０Ｄ、１３Ｄ、１４Ｄ）に多く見られる。この原因については、今回導出した計算方法では、一個の加速度センサ１２から得られる一次元の情報、時間変化に伴う変動量だけで二成分に分離しようとした点に有ると考えられる。例えば、条件１Ｄでの計算で得られた共分散行列では、二つの成分が完全に独立に取り扱われていれば「０」になるべき非対角項は約７となっており、二つの成分が混交した結果が得られていることが予測される。
【０１０３】
よって、実際の使用時には、低周波数成分の振動やドリフト量の絶対値については、それほど高い信頼性は得られていないものと考える必要がある。
【０１０４】
これまで検討してきた計算方法を、今回の事例で想定したモデルデータ（図４）に適用した例を図９（ａ）〜（ｄ）に示す。原理的に不一致が大きくなる時間領域の両端で、人工データ（点線）との僅かな不一致が見られるものの、全体的な再現性は非常に高い。よって、実測データへ適用するのに十分な性能を有する。
【０１０５】
続いて、測定データへ適用した場合について説明する。
【０１０６】
今回検討した推定計算を、実測データである図２（ａ）に適用した結果を図１０（ａ）〜（ｄ）に示す。なお、全変位量は低周波数成分の図とほぼ一致することから省略した。
【０１０７】
何れの推定結果にもベースラインの曲がりは見られず、筐体やフレームに強固に固定された部材の運動を推定した結果として妥当な結果が得られている。この結果を、上述した積分法による解析結果と比較すると、以下のような違いが見られる。
【０１０８】
（１）直流的なドリフト（原点からのずれ）
ドリフト成分はほぼ一定の値となり、その平均値は−０．０３６であった。この値は、完全な直流的ドリフトを仮定した全データの平均（−０．０３０４６８）に近く、良好な整合性を表している。ベースラインを二次式に近似して得られる値は二次の項の係数である−０．０１５３０４であり、この近似は不適当であったと思われる。
（２）５秒を越える時間領域での±１〜２ｍｍ程度の大きな変位
このような変位変動はカルマンフィルタによる推定結果には全く見られていない。ベースラインの補正方法が不適切だったために現れたゴーストと考えられる。
（３）約１秒おきに発生しているスパイク状（下向き数百μｍ）の変位
図１０（ｃ）に明瞭に現れた。ベースラインの曲がりが消えたことで、変位量も含め、高い周期性が確認される。
（４）±３μｍ程度のミリ秒周期での微小変位変動
（５）±５〜１０μｍ程度の約１秒周期での単発的変位変動
高周波数成分の推定結果（図１０（ａ），（ｂ））には、フーリエ変換して求めた周波数が２００Ｈｚ、振幅（ｒｍｓ値）が２．７２μｍの波形が得られ、積分法で得られた（４）と一致するものと考えられる。しかし、（５）に相当する波形は全く見られない。これは、（３）の信号を正確に抽出できたため、図中から消失したと考えられる。
【０１０９】
総じて言えば、カルマンフィルタによる状態推定法を用いることで、筐体やフレームに強固に固定された部材の運動を推定した結果として妥当な結果が得られる。
【０１１０】
以上説明したように、加速度測定による振動解析を行うに際して、推移行列に予め定めた係数を付加したカルマンフィルタを用いて積分量を推定することで、積分法で見られたベースライン補正等の困難を回避して、安定して加速度から変位量が推定される。
【０１１１】
なお、上記では、カルマンフィルタの推移行列Ｆの１つ項に係数δを付加した例を説明したが、これに限るものではなく、推移行列Ｆの複数の項に係数δを付加するようにしてもよい。
【０１１２】
すなわち、ｙｎを測定された加速度、ｖを運動そのものが有する雑音（変動）、ωを測定時に加えられる雑音、Ｆを推移行列、Ｈを観測行列として、上述の（１）式、（２）式を用いて、一期先予測（（５）式、（６）式）、及びフィルタ（（７）式〜（９）式）の計算を０からｎまで繰り返し計算して求める方法において、以下の条件（ｘｎ、Ｈ、Ｆ、Ｇ、Ｑ）を用いるようにしてもよい。なお、一期先予測およびフィルタは、ｘｉ｜ｊを時刻ｊまでの推定値から推定した時刻ｊでの値（ばらつきの値）、Ｖをばらつきの分散を表す行列（分散共分散行列）、Ｋをカルマンフィルタ（ベクトル）、Ｒを観測ノイズの分散とする。
【０１１３】
まず、ドリフト項を考慮しない場合（超低周波数成分が小さく、ｖ、ωで表される場合）には、ｘｎ、Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２）
Ｈ＝（１，−２，１）／Δｔ２
Δｔ＝（加速度を測定した時間間隔）
とし、推移行列Ｆに以下の何れかを用いる。この場合には、後述のドリフト項を考慮する場合よりも行列が小さく計算負荷が抑制される。
【０１１４】
【数１３】
2016061742.tif 000018
【０１１５】
一方、ドリフト項を考慮する場合（予め定めた超低周波数成分が大きく、ドリフト項を考慮しないと計算が上手く行かない場合）には、ｘｎ、Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１，ｘｎ−２，Ｄｎ）、
Ｈ＝（１／Δｔ２，−２／Δｔ２，１／Δｔ２，１）
Ｄｎ＝（時刻ｎでの予め定めた超低周波数ノイズ）
とし、推移行列Ｆに以下の何れかを用いる。この場合には、超低周波数成分の推定値も同時に求められるため、静電容量方式のセンサ等を用いて、傾きの変化を同時に求めることも可能となる。
【０１１６】
【数１４】
2016061742.tif 000019
【０１１７】
さらに、上記では、加速度から積分量として変位量を推定する例を説明したが、これに限るものではなく、速度を積分量として推定するようにしてもよい。或いは、加速度以外の物理量から積分量を推定する場合に推移行列に係数δを付加したカルマンフィルタを用いるようにしてもよい。
【０１１８】
例えば、速度を積分量として推定する際には、上述の（１）式、（２）式を用いて、一期先予測（（５）式、（６）式）、及びフィルタ（（７）式〜（９）式）の計算を０からｎまで繰り返し計算して求める。
【０１１９】
ここで、ドリフト項を考慮しない場合（予め定めた超低周波数成分が小さく、ｖ、ωで表される場合）には、ｘｎ、Ｈを、
ｘｎ＝（ｘｎ，ｘｎ−１）
Ｈ＝（１，−１）／Δｔ
Δｔ＝（加速度を測定した時間間隔）
とし、推移行列Ｆに以下の何れかを用いればよい。この場合には、以下のドリフト項を考慮する場合よりも行列が小さく計算負荷が抑制される。
【０１２０】
【数１５】
2016061742.tif 000020
【０１２１】
一方、ドリフト項を考慮する場合には、ｘｎ、Ｈを
ｘｎ＝（ｘｎ，ｘｎ−１，Ｄｎ）、
Ｈ＝（１，−１，Δｔ）／Δｔ
Ｄｎ＝（時刻ｎでの予め定めた超低周波数ノイズ）
とし、推移行列Ｆに以下の何れかを用いる。
【０１２２】
【数１６】
2016061742.tif 000021
【０１２３】
なお、推移行列Ｆの複数の項に係数δを付加する各々の場合における行列やベクトル（ｘｎ、Ｈ、Ｆ）は、共に固定区間平滑化（（１０）式〜（１２）式）や、拡大状態ベクトル法等の拡張計算に適用するようにしてもよい。また、積分定数に対応するｔ＝０でのｘ０＝（０，０，０）としておけば、繰り返し計算を数回から数十回行った段階で妥当な推定値に収束する。また、固定区間平滑化後の推定値は、ｘ０＝（０，０，０）近傍の値も正確に推定可能である。
【０１２４】
また、計算に必要となる各ノイズの分散値は、実験データをもとに最尤推定しても良いし、観測ノイズの分散を１とし、システムノイズを１〜０．００１程度にしても良い。すなわち、推定変位量の平均値は各ノイズの分散値に鈍感であるため、０や極端に大きい値（例えば、１００以上）でなければ良好な結果が得られる。
【０１２５】
また、上記の実施の形態では、積分量推定プログラムをＰＣ１８に実行することにより、加速度から変位量を推定する例を説明したが、積分量推定プログラムによって実行される処理の一部または全てをハードウエアで行うようにしてもよい。
【０１２６】
また、上記の実施形態に係る積分量推定装置１０によって行われる積分量の推定は、プログラムとして記憶媒体に記憶して流通させるようにしてもよい。
【符号の説明】
【０１２７】
１０    積分量推定装置
１２    加速度センサ
１８    ＰＣ
【図１】
2016061742.tif 000026
【図２】
2016061742.tif 000027
【図３】
2016061742.tif 000028
【図４】
2016061742.tif 000029
【図５】
2016061742.tif 000030
【図６】
2016061742.tif 000031
【図７】
2016061742.tif 000032
【図８】
2016061742.tif 000033
【図９】
2016061742.tif 000034
【図１０】
2016061742.tif 000035