特許 6980199 推定装置、振動センサシステム、推定装置により実行される方法、及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6980199

(24)【登録日】2021年11月19日

(45)【発行日】2021年12月15日

(54)【発明の名称】推定装置、振動センサシステム、推定装置により実行される方法、及びプログラム

(51)【国際特許分類】

   G01V 1/00 20060101AFI20211202BHJP

   G01H 1/00 20060101ALI20211202BHJP

   G01H 17/00 20060101ALI20211202BHJP

【ＦＩ】

   G01V1/00 D

   G01H1/00 E

   G01H17/00 Z

【請求項の数】8

【全頁数】41

(21)【出願番号】特願2020-28155(P2020-28155)

(22)【出願日】2020年2月21日

(65)【公開番号】特開2021-131361(P2021-131361A)

(43)【公開日】2021年9月9日

【審査請求日】2021年3月11日

(73)【特許権者】

【識別番号】591088504

【氏名又は名称】株式会社東京測振

(73)【特許権者】

【識別番号】597103872

【氏名又は名称】倉橋護謨工業株式会社

(74)【代理人】

【識別番号】100094569

【弁理士】

【氏名又は名称】田中 伸一郎

(74)【代理人】

【識別番号】100103610

【弁理士】

【氏名又は名称】▲吉▼田 和彦

(74)【代理人】

【識別番号】100109070

【弁理士】

【氏名又は名称】須田 洋之

(74)【代理人】

【識別番号】100067013

【弁理士】

【氏名又は名称】大塚 文昭

(74)【代理人】

【識別番号】100086771

【弁理士】

【氏名又は名称】西島 孝喜

(74)【代理人】

【識別番号】100109335

【弁理士】

【氏名又は名称】上杉 浩

(74)【代理人】

【識別番号】100120525

【弁理士】

【氏名又は名称】近藤 直樹

(74)【代理人】

【識別番号】100139712

【弁理士】

【氏名又は名称】那須 威夫

(74)【代理人】

【識別番号】100158469

【弁理士】

【氏名又は名称】大浦 博司

(72)【発明者】

【氏名】木下 繁夫

(72)【発明者】

【氏名】佐藤 慎祐

【審査官】 佐野 浩樹

(56)【参考文献】

【文献】 特開平１０−０４５３７９（ＪＰ，Ａ）

【文献】 特開２００４−３４７５０９（ＪＰ，Ａ）

【文献】 特開平０７−１０１６７３（ＪＰ，Ａ）

【文献】 特開２０１３−０８８１６２（ＪＰ，Ａ）

【文献】 特開２００３−０９０７５７（ＪＰ，Ａ）

【文献】 特開２０１９−１８９１２１（ＪＰ，Ａ）

【文献】 特開２０１７−１２３６３０（ＪＰ，Ａ）

【文献】 特開２０１６−０４５７６７（ＪＰ，Ａ）

【文献】 米国特許出願公開第２０１３／０１４４９２３（ＵＳ，Ａ１）

【文献】 米国特許出願公開第２００６／００５０６１２（ＵＳ，Ａ１）

【文献】 中国特許出願公開第１０３１４９５９１（ＣＮ，Ａ）

【文献】 足立 修一，線形カルマンフィルタの基礎，計測と制御，日本，2017年，第５６巻、第９号，p. 632-637，＜検索日：2021/04/20＞，＜DOI:https://doi.org/10.11499/sicejl.56.632＞, <https://doi.org/10.11499/sicejl.56.632>からダウンロード可能

【文献】 鏡 慎吾，石川 正俊，センサフュージョン―センサネットワークの情報処理構造―，電子情報通信学会論文誌 A，2005年，Vol. J88-A, No. 12，p. 1404-1412，<検索日：2021/04/16>, <https://search.ieice.org/bin/index.php?category=A&lang=J&vol=J88-A&num=12&abst=>からダウンロード可能

【文献】 KALMAN, R. E.，A New Approach to Linear Filtering and Prediction Problems，Transactions of the ASME--Journal of Basic Engineering，1960年，Vol. 82, No. Series D，p. 35-45，<検索日：2021/4/16>, <https://www.cs.unc.edu/~welch/kalman/kalmanPaper.html>からダウンロード可能

【文献】 木下 繁夫，状態帰還型換振器と相補フィルタリングによる特性補償，地震，2014年，第2輯，第66巻，第4号，p. 113-125，<検索日：2021/4/16>，<DOI：https://doi.org/10.4294/zisin.66.113>からダウンロード可能

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ０１Ｈ１／００−１７／００

Ｇ０１Ｖ１／００−９９／００

(57)【特許請求の範囲】

【請求項1】

振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、該状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置であって、
前記外部入力として振動に関する物理量の外部測定器により測定された測定値を取得する、測定値取得部と、
前記振動センサのセンサ値を取得する、センサ値取得部と、
前回の推定により得られた状態ベクトル推定値と、前回の演算により得られた事後共分散行列と、該前回の推定及び該前回の演算に対応する時刻に対応する前記測定値と、システムノイズの共分散行列とを格納する、データ格納部と、
前記前回の推定により得られた状態ベクトル推定値と、前記前回の推定及び前記前回の演算に対応する時刻に対応する前記測定値とを用いて、状態ベクトル事前推定値を演算する、状態ベクトル事前推定値演算部と、
前記前回の演算により得られた事後共分散行列と、前記システムノイズの共分散行列とを用いて、事前共分散行列を演算する、事前共分散行列演算部と、
前記事前共分散行列を用いてカルマンゲインを演算する、カルマンゲイン演算部と、
前記状態ベクトル事前推定値と、前記カルマンゲインと、前記センサ値とを用いて、今回の推定における状態ベクトル推定値を演算する、状態ベクトル推定値演算部と、
前記事前共分散行列と、前記カルマンゲインとを用いて、今回の演算における事後共分散行列を演算する、事後共分散行列演算部と
を備えた、推定装置。

【請求項2】

前記振動センサと、
前記外部測定器と、
請求項１に記載の推定装置と
を備えた、振動センサシステム。

【請求項3】

前記外部測定器は演算により加速度変換できるもので、前記振動センサのダイナミックレンジの上限を超えるダイナミックレンジを有し、前記測定値取得部は加速度の測定値を取得する、請求項２に記載の振動センサシステム。

【請求項4】

前記振動センサは振り子を備え、変位、速度、加速度のうち１以上の値を測定し、前記センサ値取得部は、該変位、該速度、該加速度のうち１以上の値を取得する、請求項２又は３に記載の振動センサシステム。

【請求項5】

前記カルマンゲイン演算部は、カルマンゲイン調整項を用いた演算によりカルマンゲインを調整し、該カルマンゲインの調整は、前記振動センサの変位、速度、加速度のうち１以上の値が所定値よりも大きい場合には、該１以上の値が所定値よりも小さい場合に比べて該カルマンゲイン調整項を大きくすることでカルマンゲインを小さくすることにより行われる、請求項４に記載の振動センサシステム。

【請求項6】

前記振動センサは、
振り子を備え、変位、速度、加速度のうち１以上の値を測定する振動センサであり、
振り子を備えた振動センサの動作を演算により模擬し、模擬的な振り子の変位、速度、加速度のうち１以上の値を演算により決定する模擬的センサが後段のカルマンフィルタにより備えられ、と
前記状態ベクトル推定値演算部は、前記振動センサにより測定された変位、速度、加速度のうち１以上の値の大きさに応じて、該振動センサにより測定された変位、速度、加速度のうち１以上の値と、前記模擬的センサにより決定された前記模擬的な振り子の変位、速度、加速度のうち１以上の値に係数を乗じたものとのうちいずれかを前記センサ値として用いる、請求項２又は３に記載のシステム。

【請求項7】

振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、該状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置により実行される方法であって、
前記外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、
前記振動センサのセンサ値を取得するステップと、
前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、
前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、
前記事前共分散行列を用いてカルマンゲインを演算するステップと、
前記状態ベクトル事前推定値と、前記カルマンゲインと、前記センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、
前記事前共分散行列と、前記カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップと
を含む、方法。

【請求項8】

振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、該状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する方法を、コンピュータに実行させるためのプログラムであって、
前記外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、
前記振動センサのセンサ値を取得するステップと、
前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、
前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、
前記事前共分散行列を用いてカルマンゲインを演算するステップと、
前記状態ベクトル事前推定値と、前記カルマンゲインと、前記センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、
前記事前共分散行列と、前記カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップと
をコンピュータに実行させるための、プログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、地震計等の振動センサに関する状態を推定する推定装置、これを用いた振動センサシステム、及び関連する方法、プログラムに関する。特に本発明は、複数のセンサ（測定器と呼ぶこともある。）からのデータを統合するセンサフュージョンにより、振動センサに関する状態を推定する推定装置、振動センサシステム、及び関連する方法、プログラムに関する。

【背景技術】

【0002】

センサフュージョンとは、複数のセンサからのデータを結合して、それらの利点を取り込み、欠点を補い、対象の状態推定や制御を図る技術である。ジャイロスコープと加速度計からなる慣性航法システムにグローバル航行衛星システムからの位置信号を統合した統合慣性航法システム（非特許文献１の第９章）はその代表的な成功例とされる。最近では、ＡＩ（人工知能）技術と連結して、自動車の自動運転にミリ波レーダーをはじめとする様々なセンサを統合するセンサフュージョン技術が話題となっている。

【0003】

センサフュージョンにおける技法には、Ｂａｙｅｓ推定や最尤法のような統計的推論が用いられるが（非特許文献２）、多用される方法の一つはカルマンフィルタ（非特許文献３）である。実際、非特許文献１などは、カルマンフィルタの教科書の中で、応用分野の一つとして統合慣性航法システムを扱っている。統計的推論で扱うのは、センサからのデータであり、これらを如何に統合して活用するかをセンサフュージョンの目的としている。外向指向の技術である。

【0004】

地震観測においても、これまでは、高感度短周期速度計、広帯域速度計および強震計などが目的に沿って用いられてきた。最近では、例えば、広帯域速度計と強震計を併設し、より広いダイナミックレンジで地震をリアルタイムで把握する傾向が世界的な潮流となっている。この一体化は、センサ自体の性能向上を目的とするため内的指向であり、グローバル航行衛星システムや自動車の自動運転のような華やかさとは無縁である。とは言え、複数の地震計を融合して一体化するセンサフュージョンは、測定周波数帯域とダイナミックレンジの拡大を追求する地震観測では重要な技術課題の一つである。地震観測におけるセンサフュージョンの基本は二種類の地震計の統合である。例えば、広帯域速度計と短周期速度計のセンサフュージョンにより、測定周波数帯域の拡大を図るものである。これは、相補フィルタリングにより容易に実現可能となる（非特許文献４）。しかしながら、二種類の地震計を融合することにより、ダイナミックレンジの拡大を図る種類のセンサフュージョンは、実用レベルまで至っていないのが、本発明者の知る限りでは現状と考えられる。

【先行技術文献】

【非特許文献】

【0005】

【非特許文献1】Grewal, M.S. and A. P. Andrews (2008). Kalman filtering, Wiley-IEEE Press.

【非特許文献2】鏡慎吾・石川正俊(2005). センサーフュージョン−センサーネットワークの情報処理構造−， 電子情報通信学会論文誌 A， J88-A, No. 12, 1404-1412

【非特許文献3】Kalman, R.E. (1960). A new approach to linear filtering and prediction problems, ASME. Journal of Basic Engineering, 82, 34-45

【非特許文献4】木下繁夫(2014). 状態帰還型換振器と相補フィルタリングによる特性補償，地震，66，101-113

【非特許文献5】Anderson, B.D.O. and J.B.Moore (1979). Optimal filtering, Prentice-Hall, Inc.

【発明の概要】

【発明が解決しようとする課題】

【0006】

以上に鑑み、本発明は、地震計等、複数の振動センサのセンサフュージョンを行うための実用レベルの手法を提供することを課題とする。

【課題を解決するための手段】

【0007】

上記課題を解決するべく、本発明は、振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置であって、外部入力として振動に関する物理量の外部測定器により測定された測定値を取得する、測定値取得部と、振動センサのセンサ値を取得する、センサ値取得部と、前回の推定により得られた状態ベクトル推定値と、前回の演算により得られた事後共分散行列と、前回の推定及び前回の演算に対応する時刻に対応する測定値と、システムノイズ（システム雑音）の共分散行列とを格納する、データ格納部と、前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算する、状態ベクトル事前推定値演算部と、前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算する、事前共分散行列演算部と、事前共分散行列を用いてカルマンゲインを演算する、カルマンゲイン演算部と、状態ベクトル事前推定値と、カルマンゲインと、センサ値とを用いて、今回の推定における状態ベクトル推定値を演算する、状態ベクトル推定値演算部と、事前共分散行列と、カルマンゲインとを用いて、今回の演算における事後共分散行列を演算する、事後共分散行列演算部とを備えた、推定装置を提供する。

【0008】

さらに本発明は、振動センサと、外部測定器と、上記推定装置とを備えた、振動センサシステムを提供する。

【0009】

上記外部測定器は演算により加速度変換できるもので、前記振動センサのダイナミックレンジの上限を超えるダイナミックレンジを有し測定値取得部は加速度の測定値を取得するものであってよい。

【0010】

上記振動センサは振り子を備えてよく、振動センサは変位、速度、加速度のうち１以上の値を測定してよく、センサ値取得部は、変位、速度、加速度のうち１以上の値を取得するものであってよい。

【0011】

上記カルマンゲイン演算部は、カルマンゲイン調整項を用いた演算によりカルマンゲインを調整するものであってよく、カルマンゲインの調整は、振動センサの変位、速度、加速度のうち１以上の値が所定値よりも大きい場合には、１以上の値が所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくすることでカルマンゲインを小さくすることにより行われるものであってよい。

【0012】

上記振動センサは、振り子を備え、変位、速度、加速度のうち１以上の値を測定する振動センサであり、振り子を備えた振動センサの動作を演算により模擬し、模擬的な振り子の変位、速度、加速度のうち１以上の値を演算により決定する模擬的センサが後段のカルマンフィルタにより備えられ、状態ベクトル推定値演算部は、振動センサにより測定された変位、速度、加速度のうち１以上の値の大きさに応じて、振動センサにより測定された変位、速度、加速度のうち１以上の値と、模擬的センサにより決定された模擬的な振り子の変位、速度、加速度のうち１以上の値に係数を乗じたものとのうちいずれかをセンサ値として用いるものであってよい。

【0013】

また本発明は、振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する、推定装置により実行される方法であって、外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、振動センサのセンサ値を取得するステップと、前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、事前共分散行列を用いてカルマンゲインを演算するステップと、状態ベクトル事前推定値と、カルマンゲインと、センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、事前共分散行列と、カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップとを含む、方法を提供する。

【0014】

また本発明は、振動センサに関する状態ベクトルであって外部入力に応じて時間変化する状態ベクトルを離散的な時刻に応じて推定するとともに、状態ベクトルに関する共分散行列を離散的な時刻に応じて演算する方法を、コンピュータに実行させるためのプログラムであって、外部入力として振動に関する物理量の外部測定器により測定された測定値を取得するステップと、振動センサのセンサ値を取得するステップと、前回の推定により得られた状態ベクトル推定値と、前回の推定及び前回の演算に対応する時刻に対応する測定値とを用いて、状態ベクトル事前推定値を演算するステップと、前回の演算により得られた事後共分散行列と、システムノイズの共分散行列とを用いて、事前共分散行列を演算するステップと、事前共分散行列を用いてカルマンゲインを演算するステップと、状態ベクトル事前推定値と、カルマンゲインと、センサ値とを用いて、今回の推定における状態ベクトル推定値を演算するステップと、事前共分散行列と、カルマンゲインとを用いて、今回の演算における事後共分散行列を演算するステップとをコンピュータに実行させるための、プログラムを提供する。

【発明の効果】

【0015】

本発明によれば、複数の振動センサのセンサフュージョンにより、地震計等、複数の振動センサを統合して用いることが可能となる。一例においては、高感度換振器と低感度換振器のセンサフュージョンにより高感度換振器のダイナミックレンジを拡大することが可能となる。

【図面の簡単な説明】

【0016】

【図1】センサフュージョンのフローを概念的に説明する図。

【図2】カルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図。

【図3】本発明の一実施形態に係る振動センサシステムの構成を示すブロック図。

【図4】推定装置の装置構成例を示すブロック図。

【図5】振動センサ（高感度速度計）の一例を示す図。

【図6】模擬的振動センサの構成を示すブロック線図。

【図7】高感度速度計（ＦＳＦセンサ）の連続時間系状態空間表現を示すブロック線図。

【図8】ＦＳＦセンサの復元加速度の状態ベクトルとフィードバックゲインの関係を示す図。

【図9】ＦＳＦセンサの加速度入力に対するＶＢＢ等価な周波数応答特性（モデル１ａ２）を示すグラフ。

【図10】速度計の一例の構成図（株式会社東京測振製「ＶＳＥ−１５Ｄ−６」取扱説明書より引用）。

【図11】振動センサシステムの設計時の動作フローを示すフローチャート。

【図12】振動センサシステムの運用時の動作フローを示すフローチャート。

【図13】図１１，図１２の動作フロー中、カルマンフィルタを用いた処理部分の動作フローを示すフローチャート。

【図14】振動センサ（高感度変位計）の一例を示す図。

【図15】高感度変位計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図。

【図16】高感度変位計の連続時間系状態空間表現を示すブロック線図。

【図17】振動センサ（高感度加速度計）の一例として、高感度負帰還型加速度計（力平行型）を示す図。

【図18】高感度加速度計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図。

【図19】高感度加速度計の連続時間系状態空間表現を示すブロック線図。

【発明を実施するための形態】

【0017】

以下、本発明を実施するための形態を説明する。ただし本発明が以下に説明する具体的態様に限定されるわけではなく、本発明の請求の範囲内で適宜変更可能であることに留意する。

【0018】

以下の実施形態においては、高感度換振器（高感度振動センサ）と、高感度換振器のダイナミックレンジの上限を超えるダイナミックレンジを有する低感度加速度換振器（低感度振動センサ。「外部測定器」の一例。）のセンサフュージョンにより高感度換振器のダイナミックレンジを拡大する一つの方法を提案する。最初に、広帯域速度計（ＶＢＢ）と強震計のセンサフュージョンを具体例として展開するが、この方法は、三種類の高感度地震計（変位計、速度計及び加速度計）と強震計の統合にも適用可能である。

【0019】

センサフュージョンの概念
提案するセンサフュージョンのフローの一例を図１に概念的に示す。この例は、観測された強震計の加速度記録と並行観測された高感度換振器からの記録をフュージョン（融合）してダイナミックレンジの拡大を図るものである。このために、高感度換振器と等価な仮想センサ（模擬的センサ）を、ソフトウェアを用いて構築することができる（後述のとおり、仮想センサではなく実際の高感度換振器を用いてセンサフュージョンを行うこともできるし、仮想センサと実際のセンサを組み合わせて使い分けてもよい。）。この仮想センサは、全状態帰還型（ｆｕｌｌ−ｓｔａｔｅ−ｆｅｅｄｂａｃｋ，略してＦＳＦ）センサ（非特許文献４）と呼ばれるものである。このＦＳＦセンサは、高感度換振器と等価であるから、これを状態空間で表現し、入力に観測された加速度を確定された信号として与え、出力に観測された高感度換振器からの記録に正規性雑音を加えて宛がうことが可能となる。ここにおいて「状態」とは、ＦＳＦセンサ内部での振り子運動の状態である。よって、高感度換振器の出力、すなわちＦＳＦセンサの出力がカルマンフィルタを用いて得られた振り子の状態から推定されることとなる。仮想センサとしてのＦＳＦセンサは実機ではないので、電源電圧の制限から逃れられる。このため推定されるＦＳＦセンサの出力は、実観測された高感度換振器からの記録がクリップしても（所定値を超えて飽和しても）、加速度計が正常動作していれば、高感度換振器の実記録に加える観測雑音（観測ノイズ）を大きくしてカルマンゲインを小さくすることで、クリップしない状態で出力される。また、観測雑音を小さくしてカルマンゲインを大きくすれば実観測された波形がカルマンフィルタから忠実に出力される。以降において、「高感度換振器」というときには、実機としての高感度換振器を指す。

【0020】

仮想センサ（ソフトウェアセンサ）を用いたセンサフュージョンの構成
図２は、カルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図である。説明を簡潔にするべく、ここではＦＳＦセンサが模擬的な速度計であるとするが、後述のとおりＦＳＦセンサは模擬的な変位計、或いは加速度計であってもよいし、またＦＳＦセンサが模擬的な変位計と速度計の両方として動作することにより、低感度換振器（加速度計）と変位計及び速度計とを融合させるなど、多数の振動センサによりセンサフュージョンを行うこともできる。

【0021】

図２のブロック線図中、

【数1】

は低感度加速度計の加速度記録（単位はｍ（メートル）／ｓ²（秒の二乗）とする。以下においても同様。）とし、ｙ_m（ｎ）は正規性雑音を加えた高感度地震計の記録とする。具体的に、図２の例ではｍ＝２とし、高感度速度計からの速度記録であるとするが（単位はｍ（メートル）／ｓ（秒）とする。以下においても同様。）、ｍ＝３のときは高感度加速度計からの加速度記録であるし（単位はｍ（メートル）／ｓ²（秒の二乗）とする。以下においても同様。）、ｍ＝１のときのｙ₁（ｎ）は高感度変位計からの変位記録である（単位はｍ（メートル）とする。以下においても同様。）。ここでｎは整数であり離散時刻を表す。標本化時間をΔＴとすれば、離散時刻ｎに対応する実際の時刻はｎΔＴ（単位はｓとする。以下においても同様。）となり、この時刻ｎΔＴで本来は表記すべきであるが、簡略化のためにΔＴを省略する。また図２中、Σは、そこに向かう信号値を足し合わせて、そこから出る方向に出力することを示す（以下においても同様。なおマイナスの符号が付いている場合には、その信号値は足されるのではなく引かれる。）。ｚは、離散時刻をｎからｎ＋１まで推移させる時間推移演算子であり、ｚ^-1は、その逆の演算（離散時刻をｎ＋１からｎまで推移させる）をする演算子である。

【数2】

は状態ベクトルを観測データに変換する演算子であり（観測対象に応じた行列又はベクトルで表現される）、離散時刻ｎにおける状態（状態ベクトル）を

【数3】

とし（列ベクトル）、離散時刻ｎにおける状態ベクトルのカルマンフィルタによる推定値を

【数4】

とすると（列ベクトル）、以下の観測方程式

【数5】

により、観測量の推定値が与えられる。

【数6】

は３次の単位行列である。

【数7】

は状態ベクトルの推移を規定する遷移行列であり、

【数8】

は入力信号を選択する列ベクトルであり、それぞれ、モデルに応じて具体的にその値が定められる。

【数9】

はカルマンゲインであり、この例においてはベクトル（列ベクトル）として後述のとおりフィルタリングステップ内で決定される。

【0022】

図２に示すとおり、カルマンフィルタには、低感度加速度計の加速度記録

【数10】

と、高感度速度計の記録

【数11】

に観測雑音ｒ（ｎ）を加えた記録ｙ_m（ｎ）が入力される。観測雑音ｒ（ｎ）は正規性ノイズであり、その平均値Ｅ［ｒ（ｎ）］が０であるとし、その分散Ｅ［ｒ²（ｎ）］をＲ（ｎ）とする（Ｅは統計的平均を表す。以下においても同様。）。Ｒ（ｎ）は、外部から（振動センサシステムの使用者等により）設定される量であり、Ｒ（ｎ）を調整することによりカルマンゲインの大きさを調整できる。

【0023】

図２に示すセンサフュージョンにおいては、

【数12】

を外部からの制御入力とみなし、ｙ_m（ｎ）を観測値とし、これらを用いて線形カルマンフィルタの予測ステップとフィルタリングステップとを実行することにより、ＦＳＦセンサに関する状態が推定される。高感度換振器により実際に観測された観測値に正規性ノイズを加えたｙ_m（ｎ）を、図２中のｙ_m（ｎ）として（図２ではｍ＝２としているが、上述のとおりｍは任意。）入力してカルマンフィルタによる推定を行う。図２中のカルマンフィルタによる状態推定は、

【数13】

【数14】

により記述される線形システムに対して行われる。ここで、

【数15】

【数16】

である。
なお、

【数17】

は正規性のシステム雑音（平均がゼロベクトルの列ベクトル）であり、

【数18】

は、

【数19】

の転置行列として与えられる行ベクトルである。本明細書において、［＊］’は、行列（又はベクトル）［＊］の転置行列（ｔｒａｎｓｐｏｓｅｄ ｍａｔｒｉｘ）又は転置行列として与えられるベクトルであるとする。「’」記号により転置を表す。また、（式４）で表される

【数20】

は、システム雑音の共分散行列である。

【0024】

振動センサシステムの具体的構成
次に、本発明によるセンサフュージョンを行うための振動センサシステムの具体的構成を説明する。図３は、本発明の一実施形態に係る振動センサシステムの構成を示すブロック図である。振動センサシステム１は、推定装置２と、外部測定器３と、振動センサ４とを備え、推定装置２は、疑似的センサのパラメータを持つカルマンフィルタ演算部２００を備える。

【0025】

推定装置２は、上述のカルマンフィルタによる推定を行うための装置であり、測定値取得部２０１、データ格納部２０２、センサ値取得部２０３を備える。またカルマンフィルタ演算部２００は測定帯域内で振動センサと等価な疑似的センサを用いてカルマンフィルタを動作させる装置であり、状態ベクトル事前推定値演算部２０４、事前共分散行列演算部２０５、カルマンゲイン演算部２０６、状態ベクトル推定値演算部２０７、事後共分散行列演算部２０８を備える。これらの各種機能部は、状態推定プログラムを実行するためのコンピュータにより実現されてもよい。そのような例としての、推定装置の装置構成例を図４のブロック図に示す。ただし、推定装置２をコンピュータにより実現することは必須ではなく、例えばＡＳＩＣ（ａｐｐｌｉｃａｔｉｏｎ ｓｐｅｃｉｆｉｃ ｉｎｔｅｇｒａｔｅｄ ｃｉｒｃｕｉｔ、特定用途向け集積回路）、ＦＰＧＡ（ｆｉｅｌｄ−ｐｒｏｇｒａｍｍａｂｌｅ ｇａｔｅ ａｒｒａｙ）等の集積回路を用いて推定装置２を構成してもよいし、各種演算回路や記憶回路等を組み合わせて推定装置２の各種機能部を構成してもよい。

【0026】

図４に示す推定装置２は、中央処理装置（ＣＰＵ： Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、一時記憶メモリ（ＲＡＭ：Ｒａｎｄｏｍ ａｃｃｅｓｓ ｍｅｍｏｒｙ）を備えた処理部２０９と、ハードディスクドライブ、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）等の記憶デバイスを用いて構成される記憶部２１０と、ＵＳＢ（Ｕｎｉｖｅｒｓａｌ Ｓｅｒｉａｌ Ｂｕｓ）規格等に従ってデータを入出力するためのデータ入出力機器、ディスプレイ装置、キーボード、マウス等を備えた入出力部２１１とを備える。記憶部２１０には、状態推定プログラム、ＯＳ（Ｏｐｅｒａｔｉｎｇ ｓｙｓｔｅｍ：オペレーティングシステム）等の各種プログラム、及び、状態推定プログラム用データ、その他の各種データが記憶されており、ＲＡＭに適宜プログラム、データを読み込みつつ、ＣＰＵが状態推定プログラム等のプログラムを実行することにより、測定値取得部２０１、データ格納部２０２、センサ値取得部２０３、状態ベクトル事前推定値演算部２０４、事前共分散行列演算部２０５、カルマンゲイン演算部２０６、状態ベクトル推定値演算部２０７、事後共分散行列演算部２０８が実現される。

【0027】

具体的に、まず測定値取得部２０１は、外部測定器３（低感度換振器としての加速度計）から加速度値を取得する機能部であり、状態推定プログラム、各種プログラムを実行するＣＰＵによって制御される入出力部２１１のデータ入出力機器により実現される。データ格納部２０２は、状態推定プログラム、各種プログラムを実行するＣＰＵによって制御される記憶部２１０によって実現され、後述のとおり、繰り返し行われる推定に関し、前回の推定により得られた状態ベクトル推定値と、前回の演算により得られた事後共分散行列と、前回の推定及び前回の演算に対応する時刻に対応する測定値と、システムノイズの共分散行列とを格納する。センサ値取得部２０３は、振動センサ４のセンサ値を取得する機能部であり、状態推定プログラム、各種プログラムを実行するＣＰＵによって制御される入出力部２１１のデータ入出力機器により実現される。状態ベクトル事前推定値演算部２０４は、状態推定プログラム（特に状態ベクトル事前推定値演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、状態ベクトル事前推定値の演算を行う。事前共分散行列演算部２０５は、状態推定プログラム（特に事前共分散行列演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、事前共分散行列の演算を行う。カルマンゲイン演算部２０６は、状態推定プログラム（特にカルマンゲイン演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、カルマンゲインの演算を行う。状態ベクトル推定値演算部２０７は、状態推定プログラム（特に状態ベクトル推定値演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、状態ベクトル推定値の演算を行う。事後共分散行列演算部２０８は、状態推定プログラム（特に事後共分散行列演算モジュール）、各種プログラムを実行するＣＰＵによって実現される機能部であり、事後共分散行列の演算を行う。

【0028】

外部測定器３は、ここにおいては低感度換振器としての加速度計（演算により加速度変換できるもので、高感度換振器のダイナミックレンジの上限を超えるダイナミックレンジを有する。）であり、例えばＭＥＭＳ（Ｍｉｃｒｏ Ｅｌｅｃｔｒｏ Ｍｅｃｈａｎｉｃａｌ Ｓｙｓｔｅｍｓ）型の加速度センサを用いて構成される。外部測定器３は、命令信号の外部からの入力、加速度の測定記録の外部への出力等をするための入出力機器も備えている。

【0029】

（第１の実施形態）
高感度速度計のモデル
次に、本発明の第１の実施形態として、高感度速度計のモデルを用いる振動センサシステムの実施形態を具体的に説明する。図５は、振動センサ４（高感度速度計としてのＶＢＢ型地震計）の一例を示す図である。

【0030】

ＦＳＦセンサは古典的なＰＩＤ制御（Ｐｒｏｐｏｒｔｉｏｎａｌ−Ｉｎｔｅｇｒａｌ−Ｄｉｆｆｅｒｅｎｔｉａｌ Ｃｏｎｔｒｏｌｌｅｒ）理論を用いて設計される負帰還型地震計と同様に組み立てられ、図５に示すＶＢＢ型地震計の動作を模擬する。図５の振動センサは、サーボコイルに接続された振り子４００１（振り子とサーボコイルの具体的な接続は、例えば図１０を参照。なお、検定コイル、駆動用コイルは図１０紙面内では断面として示されており、実際は輪状の形状を有している。マグネットも図１０紙面内では断面として示されており、実際は内部が部分的にくりぬかれた円柱のような形状を有している。）、減衰要素４００２（具体例としては空気抵抗によるダンパーであり実体はない）、剛性要素４００３（具体例としては振り子を支えるバネ）、３枚の電気容量板（極板）を備える変位検出器４００４を備え、これらは容器４００５に収納されている。図５に示すとおり、振り子４００１は、変位検出器４００４に含まれる３枚の容量板のうち、中央の容量板に接続されている。また振動センサ４は、振り子の変位−電圧変換増幅器Ａ_D（便宜上、増幅器Ａ_Dの変換係数をＡ_D［Ｖ／ｍ］で表す。），微分器Ｇ_S（便宜上、微分器Ｇ_Sの係数をＧ_Sで表す。微分器Ｇ_Sの作用は、数式では

【数21】

で表される。），インピーダンス要素Ｒ_I（抵抗Ｒ_I［Ω］で表す。），インピーダンス要素Ｒ_D（抵抗Ｒ_D［Ω］で表す。），インピーダンス要素Ｒ_V（抵抗Ｒ_V［Ω］で表す。）を備え、また積分記号

【数22】

で表される積分器（Ｔ（ｓ）は時定数とする。積分器の作用は、数式では

【数23】

で表される。）を備える。なお、図５中では明示していないが、振り子４００１が接続されるサーボコイルのまわりには磁石が配置されており（図１０を参照。）サーボコイルは磁界の中で運動する。

【0031】

振動センサ４が設置されている地面等が振動した場合に対応して、振動センサ４に入力変位ｗ（ｔ）が与えられると（単位はｍとする。本明細書において、単位について明記しない量はＭＫＳＡ単位系で表される量であるとする。）、振り子４００１が振動する（格納容器４００５内の振り子の、容器に対する相対変位をｘ（ｔ）（単位はｍとする。）とする）。振り子４００１の変位により、変位検出器４００４に含まれる３枚の容量板のそれぞれの間隔が変化する。具体的に、静止状態において、上の容量板と中央の容量板の間隔と、中央の容量板と下の容量板の間隔とが、それぞれｄ［ｍ］であったとすると、図６中に示す方向で振り子４００１がｘだけ変位することにより、振り子４００１に接続された中央の容量板も同様に変位する。これにより、上の容量板と中央の容量板の間隔はｄ−ｘとなり、中央の容量板と下の容量板の間隔はｄ＋ｘとなる。したがって、上の容量板と中央の容量板の間の電気容量と、中央の容量板と下の容量板との間の電気容量がそれぞれ変化し、振り子４００１の変位ｘが電気回路Ａ_Dによって変位に比例した電圧として検出される。この電圧に対して上述の回路要素により微積分等を行い、サーボコイルに電流ｉ（ｔ）を流し、振り子４００１が静止するように制御する。すなわち、磁石（図１０参照。）による磁界中のサーボコイルに電流ｉ（ｔ）が流れる（フィードバック）ことにより力が発生し、この力により振り子が静止するように制御される。具体的には、入力変位ｗ（ｔ）に対して、

【数24】

を復元加速度として発生させ（Ｇ［Ｎ／Ａ］はサーボコイルのモータジェネレータ定数（ｍｏｔｏｒ ｇｅｎｅｒａｔｏｒ ｃｏｎｓｔａｎｔ）であり、ｍ［ｋｇ］は振り子４００１の質量である。）、復元加速度

【数25】

を実現する（ｗの上に「・・」を付しているのは、時間による２回の微分を意味する。上に「・」を付す場合には時間による１回の微分を意味する。時間変化するｗ以外の量についても、適宜同様に時間微分を表記する。）。この制御に要した電流から、加速度、速度、変位がセンサ値として得られる。具体的には、微分器Ｇ_Sの出力端から振り子４００１の加速度出力が得られ、変位−電圧変換増幅器Ａ_Dの出力端から振り子４００１の速度出力が得られ、積分器の出力端から振り子４００１の変位出力が得られる。

【0032】

次に、仮想センサであるＦＳＦセンサによって図５に示すＶＢＢ型地震計の動作を模擬するためのモデルを、非特許文献４にならって説明する。まず、仕様として、機構部の振り子系のパラメータとＦＳＦセンサの設計目標を以下のように与える。
［機構部の振り子系］
振り子の質量＝ｍ［ｋｇ］，減衰定数＝ｈ［無次元量］，固有円振動数＝ω₀［ｒａｄ］、サーボコイルのモータジェネレータ定数＝Ｇ［Ｎ／Ａ］
［ＦＳＦセンサの特性］
速度出力感度Ｓ_V［Ｖ／（ｍ／ｓ）］，低域遮断円振動数＝ω_C［ｒａｄ］，及び、

【数26】

【0033】

上記仕様に基づき、ＦＳＦセンサは４段階のステップで実現される。まずステップ１では、フィードバックゲインベクトル（列ベクトル）

【数27】

を下記のとおり決定する。ここでは、発振防止用の高域遮断円振動数ω₃を制御パラメータとして試みに与えて計算する。通常、安全側に見て１０³［Ｈｚ］程度とする。
［ステップ１：フィードバックゲインベクトル］

【数28】

【数29】

【0034】

次に、ステップ２では、加速度及び変位信号の感度、Ｓ_A及びＳ_Dを決定する。このとき、積分器の時定数Ｔと微分器の係数Ｇ_Sが制御パラメータとして与えられる。逆に、Ｓ_A及びＳ_Dを与えて、Ａ_D，Ｔ，Ｇ_Sを求めても問題はない。
［ステップ２：感度（Ｓｅｎｓｉｔｉｖｉｔｙ）］

【数30】

【0035】

最後に、ステップ３で、図５のフィードバックインピーダンスＲ_V，Ｒ_D，Ｒ_Iを決定する。
［ステップ３：フィードバックインピーダンス］

【数31】

【0036】

ステップ４は確認作業であり、得られたＦＳＦセンサの伝達関数の計算である。ここで、Ｈ_A（ｓ），Ｈ_V（ｓ），Ｈ_D（ｓ）が加速度、速度、及び、変位入力に対する、ＦＳＦセンサの加速度、速度、及び、変位出力に関する伝達関数となる。これらの伝達関数の極構造は、共通の特性方程式Δ（ｓ）＝０から求められる。
［ステップ４：伝達関数（Ｔｒａｎｓｆｅｒ ｆｕｎｃｔｉｏｎ）］

【数32】

【0037】

上記ステップで構築されたＦＳＦセンサのブロック線図（ｂｌｏｃｋ ｄｉａｇｒａｍ）は、状態ベクトル：

【数33】

を用いれば、図７に示すとおりとなる。図７において、ＦＳＦセンサの出力端子Ｄ，Ｖ，Ａでは、

【数34】

が出力として得られる（この時点では観測ノイズを導入していないことに留意する。）。また、復元加速度

【数35】

は、状態フィードバックゲイン

【数36】

を用いれば、

【数37】

（内積演算）となる（これが全状態帰還の意味である。）。その詳細は、図８に示すとおりである。図８は図７の復元加速度ｂ（ｔ）と状態ベクトルｘ（ｔ）の部分のみを抜き出した図であり、式１８を視覚化したものである。

【0038】

図７及び図８のブロック線図を数式表現、すなわち状態空間表現すれば、以下の式のとおりとなる。

【数38】

【数39】

【0039】

観測ベクトルにおいて速度出力のみの観測をするとすれば、（式２０）は次式となる。

【数40】

【0040】

ＦＳＦセンサの設計例として、ＶＢＢ型を対象とすると以下のようになる。まず振り子の質量をｍ＝０．０９［ｋｇ］，固有振動数をｆ₀＝１．５［Ｈｚ］，減衰定数を０．４，モータジェネレータ定数Ｇ＝５６［Ｎ／Ａ］とする換振器要素を考える。次に、ＶＢＢの特性として、Ｓ_V＝７５０［Ｖ／（ｍ／ｓ）］及び、

【数41】

とする。この条件で、ｆ₃＝１．０６２［ｋＨｚ］を与えると、フィードバックゲインが

【数42】

となり、

【数43】

が得られる。更に、

【数44】

とすれば、加速度出力の感度

【数45】

と、変位出力の感度

【数46】

が得られる。設計されたＦＳＦセンサの加速度入力に対する加速度、速度及び変位出力の周波数応答特性は、図９に示すとおりとなる。これは、標準的なＶＢＢの特性である。ＶＢＢの実機と同様に変位（Ｄ）端からの出力は加速度入力に対して直流感度を有する特性となる。図中のゲイン特性のｄＢ表示では、

【数47】

を０［ｄＢ］としている。

【0041】

以上説明した高感度速度計のモデルをまとめれば、以下のとおりである。
状態ベクトルを

【数48】

としたとき、

【数49】

【数50】

ここで、

【数51】

は観測ノイズを含まない量である。

【0042】

（式２８）〜（式３０）は連続時間系の表現であるが、標本化時間ΔＴをパラメータとして離散系に変換すると、

【数52】

で

【数53】

の条件下で次式となる。

【数54】

ここで

【数55】

及び

【数56】

【数57】

【0043】

（式３３）〜（式３６）の線形モデルに対して（システムノイズも観測ノイズもこの段階では導入していない。システムのブロック線図は図６を参照。）、（式２）を用いて既に説明したシステムノイズ

【数58】

を導入する（平均０、共分散行列Ｑ（ｎ）の正規性雑音とする。）。さらに、

【数59】

に観測ノイズｒ（ｎ）を加えた量を、以下のとおり新たにｙ₂（ｎ）と定義する。

【数60】

観測ノイズｒ（ｎ）は、平均０，分散Ｒ（ｎ）の正規性雑音であるとする。上記（式３７）が、観測ノイズがある場合の高感度速度計の信号モデルである。

【0044】

以上から、図２のセンサフュージョンの構成は、システムノイズ、観測ノイズを含む以下の線形システム

【数61】

により記述される。

【0045】

このモデルにおいて、（式３８）を

【数62】

と記述したとき、

【数63】

であるから、
ＶＢＢがクリップしない範囲では、

【数64】

が充分な精度で得られているとすると、

【数65】

となるため、

【数66】

となる。ＶＢＢがクリップする範囲では、

【数67】

ならば

【数68】

となるが、

【数69】

ではなく、

【数70】

を意味する。ＶＢＢで観測された信号は、このような性質の観測雑音ｒ（ｎ）によってモデル化される。（式３８）では、更に解釈を拡大して、ｒ（ｎ）を平均零、分散Ｒ（ｎ）の時変正規性雑音としてモデル化している。

【0046】

カルマンフィルタによる状態推定
（式３８）のシステムに、制御入力がある場合の線形カルマンフィルタを適用可能である。この場合のカルマンフィルタの予測ステップ、フィルタリングステップは、以下のとおりである。

【数71】

【0047】

（式４７）中、

【数72】

は、状態ベクトル

【数73】

の事前状態推定値であり、

【数74】

は、状態ベクトル

【数75】

に関する共分散行列（誤差共分散行列）

【数76】

の事前値である、事前共分散行列（事前誤差共分散行列）である。

【数77】

は、離散時刻ｎ−１における状態ベクトル

【数78】

の推定値であり、

【数79】

は、離散時刻ｎ−１における共分散行列

【数80】

の、カルマンフィルタにより演算される値であり、事後共分散行列である。

【数81】

は、離散時刻ｎにおける状態ベクトル

【数82】

の推定値であり、

【数83】

は、離散時刻ｎにおける共分散行列

【数84】

のカルマンフィルタにより演算される値であり、事後共分散行列である。また、

【数85】

は単位行列（現在の実施形態においては３行３列の行列）である。

【0048】

以上のとおり、システムノイズと観測ノイズとを導入することにより、カルマンフィルタを用いたセンサフュージョン（ＦＳＦセンサ、或いは高感度換振器の振り子の状態推定）が可能となる。

【0049】

振動センサシステムの動作フロー
次に振動センサシステム１の全体的な動作フローを設計時と運用時とに分けて説明する。

【0050】

設計時
図１１は、振動センサシステムの設計時の動作フローを示すフローチャートである。ステップＳ１１０１において、推定装置２の処理部２０９は、記憶部２１０に記憶された疑似的センサ作成プログラムを実行することにより、上記（式６）〜（式３６）を用いて説明した、高感度換振器と等価特性を有する（全）状態帰還型換振器のパラメータ計算を行い模擬的センサを作成する。これにより、（式３３）〜（式３６）で記述される状態帰還型換振器の離散化モデルが生成される（ステップＳ１１０２）。このモデルにおいて用いられる各種パラメータ等のデータは、推定装置２によるカルマンフィルタを用いた処理において用いられるので、各種パラメータ等のデータは推定装置２の記憶部２１０に記憶される（図４参照）。

【0051】

引き続き、推定装置２の入出力部２１１は、外部測定器３から低感度換振器離散信号（加速度記録）を、そして振動センサ４から高感度換振器離散信号（速度記録。）を、それぞれ時々刻々と取得し（ステップＳ１１０３）、これら取得した記録のデータを記憶部２１０に記憶させる。模擬的センサ（処理部２０９）は、（式３７）に記述されるとおり、高感度換振器離散信号に観測雑音モデルからの観測信号を加え（ステップＳ１１０４）、またシステムノイズも含む（式３８）の線形システムを生成する。以降、模擬的センサは、外部測定器３からの加速度記録を時々刻々と取得し、システムノイズ、観測ノイズをそれぞれの雑音モデル（予め定義されて記憶部２１０に記憶されているものとする。後述のとおり、雑音モデルは事後的に調整可能。）から時々刻々と生成し、離散的な各々の時刻ｎについて、（式３８）のモデルに従い高感度速度記録ｙ₂（ｎ）を生成し続ける。なお、（式３８）中の状態ベクトル

【数86】

における、ｎ＝０での初期値は、予め与えられて記憶部２１０内の状態推定プログラム用データとして記憶されているものとする（システムノイズ等、その他に初期値が必要な変数についても同様。）。

【0052】

次に、推定装置２の処理部２０９は、取得したそれら記録データと、記憶部２１０に記憶された各種データを用いて状態推定プログラムを実行することにより、（式４７）に従い、外生入力を持つカルマンフィルタによる高感度換振器の振り子運動推定を行い、状態ベクトル

【数87】

の推定値

【数88】

の演算と、推定誤差の事後共分散行列

【数89】

の、カルマンフィルタによる値

【数90】

の演算とを、時々刻々と行う（ステップＳ１１０５）。なお、（式４７）中のカルマンゲインにおける、

【数91】

は、観測ノイズの分散であり、外部から設定可能な量（「カルマンゲイン調整項」）である（一例においては、使用者が入出力部２１１を介して入力し、状態推定プログラムを実行する処理部２０９によって、記憶部２１０に状態推定プログラム用データとして記憶される。）。状態推定プログラムを実行する処理部２０９は、振動センサ４から入力される速度記録の示す速度が所定値（一例においては、使用者が入出力部２１１を介して入力し、状態推定プログラムを実行する処理部２０９によって、記憶部２１０に状態推定プログラム用データとして記憶される。）よりも大きい場合には（高感度速度計の速度記録が飽和している場合など）、所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくすることで、カルマンゲインを小さくすることができる。

【0053】

状態推定プログラムを実行する処理部２０９による、上記カルマンフィルタによる推定により、カルマンフィルタから推定された、高感度換振器の振り子運動から計算される状態帰還型換振器の出力が得られる（Ｓ１１０６）。この出力は、高感度（振動センサ４）及び低感度（外部測定器３）の融合信号であり（Ｓ１１０７）、推定された状態ベクトルの各成分を入出力部２１１のディスプレイ装置によって表示することにより、使用者は推定された振り子の状態を知ることができる。

【0054】

さらに設計時においては、状態推定プログラムを実行する処理部２０９が、外部測定器３、振動センサ４からの信号により示される加速度、速度等の記録値と、カルマンフィルタを用いて推定された状態ベクトルから演算される加速度、速度等の値との比較判定を行う（ステップＳ１１０８）。一例においては、ある離散時刻において、外部測定器３、振動センサ４からの信号により示される記録値と、カルマンフィルタを用いて推定された状態ベクトルから演算される値との間の差異が所定値を超える場合、状態推定プログラムを実行する処理部２０９は、記録値と演算値との一致性が低いと判断し（ステップＳ１１０８の条件分岐における「ＮＯ」）、観測ノイズの分散の大きさを変更するなどして観測雑音モデルの調整を行う（ステップＳ１１０９）。調整後の観測雑音モデルからの観測雑音を用いて、ステップＳ１１０４〜ステップＳ１１０８までの処理が再度行われる。ステップＳ１１０８において一致性が高いと判断されるまで、これらの処理が繰り返される。ステップＳ１１０８において一致性が高いと判断されると（一例においては、外部測定器３、振動センサ４からの信号により示される記録値と、カルマンフィルタを用いて推定された状態ベクトルから演算される値との間の差異が所定値を超えない場合。ステップＳ１１０８の条件分岐における「ＹＥＳ」）、設計は終了する。

【0055】

運用時
図１１を用いて説明した設計処理により、（式３８）の線形システム、（式４７）のカルマンフィルタが、各種演算パラメータや雑音モデルを含めて決定される。運用時には、これらを用いて、振動センサ４に関する状態ベクトルの推定が離散的な時刻について時々刻々と行われる（ｎ＝１，２，…とインクリメントしながら推定を繰り返す。）。運用時のフローチャートは図１２に示すとおりである。ステップＳ１２０１，Ｓ１２０２，Ｓ１２０３，Ｓ１２０４，Ｓ１２０５は、それぞれ、設計時のステップＳ１１０３，Ｓ１１０４，Ｓ１１０５，Ｓ１１０６，Ｓ１１０７と同様であってよい。

【0056】

カルマンフィルタによる演算の詳細
図１１のフローチャート中のステップＳ１１０５，図１２のフローチャート中のステップＳ１２０３において行われるカルマンフィルタによる演算の詳細なフローを、図１３のフローチャートに示す。

【0057】

まず、状態ベクトル事前推定値演算部２０４は、記憶部２１０から、前回の状態ベクトル推定値データ、前回の事後共分散行列データ、前回の測定値データ、演算パラメータの読み出しを行う（ステップＳ１３０１）。読み出したデータは処理部２０９のＲＡＭに一時記憶される（ＲＡＭへの各種データの一時記憶は、以降の記載、そしてこれまでの記載において適宜省略する。記憶部からの各種データの読み出しも適宜記載を省略する。）。初回の推定時、すなわち（式４７）においてｎ＝１として処理を行う場合、前回の状態ベクトル推定値データ、前回の事後共分散行列データ、前回の測定値データ（ｎ＝０におけるそれぞれのデータ）とは、それぞれ事前に使用者が入出力部２１１を介して記憶部２１０に記憶させる、外部から与えられたそれぞれの初期値データである。

【0058】

引き続き、状態ベクトル事前推定値演算部２０４は、（式４７）中、予測ステップの式

【数92】

に従い、状態ベクトルの事前推定値

【数93】

を演算し（ステップＳ１３０２）、処理部２０９のＲＡＭに一時記憶する。

【0059】

次に、事前共分散行列演算部２０５は、（式４７）中、予測ステップの式

【数94】

に従い、事前共分散行列

【数95】

を演算し（ステップＳ１３０３）、処理部２０９のＲＡＭに一時記憶する。

【0060】

次に、カルマンゲイン演算部２０６は、（式４７）中、フィルタリングステップの式

【数96】

に従い、カルマンゲインを演算し（ステップＳ１３０４）、処理部２０９のＲＡＭに一時記憶する。ここで、カルマンゲイン演算部２０６は、カルマンゲイン調整項

【数97】

を用いた演算によりカルマンゲインを調整する。具体的に、図１１のフロー中のステップＳ１１０３、又は図１２のフロー中のステップＳ１２０１で取得した高感度換振器離散信号の値（ｖ（ｎ））が所定値よりも大きい場合には、その値が所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくすることでカルマンゲインを小さくすることにより、カルマンゲインを調整する。振動センサ４からのセンサ値が、速度記録以外に変位記録、加速度記録である場合も、同様にカルマンゲイン演算部２０６により、各々の所定値と当該センサ値との比較結果に応じてカルマンゲイン調整項を介してカルマンゲインが調整される（振動センサ４からのセンサ値としての振り子の変位、速度、加速度のうち１以上の値が所定値よりも大きい場合には、１以上の値が所定値よりも小さい場合に比べてカルマンゲイン調整項を大きくする）。

【0061】

次に、状態ベクトル推定値演算部２０７が、（式４７）中、フィルタリングステップの式

【数98】

に従い、状態ベクトル推定値

【数99】

を演算し（ステップＳ１３０５）、処理部２０９のＲＡＭに一時記憶する。

【0062】

次に、事後共分散行列演算部２０８が、（式４７）中、フィルタリングステップの式

【数100】

に従い、事後共分散行列

【数101】

を演算し（ステップＳ１３０６）、処理部２０９のＲＡＭに一時記憶する。

【0063】

次に、ステップＳ１３０７において、状態ベクトル推定値演算部２０７は、ステップＳ１３０５で推定された状態ベクトル推定値を記憶部２１０に記憶し（離散時刻ｎにおける、この状態ベクトル推定値は、離散時刻ｎ＋１のカルマンフィルタ処理において、「前回の状態ベクトル推定値データ」として用いられる。）、また、事後共分散行列演算部２０８は、ステップＳ１３０６で演算した事後共分散行列を記憶部２１０に記憶する（離散時刻ｎにおける、この事後共分散行列は、離散時刻ｎ＋１のカルマンフィルタ処理において、「前回の事後共分散行列データ」として用いられる。）。その他、状態推定プログラムを実行する処理部２０９は、外部測定器３から取得した測定値（離散時刻ｎにおける、この測定値は、離散時刻ｎ＋１のカルマンフィルタ処理において、「前回の測定値データ」として用いられる。）等、種々のデータを記憶部２１０に記憶させる。ステップＳ１３０１〜Ｓ１３０７の処理は、ステップＳ１１０３，ステップ１２０１のデータ信号取得とともに時々刻々と繰り返され（データ取得を先にまとめて行い、その後にカルマンフィルタ処理をまとめて行ってもよいし、ある離散時刻についてデータ取得をした直後にその離散時刻についてカルマンフィルタ処理を行い、離散時刻に１を加えてから（インクリメント）、次の離散時刻についてデータ取得を行い、カルマンフィルタ処理を行うという流れで、データ取得とカルマンフィルタ処理とを同期させてもよい。）、低感度換振器離散信号と高感度換振器離散信号とが取得されなくなるか、或いは離散時刻が所定の終了値に達する等、所定の終了条件が満たされると、カルマンフィルタによる処理は終了する。

【0064】

図１３で示されるフローのカルマンフィルタ処理を含め、図１２で示される運用時のフローに従う処理が時々刻々と繰り返されることにより、時々刻々と変化する、振動センサ４に関する状態ベクトル（振り子の状態ベクトル）が推定される。状態ベクトル推定値は、状態推定プログラムを実行する処理部２０９により、入出力部２１１のディスプレイ装置に表示される等し、これにより使用者は状態ベクトル推定値を確認することができる。

【0065】

（第２の実施形態）
高感度変位計のモデル
次に、本発明の第２の実施形態として、振動センサ４として高感度変位計を用いる場合の実施形態を説明する。ただし、第１の実施形態と同様の部分については説明を省略する。

【0066】

図１４は、振動センサ（高感度変位計としてのＶＢＢ型地震計）の一例を示す図である。図１５は、高感度変位計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図である。高感度変位計は、図５の速度計と同じ構成で実現される。

【0067】

図１４の振動センサにおいて、変位信号は図１４の＃１の位置から取り出せるが、＃２の位置のように測定帯域外で減衰特性を有する外部積分器を通して得るのが一般的である。測定帯域内での連続時間系状態空間表現は、図１６に示すとおりとなるが、変位信号

【数102】

とともに、速度信号

【数103】

も同時に出力されるため、これらを併せて利用することが可能である。

【0068】

このとき、状態空間表示は、振動センサ４に関する振り子の状態ベクトルを

【数104】

とすると、

【数105】

【数106】

となる。この時間連続系の表現は、標本化時間ΔTをパラメータとして離散系に変換すると、

【数107】

で

【数108】

の条件下で次式となる：

【数109】

【数110】

【0069】

したがって、平均ベクトルがゼロベクトル（２次元の列ベクトル）、共分散行列が

【数111】

の正規性雑音

【数112】

を考えると、高感度変位計の信号モデルは、

【数113】

となり、以下のカルマンフィルタ

【数114】

によるセンサフュージョンが可能となる。振動センサシステム１の装置構成、設計時、運用時の動作フローは、システムノイズ

【数115】

を第１の実施形態と同様に導入した上で、線形システムの式として

【数116】

を用い、カルマンフィルタの式として（式６２）の式を用いる以外は、第１の実施形態と同様に実現できる。

【0070】

（第３の実施形態）
高感度加速度計のモデル
次に、本発明の第３の実施形態として、振動センサ４として高感度加速度計を用いる場合の実施形態を説明する。ただし、第１の実施形態と同様の部分については説明を省略する。

【0071】

図１７は、振動センサ（高感度加速度計）の一例を示す図である。図１８は、高感度加速度計と強震計のカルマンフィルタを用いたセンサフュージョンの構成を示すブロック線図である。高感度加速度計は、図５の速度計と同じ構成でも実現可能である。

【0072】

図５のＶＢＢ型加速度計の場合、加速度信号は図７のＡ端子から

【数117】

として出力される。高感度速度計の

【数118】

を、

【数119】

で置き換えれば、センサフュージョンが行える。しかしながら、図５の構成の地震計は周波数０での感度も零となり、常用される力平衡型の負帰還型加速度計と比較して不利となる。そこで、ここでは、高感度加速度計として、図１７に示す通常の負帰還型加速度計を扱うことにする。図１７の加速度計に対する連続時間系の状態空間表示は図１９となる。

【0073】

図１９の状態空間表示は、振動センサ４に関する振り子の状態ベクトルを

【数120】

としたとき、

【数121】

【数122】

となる。この時間連続系の表現は、標本化時間ΔＴをパラメータとして離散系に変換すると、

【数123】

で

【数124】

の条件下で次式となる：

【数125】

【数126】

【0074】

したがって、平均０、分散が

【数127】

の正規性雑音

【数128】

を考えると、高感度加速度計の信号モデルは、

【数129】

となり、図１８に示すカルマンフィルタ

【数130】

によるセンサフュージョンが可能となる。振動センサシステム１の装置構成、設計時、運用時の動作フローは、システムノイズ

【数131】

を第１の実施形態と同様に導入した上で、線形システムの式として

【数132】

を用い、カルマンフィルタの式として（式７２）の式を用いる以外は、第１の実施形態と同様に実現できる。

【産業上の利用可能性】

【0075】

本発明は、地震計をはじめとするあらゆる振動センサに関する状態推定に利用することができる。

【符号の説明】

【0076】

１ 振動センサシステム
２ 推定装置
３ 外部測定器（低感度換振器）
４ 振動センサ（高感度換振器）
２００ カルマンフィルタ演算部
２０１ 測定値取得部
２０２ データ格納部
２０３ センサ値取得部
２０４ 状態ベクトル事前推定値演算部
２０５ 事前共分散行列演算部
２０６ カルマンゲイン演算部
２０７ 状態ベクトル推定値演算部
２０８ 事後共分散行列演算部
２０９ 処理部
２１０ 記憶部
２１１ 入出力部
４００１ サーボコイルに接続された振り子
４００２ 減衰要素
４００３ 剛性要素
４００４ 変位検出器（３枚の容量板）
４００５ 容器

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】