JP6132990B2

JP6132990B2 - 状態推定装置

Info

Publication number: JP6132990B2
Application number: JP2016559757A
Authority: JP
Inventors: 哲太郎山田; 良晃横山; 小幡　康; 康小幡
Original assignee: Mitsubishi Electric Corp
Current assignee: Mitsubishi Electric Corp
Priority date: 2014-11-20
Filing date: 2014-11-20
Publication date: 2017-05-24
Anticipated expiration: 2034-11-20
Also published as: US9964408B2; DE112014007193B4; US20170314920A1; JPWO2016079848A1; MY167374A; WO2016079848A1; DE112014007193T5

Description

この発明は、海面の状態推定を行う状態推定装置に関するものである。

遠方で津波を検知し早期に警戒を促すため、可能な限り早く且つ精度よく津波を予測することが求められている。そのため、正確な津波の波高値と正確な流速値をリアルタイムに知る必要がある。この対策としては、以下のような発明が先行技術として存在している。（例えば非特許文献１−４及び特許文献１参照。）

非特許文献１に開示された技術は、津波の運動モデルである非線形浅水方程式を用いて、レーダから得られる流速値から波高値の予測をリアルタイムに行うものである。
また、非特許文献２に開示された技術は、沿岸から沖合に向かって接地されたＧＰＳブイ等の津波計の波高値から観測誤差を考慮して波高値のフィルタリングを行うことで、波高値の予測をリアルタイムに行うものである。
また、特許文献１に開示された技術は、観測位置において観測される津波波高に基づいて津波波源を推定し、推定された津波波源から津波波高を推定するものである。

また、非特許文献３に開示された技術では、ＨＦレーダを用いた津波の検出と、浅水方程式による津波の流速と波高の関係式が示されている。
また、非特許文献４に開示された技術では、津波が起きていない平常時の流速の推定方法として、ＨＦレーダの流速観測値を潮汐のモデルに同化する方法が示されている。

泉宮，"海洋レーダによる海表面流速場を用いた津波のリアルタイム予測に関する研究"，海岸工学論文集，第52巻，PP.546-550. 小池，"カルマンフィルタを用いたリアルタイム津波予警報に関する研究"，土木学会論文集 No.27-0293，2003年． BELINDA J. LIPA，DONALD E. BARRICK，JOHN BOURG and BRUCE B. NYDEN"HF Radar Detection of Tsunamis"，Journal of Oceanography，Vol. 62，pp. 705 to 716，2006 Oyvind Breivik, Oyvind Saetra，"Real time assimilation of HF radar currents into a coastal ocean model, Oceans"，Journal of Marine Systems，Volume 28，Issues 3-4，pages 161-182，April 2001.

特開２００８−０８９３１６号公報

しかしながら、非特許文献１に開示された技術において、レーダによる流速値には観測誤差が含まれるが、この観測誤差を考慮していない。そのため、波高予測の精度が劣化するという課題がある。また、２台以上のレーダの設置が前提となる（２次元空間上の流速ベクトルが求められることが前提となる）という課題もある。

また、非特許文献２に開示された技術において、様々な方位より由来する可能性のある津波を検知するためには、津波計を２次元空間上に複数台設置する必要があり、津波計を設置するコストが大きいという課題がある。また、ＧＰＳブイ等の津波計では、流速値を直接観測することができないため、波高値から津波の運動方程式に従って算出した流速値は誤差が大きいという課題もある。また、ＧＰＳブイを用いた津波計は、陸地から遠方の沖合には設置できないため、近傍の海の波高値しかわからないという課題もある。
また、特許文献１に開示された技術についても、沖合の波高値を直接観測できることが前提となっているため、非特許文献２と同様の課題がある。

また、非特許文献３記載の浅水方程式を元に流速値から波高値を計算する場合、非特許文献１と同様に流速値の誤差を考慮していないため、波高値の推定精度が劣化するという課題がある。
また、非特許文献４では津波の推定方式について記載がなく、また、海面の波高値の推定方法についても記載がない。

この発明は、上記のような課題を解決するためになされたもので、正確な津波の波高値と正確な流速値をリアルタイムに知ることができる状態推定装置を提供することを目的としている。

この発明に係る状態推定装置は、ビームによる海面の流速値の計測、及び沿岸の波高値の計測を行う観測部と、ビームのレンジセル毎の流速値、当該レンジセルの境界間の波高差及び沿岸の波高値から成る状態ベクトルから、次時刻における当該状態ベクトルを予測する予測部と、平滑誤差共分散行列から予測誤差共分散行列を算出する予測誤差共分散行列算出部と、観測部及び予測誤差共分散行列算出部による処理結果から、ゲイン行列を算出するゲイン行列算出部と、観測部、ゲイン行列算出部及び予測誤差共分散行列算出部による処理結果から、平滑誤差共分散行列を算出する平滑誤差共分散行列算出部と、観測部、予測部及びゲイン行列算出部による処理結果から、波高差毎に状態ベクトルを平滑化する波高差平滑部と、波高差平滑部による処理結果を用い、沖合方向に波高値と波高差を加算していくことで、レンジセル毎の波高値を算出する波高値算出部とを備えたものである。

この発明によれば、上記のように構成したので、正確な津波の波高値と正確な流速値をリアルタイムに知ることができる。

この発明の実施の形態１に係る状態推定装置の構成を示す図である。この発明の実施の形態１に係る状態推定装置におけるレーダと津波の様子を示す図である。この発明の実施の形態１に係る状態推定装置におけるレーダ覆域と状態ベクトルを示す図である。この発明の実施の形態１に係る状態推定装置における距離と波高差との関係を示す図である。従来の状態推定装置における平滑化処理の様子を示す図である。この発明の実施の形態１に係る状態推定装置における平滑化処理の様子を示す図である。

以下、この発明の実施の形態について図面を参照しながら詳細に説明する。
実施の形態１．
図１はこの発明の実施の形態１に係る状態推定装置の構成を示す図である。
状態推定装置は、図１に示すように、観測部１、初期波高値算出部２、平滑誤差共分散行列初期値設定部３、平滑値ベクトル初期値設定部４、予測部５、波高差平滑部６、ゲイン行列算出部７、平滑誤差共分散行列算出部８、予測誤差共分散行列算出部９、波高値算出部１０、尤度算出部１１及び警報部１２から構成されている。なお、状態推定装置の各機能部は、ソフトウェアに基づくＣＰＵを用いたプログラム処理によって実行される。

観測部１は、海面の流速値及び沿岸の波高値を計測するものである。ここで、海面の流速値はレーダを用いてビームによりセル単位で計測し、沿岸の波高値は沿岸に設置された潮位計を用いて計測する。また、セルとは、距離及びビームの分解能単位を表す。この観測部１により計測された観測値（流速値と波高値）を示す情報は、初期波高値算出部２、平滑値ベクトル初期値設定部４、波高差平滑部６、ゲイン行列算出部７及び平滑誤差共分散行列算出部８に出力される。

初期波高値算出部２は、目標（海面）の状態ベクトルの波高差の初期値を算出するものである。この際、初期波高値算出部２は、観測部１により計測された流速値に基づいて、レンジセル毎に、初期時刻における流速値と次の時刻における流速値の２時刻分の流速値の差から、状態ベクトルの波高差の初期値を算出する。この初期波高値算出部２により算出された波高差を示す情報は、平滑値ベクトル初期値設定部４に出力される。

平滑誤差共分散行列初期値設定部３は、初期時刻における平滑誤差共分散行列の初期値を設定するものである。この平滑誤差共分散行列初期値設定部３により設定された平滑誤差共分散行列を示す情報は、予測誤差共分散行列算出部９に出力される。

平滑値ベクトル初期値設定部４は、観測部１により計測された観測値、及び初期波高値算出部２により算出された波高差を、目標の状態ベクトルの初期値として設定するものである。この平滑値ベクトル初期値設定部４により設定された目標の状態ベクトルを示す情報は、予測部５に出力される。

予測部５は、目標の状態ベクトルの予測を行うものである。この際、予測部５は、初期フェーズでは平滑値ベクトル初期値設定部４により設定された目標の状態ベクトルに対し、また、追尾フェーズでは波高差平滑部６により平滑化された目標の状態ベクトルに対し、次の時刻における目標の状態ベクトルを予測する。この予測部５により予測された目標の状態ベクトルを示す情報は、波高差平滑部６及び尤度算出部１１に出力される。

波高差平滑部６は、予測部５により予測された目標の状態ベクトルと、観測部１により計測された観測値、及びゲイン行列算出部７により算出されたゲイン行列に基づいて、波高差毎に目標の状態ベクトルを平滑化するものである。この波高差平滑部６により平滑化された目標の状態ベクトルを示す情報は、予測部５及び波高値算出部１０に出力される。

ゲイン行列算出部７は、観測部１により計測された観測値、及び予測誤差共分散行列算出部９により算出された予測誤差共分散行列から、ゲイン行列を算出するものである。このゲイン行列算出部７により算出されたゲイン行列を示す情報は、波高差平滑部６及び平滑誤差共分散行列算出部８に出力される。

平滑誤差共分散行列算出部８は、観測部１により計測された観測値、ゲイン行列算出部７により算出されたゲイン行列、及び予測誤差共分散行列算出部９により算出された予測誤差共分散行列から、平滑誤差共分散行列を算出するものである。この平滑誤差共分散行列算出部８により算出された平滑誤差共分散行列を示す情報は、予測誤差共分散行列算出部９に出力される。

予測誤差共分散行列算出部９は、予測誤差共分散行列を算出するものである。この際、予測誤差共分散行列算出部９は、初期化フェーズでは平滑誤差共分散行列初期値設定部３により算出された平滑誤差共分散行列を用いて、追尾フェーズでは平滑誤差共分散行列算出部８により算出された平滑誤差共分散行列を用いて、予測誤差共分散行列を算出する。この予測誤差共分散行列算出部９により算出された予測誤差共分散行列を示す情報は、ゲイン行列算出部７、平滑誤差共分散行列算出部８及び尤度算出部１１に出力される。

波高値算出部１０は、波高差平滑部６により平滑化された目標の状態ベクトルに基づいて、沖合方向に波高値と波高差を加算していくことで、レンジセル毎の波高値を算出するものである。この波高値算出部１０により算出された波高値、及びその際の平滑化された目標の状態ベクトルの流速値を示す情報は警報部１２に出力される。

尤度算出部１１は、予測部５により予測された目標の状態ベクトル、及び予測誤差共分散行列算出部９により算出された予測誤差共分散行列から、ビーム毎の尤度を算出し、当該尤度が最大となるビーム方向を津波の到来方向として検出するものである。この尤度算出部１１により検出されたビーム方向を示す情報は、警報部１２に出力される。

警報部１２は、尤度算出部１１により検出されたビーム方向、及び波高値算出部１０により算出された波高値とその際の流速値の大きさから、津波の警報を出力するものである。また、警報部１２は、波高値算出部１０により算出された波高値とその際の流速値を初期条件として、津波の運動方程式を用いて流体シミュレーションを行うことで、津波の到来時間を予測する。

次に、上記のように構成された状態推定装置の動作について説明する。
状態推定装置の動作は、初期化フェーズと追尾フェーズとに分けられる。初期化フェーズでは、２時刻分の観測値（流速値、波高値）を蓄積した上で、目標の状態ベクトルの初期値設定及び平滑誤差共分散行列の初期値設定を実行することで、カルマンフィルタの初期値を設定する。また、追尾フェーズでは、通常のカルマンフィルタと同様に、時刻毎に状態ベクトルの更新をオンラインで実行する。

状態推定装置の動作では、まず、観測部１は、ビームによる海面の流速値の計測、及び沿岸の波高値の計測を行う。図２は、レーダを用いた海面の流速値の計測を示す図である。

次いで、初期波高値算出部２は、観測部１により計測された流速値に基づいて、目標の状態ベクトルの波高差の初期値を算出する。

ここで、１台のレーダによって得られた流速値は、視線方向の流速成分しか検出できないため、ビーム毎にレンジ方向の波の相互作用を考慮したカルマンフィルタを構成する。このため、図３及び次式（１）に示すように、沿岸から沖合までの各レンジセルにおける流速値（流量）、レンジセルの境界間の空間方向の波高差、及び沿岸の波高値を目標の状態ベクトルとして、ビーム毎にカルマンフィルタを実行する。

式（１）において、Ｘ（ｋ）は時刻ｋにおける目標の状態ベクトル、Ｍ_ｉ（ｋ）は時刻ｋにおけるｉ番目のレンジセルの流量、η_０（ｋ）は沿岸の波高値、Δη_{ｉ−１，ｉ}（ｋ）はｉ番目のレンジセルの波高差を表す。

以下の処理では、ビームを固定した場合のカルマンフィルタの処理について示すが、実際にはビームの本数分だけカルマンフィルタを並列して実行することを想定している。

そして、初期波高値算出部２は、レンジセル毎に、初期時刻における流速値と次の時刻における流速値の２時刻分の流速値を入力とし、次式（２），（３）を用いて、初期時刻と次の時刻における流速値の差から状態ベクトルの波高差の初期値を算出する。

式（２），（３）において、ｕ_ｉはｉ番目のレンジセルの流速、ｔは時間、ｈは水深、ｄはレンジセル数、ｇは重力加速度、Δｘはレンジセル間距離、Δｔはサンプリング間隔を表す。また、水深ｈは、地形データ等から既知であると仮定する。

次いで、平滑誤差共分散行列初期値設定部３は、次式（４）〜（７）を用いて、初期時刻における平滑誤差共分散行列を設定する。

式（４）において、Ｐ_ｋ：ｋは時刻ｋにおける平滑誤差共分散行列、Ｒは観測誤差共分散行列、Ｑは駆動雑音共分散行列、Ｇは式（５）〜（７）で表される駆動雑音変換行列、Ｉｄはｄ×ｄのサイズの単位行列を表す。ここで、式（５）〜（７）では、津波が運動する際に波高差が正規分布に従って揺らぐことを想定しているが、波高差の時間微分が正規分布に従って揺らぐとして定式化を行ってもよい。

次いで、平滑値ベクトル初期値設定部４は、観測部１により計測された観測値、及び初期波高値算出部２により算出された波高差を、目標の状態ベクトルの初期値として設定する。

次いで、予測部５は、目標の状態ベクトルの予測を行う。この際、予測部５は、初期フェーズでは平滑値ベクトル初期値設定部４により設定された目標の状態ベクトルに対し、また、追尾フェーズでは波高差平滑部６により平滑化された目標の状態ベクトルに対し、次式（８）〜（１０）を用いて、次の時刻における目標の状態ベクトルを予測する。

式（８）において、Ｘ_{ｋ：ｋ−１}は時刻ｋにおける予測された状態ベクトルを表す。また、運動モデルを表す（２ｄ＋１）×（２ｄ＋１）のサイズの遷移行列Ａは式（９），（１０）で表される。

以下に、遷移行列Ａの詳細について説明する。
まず、１次元線形浅水方程式より、流速値と波高差との関係は、次式（１１），（１２）で表すことができる。

これに、境界条件に関する関係式（１３）を考慮して、レンジセル毎の１次元線形浅水方程式を運動モデルとして行列にまとめると、運動モデルを表す遷移行列Ａは、式（９），（１０）で表すことができる。ただし、上記では運動モデルを１次元線形浅水方程式から導いたが、非線形項を考慮した浅水方程式を利用して、拡張カルマンフィルタやパーティクルフィルタ等の非線形フィルタを用いて平滑化を行ってもよい。また、１次元線形浅水方程式は、水深が浅い地点では誤差が大きいことが知られているので、駆動雑音を観測部１による計測領域内での水深の深さに応じて設定してもよい。

次いで、波高差平滑部６は、予測部５により予測された目標の状態ベクトルと、観測部１により計測された観測値、及びゲイン行列算出部７により算出されたゲイン行列に基づいて、次式（１４），（１５）を用いて、波高差毎に目標の状態ベクトルを平滑化する。

式（１４）において、Ｘ_ｋ：ｋは時刻ｋにおける平滑化された状態ベクトル、Ｈは流量と流速との関係から式（１５）で表されるｄ×（２ｄ＋１）のサイズの観測行列、Ｚ（ｋ）は沿岸の波高値と流速値の観測値ベクトルを表す。

次いで、ゲイン行列算出部７は、観測部１により計測された観測値、及び予測誤差共分散行列算出部９により算出された予測誤差共分散行列から、次式（１６）を用いて、ゲイン行列を算出する。

式（１６）において、Ｈ^ｔは行列Ｈの転置を表す。

次いで、平滑誤差共分散行列算出部８は、観測部１により計測された観測値、ゲイン行列算出部７により算出されたゲイン行列、及び予測誤差共分散行列算出部９により算出された予測誤差共分散行列から、次式（１７）を用いて、平滑誤差共分散行列を算出する。

次いで、予測誤差共分散行列算出部９は、次式（１８）を用いて、予測誤差共分散行列を算出する。この際、予測誤差共分散行列算出部９は、初期化フェーズでは平滑誤差共分散行列初期値設定部３により算出された平滑誤差共分散行列を用い、また、追尾フェーズでは平滑誤差共分散行列算出部８により算出された平滑誤差共分散行列を用いて、予測誤差共分散行列を算出する。

式（１８）において、Ｐ_{ｋ：ｋ−１}は時刻ｋ−１における予測誤差共分散行列を表す。

次いで、波高値算出部１０は、波高差平滑部６により平滑化された目標の状態ベクトルに基づいて、次式（１９）を用いて、沖合方向に波高値と波高差を加算していくことで、レンジセル毎の波高値を算出する（図４参照）。

次いで、尤度算出部１１は、予測部５により予測された目標の状態ベクトル、及び予測誤差共分散行列算出部９により算出された予測誤差共分散行列から、次式（２０），（２１）を用いて、ビーム毎の尤度を算出し、当該尤度が最大となるビーム方向を津波の到来方向として検出する。

式（２０）において、Ｐは確率密度関数を表し、式（２１）で表される正規分布とする。また、式（２１）において、Ｎ（Ａ，ｂ）は平均Ａ、分散ｂの正規分布を表す。

次いで、警報部１２は、尤度算出部１１により検出されたビーム方向、及び波高値算出部１０により算出された波高値とその際の流速値の大きさから、津波の警報を出力する。また、波高値算出部１０により算出された波高値とその際の流速値を初期条件とし、津波の運動方程式を用いて流体シミュレーションを行うことで、津波の到来時間を予測する。

以上のように、この実施の形態１によれば、波高値を空間方向の波高差に分解し、流速値（流量）、沿岸の波高値及び波高差を目標の状態ベクトルとして時系列方向に平滑化した上で、平滑化を行った波高値と波高差を沖合方向に加算するように構成したので、正確な津波の波高値と正確な流速値をリアルタイムに知ることができる。

例えば非特許文献２に示すような波高値を状態ベクトルとする構成で、沖合の波高値を推定しようとすると、流速値は観測できるが沿岸の波高値以外の波高値が直接得られない状況において、正確に波高値を推定することができない。

図５は流量（流速値）と波高値の伝搬の様子と従来方式での沖合の波高値の平滑化の概念を示す図である。図５において、波高値と流量からできるひし形内の矢印は、１次元浅水方程式による波高値と流量の伝搬の様子を表す。
従来方式では、カルマンフィルタの初期値として、式（２），（１９）を用いて算出した初期波高値を、周囲の流量を用いて平滑化している。そのため、初期波高値の誤差の影響が大きく、時系列に平滑化を実施しても波高値の推定精度は改善されない。これに加え、沖合の波高値の平滑化に利用される流量が、図５の三角形の領域の流量に限られるため、波高値の推定精度が劣化する。

これに対して、本発明では、波高値を空間方向の波高差に分解して、波高差毎に平滑化を行い、空間方向の波高差の加算処理で波高値を求めている。これにより、初期波高値の推定精度が悪い場合においても、時刻が経過する毎に波高差が平滑化されて、波高値を精度よく推定できる。

図６は本発明での波高差の平滑化の概念を示す図である。
この図６に示すように、沿岸の波高値から一番近傍の波高差の平滑化に使用される流量（流速値）は、図中の破線内の流量である。そして、沖合の波高値は波高差の加算で表されることから、沖合の波高値の算出の際に利用される流量数は、全レンジセルの過去の全時刻の流量数となる。このことから、従来方式よりも波高値の推定精度が向上することが期待できる。

なお、本願発明はその発明の範囲内において、実施の形態の任意の構成要素の変形、もしくは実施の形態の任意の構成要素の省略が可能である。

この発明に係る状態推定装置は、正確な津波の波高値と正確な流速値をリアルタイムに知ることができ、海面の状態推定を行う状態推定装置等に用いるのに適している。

１観測部、２初期波高値算出部、３平滑誤差共分散行列初期値設定部、４平滑値ベクトル初期値設定部、５予測部、６波高差平滑部、７ゲイン行列算出部、８平滑誤差共分散行列算出部、９予測誤差共分散行列算出部、１０波高値算出部、１１尤度算出部、１２警報部。

Claims

ビームによる海面の流速値の計測、及び沿岸の波高値の計測を行う観測部と、
前記ビームのレンジセル毎の流速値、当該レンジセルの境界間の波高差及び沿岸の波高値から成る状態ベクトルに対し、次時刻における当該状態ベクトルを予測する予測部と、
平滑誤差共分散行列から予測誤差共分散行列を算出する予測誤差共分散行列算出部と、
前記観測部及び前記予測誤差共分散行列算出部による処理結果から、ゲイン行列を算出するゲイン行列算出部と、
前記観測部、前記ゲイン行列算出部及び前記予測誤差共分散行列算出部による処理結果から、平滑誤差共分散行列を算出する平滑誤差共分散行列算出部と、
前記観測部、前記予測部及び前記ゲイン行列算出部による処理結果から、波高差毎に前記状態ベクトルを平滑化する波高差平滑部と、
前記波高差平滑部による処理結果を用い、沖合方向に波高値と波高差を加算していくことで、前記レンジセル毎の波高値を算出する波高値算出部と
を備えた状態推定装置。
前記予測部は、浅水方程式による運動モデルを用い、次時刻における前記状態ベクトルを予測する
ことを特徴とする請求項１記載の状態推定装置。
前記予測部及び前記予測誤差共分散行列算出部による処理結果から、前記ビーム毎の尤度を算出し、当該尤度が高いビーム方向を津波の到来方向として検出する尤度算出部を備えた
ことを特徴とする請求項１記載の状態推定装置。
前記尤度算出部により検出されたビーム方向、前記波高値算出部により算出された波高値、前記波高差平滑部により平滑化された前記状態ベクトルの流速値から、津波の到来時間を予測する警報部を備えた
ことを特徴とする請求項３記載の状態推定装置。
前記観測部により計測された流速値に基づいて、前記レンジセル毎に、初期時刻における流速値と次時刻における流速値の２時刻分の流速値との差から、前記状態ベクトルにおける波高差の初期値を算出する初期波高値算出部を備え、
前記予測部は、初期フェーズでは、前記観測部により観測された前記レンジセル毎の流速値及び沿岸の波高値、及び前記初期波高値算出部により算出された前記レンジセルの境界間の波高差から成る状態ベクトルから、次時刻における当該状態ベクトルを予測する
ことを特徴とする請求項１記載の状態推定装置。
前記予測部は、非線形浅水方程式による運動モデルを用い、次時刻における前記状態ベクトルを予測する
ことを特徴とする請求項１記載の状態推定装置。
前記予測部は、前記浅水方程式による運動モデルにおいて、前記観測部による計測領域内での水深に応じて、駆動雑音の大きさを設定する
ことを特徴とする請求項２記載の状態推定装置。