(57)【特許請求の範囲】
【請求項１】
振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、該状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置であって、
前記外部入力として振動に関する物理量の外部測定器により測定された測定値を取得する、測定値取得部と、
前記振動センサのセンサ値を取得する、センサ値取得部と、
前回の推定により得られた状態ベクトル推定値と、前回の演算により得られた事後共分散行列と、該前回の推定及び該前回の演算に対応する時刻に対応する前記測定値と、システムノイズの共分散行列とを格納する、データ格納部と、
前記前回の推定により得られた状態ベクトル推定値と、前記前回の推定及び前記前回の演算に対応する時刻に対応する前記測定値とを用いて、状態ベクトル事前推定値を演算する、状態ベクトル事前推定値演算部と、
前記前回の演算により得られた事後共分散行列と、前記システムノイズの共分散行列とを用いて、事前共分散行列を演算する、事前共分散行列演算部と、
前記事前共分散行列を用いてカルマンゲインを演算する、カルマンゲイン演算部と、
前記状態ベクトル事前推定値と、前記カルマンゲインと、前記センサ値とを用いて、今回の推定における状態ベクトル推定値を演算する、状態ベクトル推定値演算部と、
前記事前共分散行列と、前記カルマンゲインとを用いて、今回の演算における事後共分散行列を演算する、事後共分散行列演算部と
を備えた、推定装置。
【請求項２】
前記振動センサと、
前記外部測定器と、
請求項１に記載の推定装置と
を備えた、振動センサシステム。
【請求項３】
前記外部測定器は演算により加速度変換できるもので、前記振動センサのダイナミックレンジの上限を超えるダイナミックレンジを有し、前記測定値取得部は加速度の測定値を取得する、請求項２に記載の振動センサシステム。
【請求項４】
前記振動センサは振り子を備え、変位、速度、加速度のうち１以上の値を測定し、前記センサ値取得部は、該変位、該速度、該加速度のうち１以上の値を取得する、請求項２又は３に記載の振動センサシステム。
【請求項５】
前記カルマンゲイン演算部は、カルマンゲイン調整項を用いた演算によりカルマンゲインを調整し、該カルマンゲインの調整は、前記振動センサの変位、速度、加速度のうち１以上の値が所定値よりも大きい場合には、該１以上の値が所定値よりも小さい場合に比べて該カルマンゲイン調整項を大きくすることでカルマンゲインを小さくすることにより行われる、請求項４に記載の振動センサシステム。
【請求項６】
前記振動センサは、
振り子を備え、変位、速度、加速度のうち１以上の値を測定する振動センサであり、
振り子を備えた振動センサの動作を演算により模擬し、模擬的な振り子の変位、速度、加速度のうち１以上の値を演算により決定する模擬的センサが後段のカルマンフィルタにより備えられ、と
前記状態ベクトル推定値演算部は、前記振動センサにより測定された変位、速度、加速度のうち１以上の値の大きさに応じて、該振動センサにより測定された変位、速度、加速度のうち１以上の値と、前記模擬的センサにより決定された前記模擬的な振り子の変位、速度、加速度のうち１以上の値に係数を乗じたものとのうちいずれかを前記センサ値として用いる、請求項２又は３に記載のシステム。
【請求項７】
振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、該状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置により実行される方法であって、
前記外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、
前記振動センサのセンサ値を取得するステップと、
前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、
前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、
前記事前共分散行列を用いてカルマンゲインを演算するステップと、
前記状態ベクトル事前推定値と、前記カルマンゲインと、前記センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、
前記事前共分散行列と、前記カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップと
を含む、方法。
【請求項８】
振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、該状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する方法を、コンピュータに実行させるためのプログラムであって、
前記外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、
前記振動センサのセンサ値を取得するステップと、
前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、
前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、
前記事前共分散行列を用いてカルマンゲインを演算するステップと、
前記状態ベクトル事前推定値と、前記カルマンゲインと、前記センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、
前記事前共分散行列と、前記カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップと
をコンピュータに実行させるための、プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、地震計等の振動センサに関する状態を推定する推定装置、これを用いた振動センサシステム、及び関連する方法、プログラムに関する。特に本発明は、複数のセンサ（測定器と呼ぶこともある。）からのデータを統合するセンサフュージョンにより、振動センサに関する状態を推定する推定装置、振動センサシステム、及び関連する方法、プログラムに関する。
【背景技術】
【０００２】
センサフュージョンとは、複数のセンサからのデータを結合して、それらの利点を取り込み、欠点を補い、対象の状態推定や制御を図る技術である。ジャイロスコープと加速度計からなる慣性航法システムにグローバル航行衛星システムからの位置信号を統合した統合慣性航法システム（非特許文献１の第９章）はその代表的な成功例とされる。最近では、ＡＩ（人工知能）技術と連結して、自動車の自動運転にミリ波レーダーをはじめとする様々なセンサを統合するセンサフュージョン技術が話題となっている。
【０００３】
センサフュージョンにおける技法には、Ｂａｙｅｓ推定や最尤法のような統計的推論が用いられるが（非特許文献２）、多用される方法の一つはカルマンフィルタ（非特許文献３）である。実際、非特許文献１などは、カルマンフィルタの教科書の中で、応用分野の一つとして統合慣性航法システムを扱っている。統計的推論で扱うのは、センサからのデータであり、これらを如何に統合して活用するかをセンサフュージョンの目的としている。外向指向の技術である。
【０００４】
地震観測においても、これまでは、高感度短周期速度計、広帯域速度計および強震計などが目的に沿って用いられてきた。最近では、例えば、広帯域速度計と強震計を併設し、より広いダイナミックレンジで地震をリアルタイムで把握する傾向が世界的な潮流となっている。この一体化は、センサ自体の性能向上を目的とするため内的指向であり、グローバル航行衛星システムや自動車の自動運転のような華やかさとは無縁である。とは言え、複数の地震計を融合して一体化するセンサフュージョンは、測定周波数帯域とダイナミックレンジの拡大を追求する地震観測では重要な技術課題の一つである。地震観測におけるセンサフュージョンの基本は二種類の地震計の統合である。例えば、広帯域速度計と短周期速度計のセンサフュージョンにより、測定周波数帯域の拡大を図るものである。これは、相補フィルタリングにより容易に実現可能となる（非特許文献４）。しかしながら、二種類の地震計を融合することにより、ダイナミックレンジの拡大を図る種類のセンサフュージョンは、実用レベルまで至っていないのが、本発明者の知る限りでは現状と考えられる。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Grewal, M.S. and A. P. Andrews (2008). Kalman filtering, Wiley-IEEE Press.
【非特許文献２】鏡慎吾・石川正俊(2005). センサーフュージョン−センサーネットワークの情報処理構造−，  電子情報通信学会論文誌  A， J88-A, No. 12, 1404-1412
【非特許文献３】Kalman, R.E. (1960). A new approach to linear filtering and prediction problems, ASME. Journal of Basic Engineering, 82, 34-45
【非特許文献４】木下繁夫(2014). 状態帰還型換振器と相補フィルタリングによる特性補償，地震，66，101-113
【非特許文献５】Anderson, B.D.O. and J.B.Moore (1979). Optimal filtering, Prentice-Hall, Inc.
【発明の概要】
【発明が解決しようとする課題】
【０００６】
以上に鑑み、本発明は、地震計等、複数の振動センサのセンサフュージョンを行うための実用レベルの手法を提供することを課題とする。
【課題を解決するための手段】
【０００７】
上記課題を解決するべく、本発明は、振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置であって、外部入力として振動に関する物理量の外部測定器により測定された測定値を取得する、測定値取得部と、振動センサのセンサ値を取得する、センサ値取得部と、前回の推定により得られた状態ベクトル推定値と、前回の演算により得られた事後共分散行列と、前回の推定及び前回の演算に対応する時刻に対応する測定値と、システムノイズ（システム雑音）の共分散行列とを格納する、データ格納部と、前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算する、状態ベクトル事前推定値演算部と、前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算する、事前共分散行列演算部と、事前共分散行列を用いてカルマンゲインを演算する、カルマンゲイン演算部と、状態ベクトル事前推定値と、カルマンゲインと、センサ値とを用いて、今回の推定における状態ベクトル推定値を演算する、状態ベクトル推定値演算部と、事前共分散行列と、カルマンゲインとを用いて、今回の演算における事後共分散行列を演算する、事後共分散行列演算部とを備えた、推定装置を提供する。
【０００８】
さらに本発明は、振動センサと、外部測定器と、上記推定装置とを備えた、振動センサシステムを提供する。
【０００９】
上記外部測定器は演算により加速度変換できるもので、前記振動センサのダイナミックレンジの上限を超えるダイナミックレンジを有し測定値取得部は加速度の測定値を取得するものであってよい。
【００１０】
上記振動センサは振り子を備えてよく、振動センサは変位、速度、加速度のうち１以上の値を測定してよく、センサ値取得部は、変位、速度、加速度のうち１以上の値を取得するものであってよい。
【００１１】
上記カルマンゲイン演算部は、カルマンゲイン調整項を用いた演算によりカルマンゲインを調整するものであってよく、カルマンゲインの調整は、振動センサの変位、速度、加速度のうち１以上の値が所定値よりも大きい場合には、１以上の値が所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくすることでカルマンゲインを小さくすることにより行われるものであってよい。
【００１２】
上記振動センサは、振り子を備え、変位、速度、加速度のうち１以上の値を測定する振動センサであり、振り子を備えた振動センサの動作を演算により模擬し、模擬的な振り子の変位、速度、加速度のうち１以上の値を演算により決定する模擬的センサが後段のカルマンフィルタにより備えられ、状態ベクトル推定値演算部は、振動センサにより測定された変位、速度、加速度のうち１以上の値の大きさに応じて、振動センサにより測定された変位、速度、加速度のうち１以上の値と、模擬的センサにより決定された模擬的な振り子の変位、速度、加速度のうち１以上の値に係数を乗じたものとのうちいずれかをセンサ値として用いるものであってよい。
【００１３】
また本発明は、振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置により実行される方法であって、外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、振動センサのセンサ値を取得するステップと、前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、事前共分散行列を用いてカルマンゲインを演算するステップと、状態ベクトル事前推定値と、カルマンゲインと、センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、事前共分散行列と、カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップとを含む、方法を提供する。
【００１４】
また本発明は、振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する方法を、コンピュータに実行させるためのプログラムであって、外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、振動センサのセンサ値を取得するステップと、前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、事前共分散行列を用いてカルマンゲインを演算するステップと、状態ベクトル事前推定値と、カルマンゲインと、センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、事前共分散行列と、カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップとをコンピュータに実行させるための、プログラムを提供する。
【発明の効果】
【００１５】
本発明によれば、複数の振動センサのセンサフュージョンにより、地震計等、複数の振動センサを統合して用いることが可能となる。一例においては、高感度換振器と低感度換振器のセンサフュージョンにより高感度換振器のダイナミックレンジを拡大することが可能となる。
【図面の簡単な説明】
【００１６】
【図１】センサフュージョンのフローを概念的に説明する図。
【図２】カルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図。
【図３】本発明の一実施形態に係る振動センサシステムの構成を示すブロック図。
【図４】推定装置の装置構成例を示すブロック図。
【図５】振動センサ（高感度速度計）の一例を示す図。
【図６】模擬的振動センサの構成を示すブロック線図。
【図７】高感度速度計（ＦＳＦセンサ）の連続時間系状態空間表現を示すブロック線図。
【図８】ＦＳＦセンサの復元加速度の状態ベクトルとフィードバックゲインの関係を示す図。
【図９】ＦＳＦセンサの加速度入力に対するＶＢＢ等価な周波数応答特性（モデル１ａ２）を示すグラフ。
【図１０】速度計の一例の構成図（株式会社東京測振製「ＶＳＥ−１５Ｄ−６」取扱説明書より引用）。
【図１１】振動センサシステムの設計時の動作フローを示すフローチャート。
【図１２】振動センサシステムの運用時の動作フローを示すフローチャート。
【図１３】図１１，図１２の動作フロー中、カルマンフィルタを用いた処理部分の動作フローを示すフローチャート。
【図１４】振動センサ（高感度変位計）の一例を示す図。
【図１５】高感度変位計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図。
【図１６】高感度変位計の連続時間系状態空間表現を示すブロック線図。
【図１７】振動センサ（高感度加速度計）の一例として、高感度負帰還型加速度計（力平行型）を示す図。
【図１８】高感度加速度計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図。
【図１９】高感度加速度計の連続時間系状態空間表現を示すブロック線図。
【発明を実施するための形態】
【００１７】
以下、本発明を実施するための形態を説明する。ただし本発明が以下に説明する具体的態様に限定されるわけではなく、本発明の請求の範囲内で適宜変更可能であることに留意する。
【００１８】
以下の実施形態においては、高感度換振器（高感度振動センサ）と、高感度換振器のダイナミックレンジの上限を超えるダイナミックレンジを有する低感度加速度換振器（低感度振動センサ。「外部測定器」の一例。）のセンサフュージョンにより高感度換振器のダイナミックレンジを拡大する一つの方法を提案する。最初に、広帯域速度計（ＶＢＢ）と強震計のセンサフュージョンを具体例として展開するが、この方法は、三種類の高感度地震計（変位計、速度計及び加速度計）と強震計の統合にも適用可能である。
【００１９】
センサフュージョンの概念
提案するセンサフュージョンのフローの一例を図１に概念的に示す。この例は、観測された強震計の加速度記録と並行観測された高感度換振器からの記録をフュージョン（融合）してダイナミックレンジの拡大を図るものである。このために、高感度換振器と等価な仮想センサ（模擬的センサ）を、ソフトウェアを用いて構築することができる（後述のとおり、仮想センサではなく実際の高感度換振器を用いてセンサフュージョンを行うこともできるし、仮想センサと実際のセンサを組み合わせて使い分けてもよい。）。この仮想センサは、全状態帰還型（ｆｕｌｌ−ｓｔａｔｅ−ｆｅｅｄｂａｃｋ，略してＦＳＦ）センサ（非特許文献４）と呼ばれるものである。このＦＳＦセンサは、高感度換振器と等価であるから、これを状態空間で表現し、入力に観測された加速度を確定された信号として与え、出力に観測された高感度換振器からの記録に正規性雑音を加えて宛がうことが可能となる。ここにおいて「状態」とは、ＦＳＦセンサ内部での振り子運動の状態である。よって、高感度換振器の出力、すなわちＦＳＦセンサの出力がカルマンフィルタを用いて得られた振り子の状態から推定されることとなる。仮想センサとしてのＦＳＦセンサは実機ではないので、電源電圧の制限から逃れられる。このため推定されるＦＳＦセンサの出力は、実観測された高感度換振器からの記録がクリップしても（所定値を超えて飽和しても）、加速度計が正常動作していれば、高感度換振器の実記録に加える観測雑音（観測ノイズ）を大きくしてカルマンゲインを小さくすることで、クリップしない状態で出力される。また、観測雑音を小さくしてカルマンゲインを大きくすれば実観測された波形がカルマンフィルタから忠実に出力される。以降において、「高感度換振器」というときには、実機としての高感度換振器を指す。
【００２０】
仮想センサ（ソフトウェアセンサ）を用いたセンサフュージョンの構成
図２は、カルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図である。説明を簡潔にするべく、ここではＦＳＦセンサが模擬的な速度計であるとするが、後述のとおりＦＳＦセンサは模擬的な変位計、或いは加速度計であってもよいし、またＦＳＦセンサが模擬的な変位計と速度計の両方として動作することにより、低感度換振器（加速度計）と変位計及び速度計とを融合させるなど、多数の振動センサによりセンサフュージョンを行うこともできる。
【００２１】
図２のブロック線図中、
【数１】
0006980199.tif 000002
は低感度加速度計の加速度記録（単位はｍ（メートル）／ｓ2（秒の二乗）とする。以下においても同様。）とし、ｙm（ｎ）は正規性雑音を加えた高感度地震計の記録とする。具体的に、図２の例ではｍ＝２とし、高感度速度計からの速度記録であるとするが（単位はｍ（メートル）／ｓ（秒）とする。以下においても同様。）、ｍ＝３のときは高感度加速度計からの加速度記録であるし（単位はｍ（メートル）／ｓ2（秒の二乗）とする。以下においても同様。）、ｍ＝１のときのｙ1（ｎ）は高感度変位計からの変位記録である（単位はｍ（メートル）とする。以下においても同様。）。ここでｎは整数であり離散時刻を表す。標本化時間をΔＴとすれば、離散時刻ｎに対応する実際の時刻はｎΔＴ（単位はｓとする。以下においても同様。）となり、この時刻ｎΔＴで本来は表記すべきであるが、簡略化のためにΔＴを省略する。また図２中、Σは、そこに向かう信号値を足し合わせて、そこから出る方向に出力することを示す（以下においても同様。なおマイナスの符号が付いている場合には、その信号値は足されるのではなく引かれる。）。ｚは、離散時刻をｎからｎ＋１まで推移させる時間推移演算子であり、ｚ-1は、その逆の演算（離散時刻をｎ＋１からｎまで推移させる）をする演算子である。
【数２】
0006980199.tif 000003
は状態ベクトルを観測データに変換する演算子であり（観測対象に応じた行列又はベクトルで表現される）、離散時刻ｎにおける状態（状態ベクトル）を
【数３】
0006980199.tif 000004
とし（列ベクトル）、離散時刻ｎにおける状態ベクトルのカルマンフィルタによる推定値を
【数４】
0006980199.tif 000005
とすると（列ベクトル）、以下の観測方程式
【数５】
0006980199.tif 000006
により、観測量の推定値が与えられる。
【数６】
0006980199.tif 000007
は３次の単位行列である。
【数７】
0006980199.tif 000008
は状態ベクトルの推移を規定する遷移行列であり、
【数８】
0006980199.tif 000009
は入力信号を選択する列ベクトルであり、それぞれ、モデルに応じて具体的にその値が定められる。
【数９】
0006980199.tif 000010
はカルマンゲインであり、この例においてはベクトル（列ベクトル）として後述のとおりフィルタリングステップ内で決定される。
【００２２】
図２に示すとおり、カルマンフィルタには、低感度加速度計の加速度記録
【数１０】
0006980199.tif 000011
と、高感度速度計の記録
【数１１】
0006980199.tif 000012
に観測雑音ｒ（ｎ）を加えた記録ｙm（ｎ）が入力される。観測雑音ｒ（ｎ）は正規性ノイズであり、その平均値Ｅ［ｒ（ｎ）］が０であるとし、その分散Ｅ［ｒ2（ｎ）］をＲ（ｎ）とする（Ｅは統計的平均を表す。以下においても同様。）。Ｒ（ｎ）は、外部から（振動センサシステムの使用者等により）設定される量であり、Ｒ（ｎ）を調整することによりカルマンゲインの大きさを調整できる。
【００２３】
図２に示すセンサフュージョンにおいては、
【数１２】
0006980199.tif 000013
を外部からの制御入力とみなし、ｙm（ｎ）を観測値とし、これらを用いて線形カルマンフィルタの予測ステップとフィルタリングステップとを実行することにより、ＦＳＦセンサに関する状態が推定される。高感度換振器により実際に観測された観測値に正規性ノイズを加えたｙm（ｎ）を、図２中のｙm（ｎ）として（図２ではｍ＝２としているが、上述のとおりｍは任意。）入力してカルマンフィルタによる推定を行う。図２中のカルマンフィルタによる状態推定は、
【数１３】
0006980199.tif 000014
【数１４】
0006980199.tif 000015
により記述される線形システムに対して行われる。ここで、
【数１５】
0006980199.tif 000016
【数１６】
0006980199.tif 000017
である。
なお、
【数１７】
0006980199.tif 000018
は正規性のシステム雑音（平均がゼロベクトルの列ベクトル）であり、
【数１８】
0006980199.tif 000019
は、
【数１９】
0006980199.tif 000020
の転置行列として与えられる行ベクトルである。本明細書において、［＊］’は、行列（又はベクトル）［＊］の転置行列（ｔｒａｎｓｐｏｓｅｄ  ｍａｔｒｉｘ）又は転置行列として与えられるベクトルであるとする。「’」記号により転置を表す。また、（式４）で表される
【数２０】
0006980199.tif 000021
は、システム雑音の共分散行列である。
【００２４】
振動センサシステムの具体的構成
次に、本発明によるセンサフュージョンを行うための振動センサシステムの具体的構成を説明する。図３は、本発明の一実施形態に係る振動センサシステムの構成を示すブロック図である。振動センサシステム１は、推定装置２と、外部測定器３と、振動センサ４とを備え、推定装置２は、疑似的センサのパラメータを持つカルマンフィルタ演算部２００を備える。
【００２５】
推定装置２は、上述のカルマンフィルタによる推定を行うための装置であり、測定値取得部２０１、データ格納部２０２、センサ値取得部２０３を備える。またカルマンフィルタ演算部２００は測定帯域内で振動センサと等価な疑似的センサを用いてカルマンフィルタを動作させる装置であり、状態ベクトル事前推定値演算部２０４、事前共分散行列演算部２０５、カルマンゲイン演算部２０６、状態ベクトル推定値演算部２０７、事後共分散行列演算部２０８を備える。これらの各種機能部は、状態推定プログラムを実行するためのコンピュータにより実現されてもよい。そのような例としての、推定装置の装置構成例を図４のブロック図に示す。ただし、推定装置２をコンピュータにより実現することは必須ではなく、例えばＡＳＩＣ（ａｐｐｌｉｃａｔｉｏｎ  ｓｐｅｃｉｆｉｃ  ｉｎｔｅｇｒａｔｅｄ  ｃｉｒｃｕｉｔ、特定用途向け集積回路）、ＦＰＧＡ（ｆｉｅｌｄ−ｐｒｏｇｒａｍｍａｂｌｅ  ｇａｔｅ  ａｒｒａｙ）等の集積回路を用いて推定装置２を構成してもよいし、各種演算回路や記憶回路等を組み合わせて推定装置２の各種機能部を構成してもよい。
【００２６】
図４に示す推定装置２は、中央処理装置（ＣＰＵ：  Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）、一時記憶メモリ（ＲＡＭ：Ｒａｎｄｏｍ  ａｃｃｅｓｓ  ｍｅｍｏｒｙ）を備えた処理部２０９と、ハードディスクドライブ、ＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒｉｖｅ）等の記憶デバイスを用いて構成される記憶部２１０と、ＵＳＢ（Ｕｎｉｖｅｒｓａｌ  Ｓｅｒｉａｌ  Ｂｕｓ）規格等に従ってデータを入出力するためのデータ入出力機器、ディスプレイ装置、キーボード、マウス等を備えた入出力部２１１とを備える。記憶部２１０には、状態推定プログラム、ＯＳ（Ｏｐｅｒａｔｉｎｇ  ｓｙｓｔｅｍ：オペレーティングシステム）等の各種プログラム、及び、状態推定プログラム用データ、その他の各種データが記憶されており、ＲＡＭに適宜プログラム、データを読み込みつつ、ＣＰＵが状態推定プログラム等のプログラムを実行することにより、測定値取得部２０１、データ格納部２０２、センサ値取得部２０３、状態ベクトル事前推定値演算部２０４、事前共分散行列演算部２０５、カルマンゲイン演算部２０６、状態ベクトル推定値演算部２０７、事後共分散行列演算部２０８が実現される。
【００２７】
具体的に、まず測定値取得部２０１は、外部測定器３（低感度換振器としての加速度計）から加速度値を取得する機能部であり、状態推定プログラム、各種プログラムを実行するＣＰＵによって制御される入出力部２１１のデータ入出力機器により実現される。データ格納部２０２は、状態推定プログラム、各種プログラムを実行するＣＰＵによって制御される記憶部２１０によって実現され、後述のとおり、繰り返し行われる推定に関し、前回の推定により得られた状態ベクトル推定値と、前回の演算により得られた事後共分散行列と、前回の推定及び前回の演算に対応する時刻に対応する測定値と、システムノイズの共分散行列とを格納する。センサ値取得部２０３は、振動センサ４のセンサ値を取得する機能部であり、状態推定プログラム、各種プログラムを実行するＣＰＵによって制御される入出力部２１１のデータ入出力機器により実現される。状態ベクトル事前推定値演算部２０４は、状態推定プログラム（特に状態ベクトル事前推定値演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、状態ベクトル事前推定値の演算を行う。事前共分散行列演算部２０５は、状態推定プログラム（特に事前共分散行列演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、事前共分散行列の演算を行う。カルマンゲイン演算部２０６は、状態推定プログラム（特にカルマンゲイン演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、カルマンゲインの演算を行う。状態ベクトル推定値演算部２０７は、状態推定プログラム（特に状態ベクトル推定値演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、状態ベクトル推定値の演算を行う。事後共分散行列演算部２０８は、状態推定プログラム（特に事後共分散行列演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、事後共分散行列の演算を行う。
【００２８】
外部測定器３は、ここにおいては低感度換振器としての加速度計（演算により加速度変換できるもので、高感度換振器のダイナミックレンジの上限を超えるダイナミックレンジを有する。）であり、例えばＭＥＭＳ（Ｍｉｃｒｏ  Ｅｌｅｃｔｒｏ  Ｍｅｃｈａｎｉｃａｌ  Ｓｙｓｔｅｍｓ）型の加速度センサを用いて構成される。外部測定器３は、命令信号の外部からの入力、加速度の測定記録の外部への出力等をするための入出力機器も備えている。
【００２９】
（第１の実施形態）
高感度速度計のモデル
次に、本発明の第１の実施形態として、高感度速度計のモデルを用いる振動センサシステムの実施形態を具体的に説明する。図５は、振動センサ４（高感度速度計としてのＶＢＢ型地震計）の一例を示す図である。
【００３０】
ＦＳＦセンサは古典的なＰＩＤ制御（Ｐｒｏｐｏｒｔｉｏｎａｌ−Ｉｎｔｅｇｒａｌ−Ｄｉｆｆｅｒｅｎｔｉａｌ  Ｃｏｎｔｒｏｌｌｅｒ）理論を用いて設計される負帰還型地震計と同様に組み立てられ、図５に示すＶＢＢ型地震計の動作を模擬する。図５の振動センサは、サーボコイルに接続された振り子４００１（振り子とサーボコイルの具体的な接続は、例えば図１０を参照。なお、検定コイル、駆動用コイルは図１０紙面内では断面として示されており、実際は輪状の形状を有している。マグネットも図１０紙面内では断面として示されており、実際は内部が部分的にくりぬかれた円柱のような形状を有している。）、減衰要素４００２（具体例としては空気抵抗によるダンパーであり実体はない）、剛性要素４００３（具体例としては振り子を支えるバネ）、３枚の電気容量板（極板）を備える変位検出器４００４を備え、これらは容器４００５に収納されている。図５に示すとおり、振り子４００１は、変位検出器４００４に含まれる３枚の容量板のうち、中央の容量板に接続されている。また振動センサ４は、振り子の変位−電圧変換増幅器ＡD（便宜上、増幅器ＡDの変換係数をＡD［Ｖ／ｍ］で表す。），微分器ＧS（便宜上、微分器ＧSの係数をＧSで表す。微分器ＧSの作用は、数式では
【数２１】
0006980199.tif 000022
で表される。），インピーダンス要素ＲI（抵抗ＲI［Ω］で表す。），インピーダンス要素ＲD（抵抗ＲD［Ω］で表す。），インピーダンス要素ＲV（抵抗ＲV［Ω］で表す。）を備え、また積分記号
【数２２】
0006980199.tif 000023
で表される積分器（Ｔ（ｓ）は時定数とする。積分器の作用は、数式では
【数２３】
0006980199.tif 000024
で表される。）を備える。なお、図５中では明示していないが、振り子４００１が接続されるサーボコイルのまわりには磁石が配置されており（図１０を参照。）サーボコイルは磁界の中で運動する。
【００３１】
振動センサ４が設置されている地面等が振動した場合に対応して、振動センサ４に入力変位ｗ（ｔ）が与えられると（単位はｍとする。本明細書において、単位について明記しない量はＭＫＳＡ単位系で表される量であるとする。）、振り子４００１が振動する（格納容器４００５内の振り子の、容器に対する相対変位をｘ（ｔ）（単位はｍとする。）とする）。振り子４００１の変位により、変位検出器４００４に含まれる３枚の容量板のそれぞれの間隔が変化する。具体的に、静止状態において、上の容量板と中央の容量板の間隔と、中央の容量板と下の容量板の間隔とが、それぞれｄ［ｍ］であったとすると、図６中に示す方向で振り子４００１がｘだけ変位することにより、振り子４００１に接続された中央の容量板も同様に変位する。これにより、上の容量板と中央の容量板の間隔はｄ−ｘとなり、中央の容量板と下の容量板の間隔はｄ＋ｘとなる。したがって、上の容量板と中央の容量板の間の電気容量と、中央の容量板と下の容量板との間の電気容量がそれぞれ変化し、振り子４００１の変位ｘが電気回路ＡDによって変位に比例した電圧として検出される。この電圧に対して上述の回路要素により微積分等を行い、サーボコイルに電流ｉ（ｔ）を流し、振り子４００１が静止するように制御する。すなわち、磁石（図１０参照。）による磁界中のサーボコイルに電流ｉ（ｔ）が流れる（フィードバック）ことにより力が発生し、この力により振り子が静止するように制御される。具体的には、入力変位ｗ（ｔ）に対して、
【数２４】
0006980199.tif 000025
を復元加速度として発生させ（Ｇ［Ｎ／Ａ］はサーボコイルのモータジェネレータ定数（ｍｏｔｏｒ  ｇｅｎｅｒａｔｏｒ  ｃｏｎｓｔａｎｔ）であり、ｍ［ｋｇ］は振り子４００１の質量である。）、復元加速度
【数２５】
0006980199.tif 000026
を実現する（ｗの上に「・・」を付しているのは、時間による２回の微分を意味する。上に「・」を付す場合には時間による１回の微分を意味する。時間変化するｗ以外の量についても、適宜同様に時間微分を表記する。）。この制御に要した電流から、加速度、速度、変位がセンサ値として得られる。具体的には、微分器ＧSの出力端から振り子４００１の加速度出力が得られ、変位−電圧変換増幅器ＡDの出力端から振り子４００１の速度出力が得られ、積分器の出力端から振り子４００１の変位出力が得られる。
【００３２】
次に、仮想センサであるＦＳＦセンサによって図５に示すＶＢＢ型地震計の動作を模擬するためのモデルを、非特許文献４にならって説明する。まず、仕様として、機構部の振り子系のパラメータとＦＳＦセンサの設計目標を以下のように与える。
［機構部の振り子系］
振り子の質量＝ｍ［ｋｇ］，減衰定数＝ｈ［無次元量］，固有円振動数＝ω0［ｒａｄ］、サーボコイルのモータジェネレータ定数＝Ｇ［Ｎ／Ａ］
［ＦＳＦセンサの特性］
速度出力感度ＳV［Ｖ／（ｍ／ｓ）］，低域遮断円振動数＝ωC［ｒａｄ］，及び、
【数２６】
0006980199.tif 000027
【００３３】
上記仕様に基づき、ＦＳＦセンサは４段階のステップで実現される。まずステップ１では、フィードバックゲインベクトル（列ベクトル）
【数２７】
0006980199.tif 000028
を下記のとおり決定する。ここでは、発振防止用の高域遮断円振動数ω3を制御パラメータとして試みに与えて計算する。通常、安全側に見て１０3［Ｈｚ］程度とする。
［ステップ１：フィードバックゲインベクトル］
【数２８】
0006980199.tif 000029
【数２９】
0006980199.tif 000030
【００３４】
次に、ステップ２では、加速度及び変位信号の感度、ＳA及びＳDを決定する。このとき、積分器の時定数Ｔと微分器の係数ＧSが制御パラメータとして与えられる。逆に、ＳA及びＳDを与えて、ＡD，Ｔ，ＧSを求めても問題はない。
［ステップ２：感度（Ｓｅｎｓｉｔｉｖｉｔｙ）］
【数３０】
0006980199.tif 000031
【００３５】
最後に、ステップ３で、図５のフィードバックインピーダンスＲV，ＲD，ＲIを決定する。
［ステップ３：フィードバックインピーダンス］
【数３１】
0006980199.tif 000032
【００３６】
ステップ４は確認作業であり、得られたＦＳＦセンサの伝達関数の計算である。ここで、ＨA（ｓ），ＨV（ｓ），ＨD（ｓ）が加速度、速度、及び、変位入力に対する、ＦＳＦセンサの加速度、速度、及び、変位出力に関する伝達関数となる。これらの伝達関数の極構造は、共通の特性方程式Δ（ｓ）＝０から求められる。
［ステップ４：伝達関数（Ｔｒａｎｓｆｅｒ  ｆｕｎｃｔｉｏｎ）］
【数３２】
0006980199.tif 000033
【００３７】
上記ステップで構築されたＦＳＦセンサのブロック線図（ｂｌｏｃｋ  ｄｉａｇｒａｍ）は、状態ベクトル：
【数３３】
0006980199.tif 000034
を用いれば、図７に示すとおりとなる。図７において、ＦＳＦセンサの出力端子Ｄ，Ｖ，Ａでは、
【数３４】
0006980199.tif 000035
が出力として得られる（この時点では観測ノイズを導入していないことに留意する。）。また、復元加速度
【数３５】
0006980199.tif 000036
は、状態フィードバックゲイン
【数３６】
0006980199.tif 000037
を用いれば、
【数３７】
0006980199.tif 000038
（内積演算）となる（これが全状態帰還の意味である。）。その詳細は、図８に示すとおりである。図８は図７の復元加速度ｂ（ｔ）と状態ベクトルｘ（ｔ）の部分のみを抜き出した図であり、式１８を視覚化したものである。
【００３８】
図７及び図８のブロック線図を数式表現、すなわち状態空間表現すれば、以下の式のとおりとなる。
【数３８】
0006980199.tif 000039
【数３９】
0006980199.tif 000040
【００３９】
観測ベクトルにおいて速度出力のみの観測をするとすれば、（式２０）は次式となる。
【数４０】
0006980199.tif 000041
【００４０】
ＦＳＦセンサの設計例として、ＶＢＢ型を対象とすると以下のようになる。まず振り子の質量をｍ＝０．０９［ｋｇ］，固有振動数をｆ0＝１．５［Ｈｚ］，減衰定数を０．４，モータジェネレータ定数Ｇ＝５６［Ｎ／Ａ］とする換振器要素を考える。次に、ＶＢＢの特性として、ＳV＝７５０［Ｖ／（ｍ／ｓ）］及び、
【数４１】
0006980199.tif 000042
とする。この条件で、ｆ3＝１．０６２［ｋＨｚ］を与えると、フィードバックゲインが
【数４２】
0006980199.tif 000043
となり、
【数４３】
0006980199.tif 000044
が得られる。更に、
【数４４】
0006980199.tif 000045
とすれば、加速度出力の感度
【数４５】
0006980199.tif 000046
と、変位出力の感度
【数４６】
0006980199.tif 000047
が得られる。設計されたＦＳＦセンサの加速度入力に対する加速度、速度及び変位出力の周波数応答特性は、図９に示すとおりとなる。これは、標準的なＶＢＢの特性である。ＶＢＢの実機と同様に変位（Ｄ）端からの出力は加速度入力に対して直流感度を有する特性となる。図中のゲイン特性のｄＢ表示では、
【数４７】
0006980199.tif 000048
を０［ｄＢ］としている。
【００４１】
以上説明した高感度速度計のモデルをまとめれば、以下のとおりである。
状態ベクトルを
【数４８】
0006980199.tif 000049
としたとき、
【数４９】
0006980199.tif 000050
【数５０】
0006980199.tif 000051
ここで、
【数５１】
0006980199.tif 000052
は観測ノイズを含まない量である。
【００４２】
（式２８）〜（式３０）は連続時間系の表現であるが、標本化時間ΔＴをパラメータとして離散系に変換すると、
【数５２】
0006980199.tif 000053
で
【数５３】
0006980199.tif 000054
の条件下で次式となる。
【数５４】
0006980199.tif 000055
ここで
【数５５】
0006980199.tif 000056
及び
【数５６】
0006980199.tif 000057
【数５７】
0006980199.tif 000058
【００４３】
（式３３）〜（式３６）の線形モデルに対して（システムノイズも観測ノイズもこの段階では導入していない。システムのブロック線図は図６を参照。）、（式２）を用いて既に説明したシステムノイズ
【数５８】
0006980199.tif 000059
を導入する（平均０、共分散行列Ｑ（ｎ）の正規性雑音とする。）。さらに、
【数５９】
0006980199.tif 000060
に観測ノイズｒ（ｎ）を加えた量を、以下のとおり新たにｙ2（ｎ）と定義する。
【数６０】
0006980199.tif 000061
観測ノイズｒ（ｎ）は、平均０，分散Ｒ（ｎ）の正規性雑音であるとする。上記（式３７）が、観測ノイズがある場合の高感度速度計の信号モデルである。
【００４４】
以上から、図２のセンサフュージョンの構成は、システムノイズ、観測ノイズを含む以下の線形システム
【数６１】
0006980199.tif 000062
により記述される。
【００４５】
このモデルにおいて、（式３８）を
【数６２】
0006980199.tif 000063
と記述したとき、
【数６３】
0006980199.tif 000064
であるから、
ＶＢＢがクリップしない範囲では、
【数６４】
0006980199.tif 000065
が充分な精度で得られているとすると、
【数６５】
0006980199.tif 000066
となるため、
【数６６】
0006980199.tif 000067
となる。ＶＢＢがクリップする範囲では、
【数６７】
0006980199.tif 000068
ならば
【数６８】
0006980199.tif 000069
となるが、
【数６９】
0006980199.tif 000070
ではなく、
【数７０】
0006980199.tif 000071
を意味する。ＶＢＢで観測された信号は、このような性質の観測雑音ｒ（ｎ）によってモデル化される。（式３８）では、更に解釈を拡大して、ｒ（ｎ）を平均零、分散Ｒ（ｎ）の時変正規性雑音としてモデル化している。
【００４６】
カルマンフィルタによる状態推定
（式３８）のシステムに、制御入力がある場合の線形カルマンフィルタを適用可能である。この場合のカルマンフィルタの予測ステップ、フィルタリングステップは、以下のとおりである。
【数７１】
0006980199.tif 000072
【００４７】
（式４７）中、
【数７２】
0006980199.tif 000073
は、状態ベクトル
【数７３】
0006980199.tif 000074
の事前状態推定値であり、
【数７４】
0006980199.tif 000075
は、状態ベクトル
【数７５】
0006980199.tif 000076
に関する共分散行列（誤差共分散行列）
【数７６】
0006980199.tif 000077
の事前値である、事前共分散行列（事前誤差共分散行列）である。
【数７７】
0006980199.tif 000078
は、離散時刻ｎ−１における状態ベクトル
【数７８】
0006980199.tif 000079
の推定値であり、
【数７９】
0006980199.tif 000080
は、離散時刻ｎ−１における共分散行列
【数８０】
0006980199.tif 000081
の、カルマンフィルタにより演算される値であり、事後共分散行列である。
【数８１】
0006980199.tif 000082
は、離散時刻ｎにおける状態ベクトル
【数８２】
0006980199.tif 000083
の推定値であり、
【数８３】
0006980199.tif 000084
は、離散時刻ｎにおける共分散行列
【数８４】
0006980199.tif 000085
のカルマンフィルタにより演算される値であり、事後共分散行列である。また、
【数８５】
0006980199.tif 000086
は単位行列（現在の実施形態においては３行３列の行列）である。
【００４８】
以上のとおり、システムノイズと観測ノイズとを導入することにより、カルマンフィルタを用いたセンサフュージョン（ＦＳＦセンサ、或いは高感度換振器の振り子の状態推定）が可能となる。
【００４９】
振動センサシステムの動作フロー
次に振動センサシステム１の全体的な動作フローを設計時と運用時とに分けて説明する。
【００５０】
設計時
図１１は、振動センサシステムの設計時の動作フローを示すフローチャートである。ステップＳ１１０１において、推定装置２の処理部２０９は、記憶部２１０に記憶された疑似的センサ作成プログラムを実行することにより、上記（式６）〜（式３６）を用いて説明した、高感度換振器と等価特性を有する（全）状態帰還型換振器のパラメータ計算を行い模擬的センサを作成する。これにより、（式３３）〜（式３６）で記述される状態帰還型換振器の離散化モデルが生成される（ステップＳ１１０２）。このモデルにおいて用いられる各種パラメータ等のデータは、推定装置２によるカルマンフィルタを用いた処理において用いられるので、各種パラメータ等のデータは推定装置２の記憶部２１０に記憶される（図４参照）。
【００５１】
引き続き、推定装置２の入出力部２１１は、外部測定器３から低感度換振器離散信号（加速度記録）を、そして振動センサ４から高感度換振器離散信号（速度記録。）を、それぞれ時々刻々と取得し（ステップＳ１１０３）、これら取得した記録のデータを記憶部２１０に記憶させる。模擬的センサ（処理部２０９）は、（式３７）に記述されるとおり、高感度換振器離散信号に観測雑音モデルからの観測信号を加え（ステップＳ１１０４）、またシステムノイズも含む（式３８）の線形システムを生成する。以降、模擬的センサは、外部測定器３からの加速度記録を時々刻々と取得し、システムノイズ、観測ノイズをそれぞれの雑音モデル（予め定義されて記憶部２１０に記憶されているものとする。後述のとおり、雑音モデルは事後的に調整可能。）から時々刻々と生成し、離散的な各々の時刻ｎについて、（式３８）のモデルに従い高感度速度記録ｙ2（ｎ）を生成し続ける。なお、（式３８）中の状態ベクトル
【数８６】
0006980199.tif 000087
における、ｎ＝０での初期値は、予め与えられて記憶部２１０内の状態推定プログラム用データとして記憶されているものとする（システムノイズ等、その他に初期値が必要な変数についても同様。）。
【００５２】
次に、推定装置２の処理部２０９は、取得したそれら記録データと、記憶部２１０に記憶された各種データを用いて状態推定プログラムを実行することにより、（式４７）に従い、外生入力を持つカルマンフィルタによる高感度換振器の振り子運動推定を行い、状態ベクトル
【数８７】
0006980199.tif 000088
の推定値
【数８８】
0006980199.tif 000089
の演算と、推定誤差の事後共分散行列
【数８９】
0006980199.tif 000090
の、カルマンフィルタによる値
【数９０】
0006980199.tif 000091
の演算とを、時々刻々と行う（ステップＳ１１０５）。なお、（式４７）中のカルマンゲインにおける、
【数９１】
0006980199.tif 000092
は、観測ノイズの分散であり、外部から設定可能な量（「カルマンゲイン調整項」）である（一例においては、使用者が入出力部２１１を介して入力し、状態推定プログラムを実行する処理部２０９によって、記憶部２１０に状態推定プログラム用データとして記憶される。）。状態推定プログラムを実行する処理部２０９は、振動センサ４から入力される速度記録の示す速度が所定値（一例においては、使用者が入出力部２１１を介して入力し、状態推定プログラムを実行する処理部２０９によって、記憶部２１０に状態推定プログラム用データとして記憶される。）よりも大きい場合には（高感度速度計の速度記録が飽和している場合など）、所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくすることで、カルマンゲインを小さくすることができる。
【００５３】
状態推定プログラムを実行する処理部２０９による、上記カルマンフィルタによる推定により、カルマンフィルタから推定された、高感度換振器の振り子運動から計算される状態帰還型換振器の出力が得られる（Ｓ１１０６）。この出力は、高感度（振動センサ４）及び低感度（外部測定器３）の融合信号であり（Ｓ１１０７）、推定された状態ベクトルの各成分を入出力部２１１のディスプレイ装置によって表示することにより、使用者は推定された振り子の状態を知ることができる。
【００５４】
さらに設計時においては、状態推定プログラムを実行する処理部２０９が、外部測定器３、振動センサ４からの信号により示される加速度、速度等の記録値と、カルマンフィルタを用いて推定された状態ベクトルから演算される加速度、速度等の値との比較判定を行う（ステップＳ１１０８）。一例においては、ある離散時刻において、外部測定器３、振動センサ４からの信号により示される記録値と、カルマンフィルタを用いて推定された状態ベクトルから演算される値との間の差異が所定値を超える場合、状態推定プログラムを実行する処理部２０９は、記録値と演算値との一致性が低いと判断し（ステップＳ１１０８の条件分岐における「ＮＯ」）、観測ノイズの分散の大きさを変更するなどして観測雑音モデルの調整を行う（ステップＳ１１０９）。調整後の観測雑音モデルからの観測雑音を用いて、ステップＳ１１０４〜ステップＳ１１０８までの処理が再度行われる。ステップＳ１１０８において一致性が高いと判断されるまで、これらの処理が繰り返される。ステップＳ１１０８において一致性が高いと判断されると（一例においては、外部測定器３、振動センサ４からの信号により示される記録値と、カルマンフィルタを用いて推定された状態ベクトルから演算される値との間の差異が所定値を超えない場合。ステップＳ１１０８の条件分岐における「ＹＥＳ」）、設計は終了する。
【００５５】
運用時
図１１を用いて説明した設計処理により、（式３８）の線形システム、（式４７）のカルマンフィルタが、各種演算パラメータや雑音モデルを含めて決定される。運用時には、これらを用いて、振動センサ４に関する状態ベクトルの推定が離散的な時刻について時々刻々と行われる（ｎ＝１，２，…とインクリメントしながら推定を繰り返す。）。運用時のフローチャートは図１２に示すとおりである。ステップＳ１２０１，Ｓ１２０２，Ｓ１２０３，Ｓ１２０４，Ｓ１２０５は、それぞれ、設計時のステップＳ１１０３，Ｓ１１０４，Ｓ１１０５，Ｓ１１０６，Ｓ１１０７と同様であってよい。
【００５６】
カルマンフィルタによる演算の詳細
図１１のフローチャート中のステップＳ１１０５，図１２のフローチャート中のステップＳ１２０３において行われるカルマンフィルタによる演算の詳細なフローを、図１３のフローチャートに示す。
【００５７】
まず、状態ベクトル事前推定値演算部２０４は、記憶部２１０から、前回の状態ベクトル推定値データ、前回の事後共分散行列データ、前回の測定値データ、演算パラメータの読み出しを行う（ステップＳ１３０１）。読み出したデータは処理部２０９のＲＡＭに一時記憶される（ＲＡＭへの各種データの一時記憶は、以降の記載、そしてこれまでの記載において適宜省略する。記憶部からの各種データの読み出しも適宜記載を省略する。）。初回の推定時、すなわち（式４７）においてｎ＝１として処理を行う場合、前回の状態ベクトル推定値データ、前回の事後共分散行列データ、前回の測定値データ（ｎ＝０におけるそれぞれのデータ）とは、それぞれ事前に使用者が入出力部２１１を介して記憶部２１０に記憶させる、外部から与えられたそれぞれの初期値データである。
【００５８】
引き続き、状態ベクトル事前推定値演算部２０４は、（式４７）中、予測ステップの式
【数９２】
0006980199.tif 000093
に従い、状態ベクトルの事前推定値
【数９３】
0006980199.tif 000094
を演算し（ステップＳ１３０２）、処理部２０９のＲＡＭに一時記憶する。
【００５９】
次に、事前共分散行列演算部２０５は、（式４７）中、予測ステップの式
【数９４】
0006980199.tif 000095
に従い、事前共分散行列
【数９５】
0006980199.tif 000096
を演算し（ステップＳ１３０３）、処理部２０９のＲＡＭに一時記憶する。
【００６０】
次に、カルマンゲイン演算部２０６は、（式４７）中、フィルタリングステップの式
【数９６】
0006980199.tif 000097
に従い、カルマンゲインを演算し（ステップＳ１３０４）、処理部２０９のＲＡＭに一時記憶する。ここで、カルマンゲイン演算部２０６は、カルマンゲイン調整項
【数９７】
0006980199.tif 000098
を用いた演算によりカルマンゲインを調整する。具体的に、図１１のフロー中のステップＳ１１０３、又は図１２のフロー中のステップＳ１２０１で取得した高感度換振器離散信号の値（ｖ（ｎ））が所定値よりも大きい場合には、その値が所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくすることでカルマンゲインを小さくすることにより、カルマンゲインを調整する。振動センサ４からのセンサ値が、速度記録以外に変位記録、加速度記録である場合も、同様にカルマンゲイン演算部２０６により、各々の所定値と当該センサ値との比較結果に応じてカルマンゲイン調整項を介してカルマンゲインが調整される（振動センサ４からのセンサ値としての振り子の変位、速度、加速度のうち１以上の値が所定値よりも大きい場合には、１以上の値が所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくする）。
【００６１】
次に、状態ベクトル推定値演算部２０７が、（式４７）中、フィルタリングステップの式
【数９８】
0006980199.tif 000099
に従い、状態ベクトル推定値
【数９９】
0006980199.tif 000100
を演算し（ステップＳ１３０５）、処理部２０９のＲＡＭに一時記憶する。
【００６２】
次に、事後共分散行列演算部２０８が、（式４７）中、フィルタリングステップの式
【数１００】
0006980199.tif 000101
に従い、事後共分散行列
【数１０１】
0006980199.tif 000102
を演算し（ステップＳ１３０６）、処理部２０９のＲＡＭに一時記憶する。
【００６３】
次に、ステップＳ１３０７において、状態ベクトル推定値演算部２０７は、ステップＳ１３０５で推定された状態ベクトル推定値を記憶部２１０に記憶し（離散時刻ｎにおける、この状態ベクトル推定値は、離散時刻ｎ＋１のカルマンフィルタ処理において、「前回の状態ベクトル推定値データ」として用いられる。）、また、事後共分散行列演算部２０８は、ステップＳ１３０６で演算した事後共分散行列を記憶部２１０に記憶する（離散時刻ｎにおける、この事後共分散行列は、離散時刻ｎ＋１のカルマンフィルタ処理において、「前回の事後共分散行列データ」として用いられる。）。その他、状態推定プログラムを実行する処理部２０９は、外部測定器３から取得した測定値（離散時刻ｎにおける、この測定値は、離散時刻ｎ＋１のカルマンフィルタ処理において、「前回の測定値データ」として用いられる。）等、種々のデータを記憶部２１０に記憶させる。ステップＳ１３０１〜Ｓ１３０７の処理は、ステップＳ１１０３，ステップ１２０１のデータ信号取得とともに時々刻々と繰り返され（データ取得を先にまとめて行い、その後にカルマンフィルタ処理をまとめて行ってもよいし、ある離散時刻についてデータ取得をした直後にその離散時刻についてカルマンフィルタ処理を行い、離散時刻に１を加えてから（インクリメント）、次の離散時刻についてデータ取得を行い、カルマンフィルタ処理を行うという流れで、データ取得とカルマンフィルタ処理とを同期させてもよい。）、低感度換振器離散信号と高感度換振器離散信号とが取得されなくなるか、或いは離散時刻が所定の終了値に達する等、所定の終了条件が満たされると、カルマンフィルタによる処理は終了する。
【００６４】
図１３で示されるフローのカルマンフィルタ処理を含め、図１２で示される運用時のフローに従う処理が時々刻々と繰り返されることにより、時々刻々と変化する、振動センサ４に関する状態ベクトル（振り子の状態ベクトル）が推定される。状態ベクトル推定値は、状態推定プログラムを実行する処理部２０９により、入出力部２１１のディスプレイ装置に表示される等し、これにより使用者は状態ベクトル推定値を確認することができる。
【００６５】
（第２の実施形態）
高感度変位計のモデル
次に、本発明の第２の実施形態として、振動センサ４として高感度変位計を用いる場合の実施形態を説明する。ただし、第１の実施形態と同様の部分については説明を省略する。
【００６６】
図１４は、振動センサ（高感度変位計としてのＶＢＢ型地震計）の一例を示す図である。図１５は、高感度変位計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図である。高感度変位計は、図５の速度計と同じ構成で実現される。
【００６７】
図１４の振動センサにおいて、変位信号は図１４の＃１の位置から取り出せるが、＃２の位置のように測定帯域外で減衰特性を有する外部積分器を通して得るのが一般的である。測定帯域内での連続時間系状態空間表現は、図１６に示すとおりとなるが、変位信号
【数１０２】
0006980199.tif 000103
とともに、速度信号
【数１０３】
0006980199.tif 000104
も同時に出力されるため、これらを併せて利用することが可能である。
【００６８】
このとき、状態空間表示は、振動センサ４に関する振り子の状態ベクトルを
【数１０４】
0006980199.tif 000105
とすると、
【数１０５】
0006980199.tif 000106
【数１０６】
0006980199.tif 000107
となる。この時間連続系の表現は、標本化時間ΔTをパラメータとして離散系に変換すると、
【数１０７】
0006980199.tif 000108
で
【数１０８】
0006980199.tif 000109
の条件下で次式となる：
【数１０９】
0006980199.tif 000110
【数１１０】
0006980199.tif 000111
【００６９】
したがって、平均ベクトルがゼロベクトル（２次元の列ベクトル）、共分散行列が
【数１１１】
0006980199.tif 000112
の正規性雑音
【数１１２】
0006980199.tif 000113
を考えると、高感度変位計の信号モデルは、
【数１１３】
0006980199.tif 000114
となり、以下のカルマンフィルタ
【数１１４】
0006980199.tif 000115
によるセンサフュージョンが可能となる。振動センサシステム１の装置構成、設計時、運用時の動作フローは、システムノイズ
【数１１５】
0006980199.tif 000116
を第１の実施形態と同様に導入した上で、線形システムの式として
【数１１６】
0006980199.tif 000117
を用い、カルマンフィルタの式として（式６２）の式を用いる以外は、第１の実施形態と同様に実現できる。
【００７０】
（第３の実施形態）
高感度加速度計のモデル
次に、本発明の第３の実施形態として、振動センサ４として高感度加速度計を用いる場合の実施形態を説明する。ただし、第１の実施形態と同様の部分については説明を省略する。
【００７１】
図１７は、振動センサ（高感度加速度計）の一例を示す図である。図１８は、高感度加速度計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図である。高感度加速度計は、図５の速度計と同じ構成でも実現可能である。
【００７２】
図５のＶＢＢ型加速度計の場合、加速度信号は図７のＡ端子から
【数１１７】
0006980199.tif 000118
として出力される。高感度速度計の
【数１１８】
0006980199.tif 000119
を、
【数１１９】
0006980199.tif 000120
で置き換えれば、センサフュージョンが行える。しかしながら、図５の構成の地震計は周波数０での感度も零となり、常用される力平衡型の負帰還型加速度計と比較して不利となる。そこで、ここでは、高感度加速度計として、図１７に示す通常の負帰還型加速度計を扱うことにする。図１７の加速度計に対する連続時間系の状態空間表示は図１９となる。
【００７３】
図１９の状態空間表示は、振動センサ４に関する振り子の状態ベクトルを
【数１２０】
0006980199.tif 000121
としたとき、
【数１２１】
0006980199.tif 000122
【数１２２】
0006980199.tif 000123
となる。この時間連続系の表現は、標本化時間ΔＴをパラメータとして離散系に変換すると、
【数１２３】
0006980199.tif 000124
で
【数１２４】
0006980199.tif 000125
の条件下で次式となる：
【数１２５】
0006980199.tif 000126
【数１２６】
0006980199.tif 000127
【００７４】
したがって、平均０、分散が
【数１２７】
0006980199.tif 000128
の正規性雑音
【数１２８】
0006980199.tif 000129
を考えると、高感度加速度計の信号モデルは、
【数１２９】
0006980199.tif 000130
となり、図１８に示すカルマンフィルタ
【数１３０】
0006980199.tif 000131
によるセンサフュージョンが可能となる。振動センサシステム１の装置構成、設計時、運用時の動作フローは、システムノイズ
【数１３１】
0006980199.tif 000132
を第１の実施形態と同様に導入した上で、線形システムの式として
【数１３２】
0006980199.tif 000133
を用い、カルマンフィルタの式として（式７２）の式を用いる以外は、第１の実施形態と同様に実現できる。
【産業上の利用可能性】
【００７５】
本発明は、地震計をはじめとするあらゆる振動センサに関する状態推定に利用することができる。
【符号の説明】
【００７６】
１            振動センサシステム
２            推定装置
３            外部測定器（低感度換振器）
４            振動センサ（高感度換振器）
２００        カルマンフィルタ演算部
２０１        測定値取得部
２０２        データ格納部
２０３        センサ値取得部
２０４        状態ベクトル事前推定値演算部
２０５        事前共分散行列演算部
２０６        カルマンゲイン演算部
２０７        状態ベクトル推定値演算部
２０８        事後共分散行列演算部
２０９        処理部
２１０        記憶部
２１１        入出力部
４００１      サーボコイルに接続された振り子
４００２      減衰要素
４００３      剛性要素
４００４      変位検出器（３枚の容量板）
４００５      容器
【図１】
0006980199.tif 000134
【図２】
0006980199.tif 000135
【図３】
0006980199.tif 000136
【図４】
0006980199.tif 000137
【図５】
0006980199.tif 000138
【図６】
0006980199.tif 000139
【図７】
0006980199.tif 000140
【図８】
0006980199.tif 000141
【図９】
0006980199.tif 000142
【図１０】
0006980199.tif 000143
【図１１】
0006980199.tif 000144
【図１２】
0006980199.tif 000145
【図１３】
0006980199.tif 000146
【図１４】
0006980199.tif 000147
【図１５】
0006980199.tif 000148
【図１６】
0006980199.tif 000149
【図１７】
0006980199.tif 000150
【図１８】
0006980199.tif 000151
【図１９】
0006980199.tif 000152