JP6610898B2

JP6610898B2 - 横メタセンタ高さ推定装置及び横メタセンタ高さ推定方法

Info

Publication number: JP6610898B2
Application number: JP2016521124A
Authority: JP
Inventors: 正裕玉島; 大介寺田; 秋彦松田; 仁前野; 孝二林
Original assignee: FLUID TECHNO CO., LTD.; Furuno Electric Co Ltd; Japan Fisheries Research and Education Agency
Current assignee: FLUID TECHNO CO., LTD.; Furuno Electric Co Ltd; Japan Fisheries Research and Education Agency
Priority date: 2014-05-20
Filing date: 2015-05-20
Publication date: 2019-11-27
Anticipated expiration: 2035-05-20
Also published as: US20170088236A1; JPWO2015178410A1; WO2015178410A1

Description

参照による取り込み

本出願は、２０１４年５月２０日に出願された日本特許出願特願２０１４−１０４７８６の優先権を主張し、その内容を参照することにより、本出願に取り込む。

本発明は、横メタセンタ高さ推定装置及び横メタセンタ高さ推定方法に関する。

従来、不規則に変動する波浪上を航行する船舶では、安全性の観点からも、船体運動データ、船体状態データ及び海象条件を適切に把握することが重要である。船体運動データとは、船体の変位や加速度等の船体の運動に係るデータである。また船体状態データとは、船体の喫水、排水量、横メタセンタ高さ（以下、「ＧＭ」ともいう。）等の船体の状態に係るデータである。また海象条件とは、船舶の航行領域における波浪の波高、波周期、波向き等の海象に係る情報である。

従来、これらの情報のうちの船体運動データについては、予め設定された船体状態データと、気象庁等の情報提供機関から提供される海象条件とに基づいて力学的に求めるのが一般的であった。

しかしながら、この方法では適切に船体運動データや船体状態データを把握することができなかった。これは、気象庁等の情報提供機関から提供される海象条件が、広い海域に亘る情報量が膨大なものであり、現に船舶が航行している海域の局所的なものではなく且つ精度が低いものであったためである。

そこで近年では、船舶に搭載された各種装置によって船体運動データを計測し、計測された非定常な時系列データである船体運動データに対して統計処理をリアルタイムに施すことによって、船体状態データや海象条件を統計的に推定する方法が報告されている（例えば、非特許文献１参照）。

非特許文献１には、非定常な時系列データである船体運動データを、時変係数多変量自己回帰（ＴＶＶＡＲ：Time Varying Coefficient Vector AR）モデルを用いて解析することによって、海象条件を推定する技術が開示されている。

平山次清、井関俊夫、石田茂資、"遭遇波浪のリアルタイム把握法と最新の成果"、日本造船学会、平成１５年１２月、ｐ．７４−９６

さて、上記従来技術に示すように、非定常な時系列データとしての船体運動データに基づいて船体状態データを推定しているが、その一例として、船体運動データの一つである「横揺れ角の時系列データ」に基づいて、船体状態データの一つである「ＧＭ」を推定している。具体的には、まず横揺れ角の時系列データに基づいて横揺れ固有周波数を推定し、次に推定された横揺れ固有周波数に基づいてＧＭを推定している。

しかしながら、従来の推定方法では、精度良くＧＭを推定することができなかった。これは、横揺れ固有周波数に基づいてＧＭを推定する際に、式（１）に示すような近似式を用いているためである。

式（１）において、Ｔは横揺れ固有周期、Ｃは実験定数、Ｂは船の幅、ｄは船の喫水、Ｌは船の長さである。式（１）では、貨物の積載状態などによって異なる値を示すＣを、定数として取り扱っている。そのため、その分ＧＭの推定値の精度が低いものとなってしまっていた。

本発明は、上記のような事情に鑑みてなされたものであって、精度良く横メタセンタ高さを推定することが可能な横メタセンタ高さ推定装置及び横メタセンタ高さ推定方法を提供することを目的とする。

上記の目的を達成するために、本発明に係る横メタセンタ高さ推定装置は、船体の横揺れ角の時系列データを記憶する履歴記憶手段と、前記履歴記憶手段によって記憶された船体の横揺れ角の時系列データに基づいて、当該船体の横メタセンタ高さを推定する横メタセンタ高さ推定手段と、を有し、前記横メタセンタ高さ推定手段は、まず前記船体の横揺れ角の時系列データに基づいて横揺れ固有周波数を計算し、計算された横揺れ固有周波数を観測モデルとし、横メタセンタ高さと当該船体の慣動半径とを状態変数とする一般状態空間モデルに基づく状態推定により、前記横メタセンタ高さを推定することを特徴とする。

本発明によれば、精度良く横メタセンタ高さを推定することができる。

本実施形態に係る横メタセンタ高さ推定装置を含む海象推定システムのハードウェア構成例を示す図である。本実施形態に係る横メタセンタ高さ推定装置を含む海象推定システムの機能構成例を示す図である。本実施形態に係る横メタセンタ高さ推定装置を含む海象推定システムの制御ロジックの一例を示すフローチャートである。

以下、本発明の実施形態について説明する。

［システムのハードウェア構成］
図１は、本実施形態に係る横メタセンタ高さ推定装置を含む海象推定システムのハードウェア構成例を示す図である。

図１に示す海象推定システム１は、サテライトコンパス２、情報処理装置（横メタセンタ高さ推定装置）３、ディスプレイ４を備える。この海象推定システム１は船舶の船体に搭載される。

サテライトコンパス（ＧＰＳコンパス）２は、船体の船首方向に取り付けられた２つのＧＰＳアンテナの相対的な位置関係から船体の方位を計算する方位センサーとしての機能を有する装置である。このサテライトコンパス２は、船体の横揺れ（ロール）、縦揺れ（ピッチ）、上下揺れ（ヒーブ）を計測可能な動揺センサーとしての機能も有するものとする。なお、サテライトコンパス２の代わりに、ジャイロセンサを用いても良い。

情報処理装置３は、それぞれバス３８で相互に接続されたメモリ装置３１、演算処理装置３２、インターフェース装置３３、入力装置３４、補助記憶装置３５、ドライブ装置３６を備えるコンピュータ装置である。この情報処理装置３は、サテライトコンパス２によって計測された情報に基づいて横メタセンタ高さを推定する。また、横メタセンタ高さの推定値等に基づいて海象条件を推定する。この情報処理装置３が、特許請求の範囲の「横メタセンタ高さ推定装置」に対応する。なお、この情報処理装置３及び後述するディスプレイ４はサテライトコンパス２と一体的に構成されても良い。

メモリ装置３１は、情報処理装置３の起動時に補助記憶装置３５に記憶されたプログラム（図２の船体状態データ計算部２３や海象推定部２４の機能を実現するプログラム）等を読み出して記憶するＲＡＭ（Random Access Memory）等の記憶装置である。このメモリ装置３１は、プログラムの実行に必要なファイル、データ等も記憶する。

演算処理装置３２は、メモリ装置３１に格納されたプログラムを実行するＣＰＵ（Central Processing Unit）等の演算処理装置である。インターフェース装置３３は、サテライトコンパス２やディスプレイ４等の外部機器に接続するためのインターフェース装置である。入力装置３４は、ユーザインターフェースを提供する入力装置（例えばキーボード、マウス）である。

補助記憶装置３５は、プログラムやファイル、データ等を記憶するＨＤＤ（Hard Disk Drive）等の記憶装置である。この補助記憶装置３５には、図２の船体状態データ計算部２３や海象推定部２４の機能を実現するプログラム等が記憶される。

ドライブ装置３６は、記憶媒体３７に記憶されたプログラム（例えば、図２の船体状態データ計算部２３や海象推定部２４の機能を実現するプログラム）を読み出す装置である。ドライブ装置３６によって読み出されたプログラムは、補助記憶装置３５にインストールされる。記憶媒体３７は、上記のプログラム等を記録したＵＳＢ（Universal Serial Bus）メモリ、ＳＤメモリカード等の記憶媒体である。

ディスプレイ４は、情報処理装置３によって生成された出力データ、例えば海象条件を画面出力する出力装置である。

［システムの機能構成］
図２は、本実施形態に係る横メタセンタ高さ推定装置を含む海象推定システムの機能構成例を示す図である。なお、以下の説明においては、図１と同様の構成要素については同一の符号を付して重複する説明を適宜省略する。

図２に示す海象推定システム１は、計測部２１、履歴記憶部２２、船体状態データ計算部（横メタセンタ高さ推定部）２３、海象推定部２４、出力部２５を有する。この海象推定システム１は船舶に搭載される。なお、履歴記憶部２２、船体状態データ計算部２３及び海象推定部２４の各機能は、情報処理装置３によって実現される。

計測部２１は、当該海象推定システム１が搭載された船舶の船体運動データを計測する計測手段である。ここでいう船体運動データとは、船体の横揺れ角、縦揺れ角及び上下揺れの変位等の船体の運動に係るデータである。なお、船体の横揺れ、縦揺れの角速度及び上下揺れの加速度でも良い。この計測部２１は図１のサテライトコンパス２やジャイロセンサによって実現される。

履歴記憶部２２は、計測部２１によって計測された船体運動データの履歴を記憶する履歴記憶手段である。この履歴記憶部２２は、過去から現在までの一定期間の船体運動データの時系列データを記憶する。この履歴記憶部２２は図１のメモリ３１等によって実現される。なお、船体運動データの履歴の入力元は計測部２１に限定されるものではない。例えば、船体運動データの履歴を記憶している他の情報処理装置を入力元としても良い。

船体状態データ計算部２３は、履歴記憶部２２によって記憶された一定期間の船体運動データの履歴に基づいて船体状態データを計算する船体状態データ計算手段である。ここでいう船体状態データとは、船体の喫水、排水量、ＧＭ等の船体の状態に係るデータである。この船体状態データ計算部２３は図１の演算処理装置３２等によって実現される。この船体状態データ計算部２３が、特許請求の範囲の「横メタセンタ高さ推定手段」に対応する。

海象推定部２４は、履歴記憶部２２によって記憶された船体運動データの履歴と船体状態データ計算部２３によって計算された船体状態データとに基づいて、当該海象推定システム１が搭載された船舶が航行中の海域における局所的な海象条件を推定する海象推定手段である。ここでいう海象条件とは、船舶の航行領域における波浪の波高、波周期、波向き等の海象に係る情報である。この海象推定部２４は図１の演算処理装置３２等によって実現される。

出力部２５は、船体状態データ計算部２３によって計算された船体状態データや、海象推定部２４によって推定された海象条件を出力する出力手段である。この出力部２５は図１のディスプレイ４等によって実現される。

以上に示す構成により、本実施形態に係る海象推定システム１では、計測部２１によって計測された船体運動データに基づいて、情報処理装置３が有する船体状態データ計算部（横メタセンタ高さ推定部）２３が横メタセンタ高さを推定し、その後海象推定部２４が海象条件を推定する。出力部４は、推定された横メタセンタ高さや推定された海象条件を出力する。

［システムの制御ロジック］
図３は、本実施形態に係る横メタセンタ高さ推定装置を含む海象推定システムの制御ロジックを示すフローチャートである

+
海象推定システム１は図３に示す一連のステップＳ１〜Ｓ８の制御ロジックを繰り返し行うことによって、海象条件を繰り返し推定する。特に、ステップＳ１〜Ｓ３の処理を繰り返し行うことによって、横メタセンタ高さを推定する。なお、以下、適宜図２を参照しながら説明する。

まずステップＳ１において、計測部２１が船体運動データを計測する（Ｓ１）。具体的には、船体の横揺れ角、縦揺れ角、上下揺れの変位のデータを計測する。このステップＳ１は、一連のステップＳ１〜Ｓ８の処理の繰り返しの中で逐次的に行っても良いし、バッチ処理等によってステップＳ２〜Ｓ８とは別個の処理として繰り返し行っても良い。

なお、図３に示すステップＳ１の処理を繰り返し行うことにより、過去から現在までの一定期間の横揺れ角、縦揺れ角、上下揺れの変位の時系列データ、すなわち時間毎の角度等のデータが、履歴記憶部２２に記憶される。その後、ステップＳ２、Ｓ４及びＳ６の各処理に進む。

次にステップＳ２において、船体状態データ計算部２３は、横揺れ角の時系列データに基づき、横揺れ固有周波数を計算（推定）する（Ｓ２）。ステップＳ２に係る処理は既知の技術であるが、その一例を以下に説明する。

すなわち、まず横揺れの時系列データx(t)に関する２階の線形確率力学モデル（次式（２）参照）を考える。なお、式（２）においてa₁(=2α)は減衰係数、a₂(=ω²)は固有角周波数ωの２乗、u(t)は確率過程として取り扱われる外力項を表しており、有限な分散を有するが特に白色性は仮定しない。

また、式（２）における外力項u(t)は、次式（３）に示すｍ階の連続型自己回帰モデルで表されるものとする。なお、式（３）においてb_i(i=1,・・・,m)はモデルの係数、v(t)は平均0、分散σ^２の正規白色雑音である。

式（２）を式（３）に代入することによって、次式（４）に示す白色化された(m+2)階の連続型自己回帰モデルが得られる。なお、式（４）においてc_i(i=1,・・・,m+2)はモデルの係数である。

式（４）をベクトル形式で表示すると、次式（５）のように表される。

この式（５）を離散化した上で、状態空間モデルのシステムモデルとして取り扱う。さらに、未知係数c_i(i=1,・・・,m+2)を同時に推定するために、状態空間モデルの状態ベクトルに前述の未知係数a_i，b_iを含めて考えることで、モデルを自己組織型状態空間モデルに拡張し、アンサンブルカルマンフィルタを用いた状態推定及び未知係数の推定を同時に行う。アンサンブルカルマンフィルタを用いた状態推定については既知の技術であるため、ここでは説明を省略する。

以上に示す手順により、ステップＳ２では、履歴記憶部２２に記憶された横揺れ角の時系列データに基づいて、横揺れ固有周波数を計算（推定）する。なお、上記した方法以外の方法によって横揺れ固有周波数を計算しても良い。例えば、離散型の自己回帰モデルを用いて横揺れ固有周波数を計算しても良い。

その後ステップＳ３において、船体状態データ計算部２３は、ステップＳ２で算出された横揺れ固有周波数に基づき、ＧＭを計算（推定）する（Ｓ３）。

ステップＳ３では、ステップＳ２で算出された横揺れ固有周波数（又は逆数である横揺れ固有周期）を観測データとし、ＧＭと慣動半径を状態変数とする非線形観測モデルにおいて、状態変数が時間の経過と共に若干揺らぐと仮定してこれをシステムモデルと考え、一般状態空間モデル解析を行うことによってＧＭと慣動半径とを同時に推定する。

すなわち、ステップＳ２で算出された横揺れ固有周波数f(=ω/2π)とＧＭとの間には、次式（６）に示す関係が成立する。なお、式（６）においてTは横揺れ固有周期（単位：ｓ）、ｆは横揺れ固有周波数（単位：Ｈｚ）、kは慣動半径（単位：ｍ）、gは地球の重力加速度（単位：ｍ／ｓ^２）である。式（６）においてＧＭ（単位：ｍ）及びkが未知数である。

そこで、ステップＳ３では、二つの未知数ＧＭ及び慣動半径ｋを状態変数とする状態空間モデルを考える。すなわち、時刻ｎにおいてＧＭ_ｎ及びｋ_ｎの状態推定量から横揺れ固有周期Ｔ_ｎが観測される非線形観測モデルを考える。状態変数が時間の経過と共に若干揺らぐと仮定してこれをシステムモデルと考えることによって、次式（７）に示す一般状態空間モデルを構成する。なお、式（７）においてv_nはシステムノイズであり、w_nは観測ノイズである。いずれも簡単のため正規白色雑音とする。

ステップＳ３では、式（７）で表される一般状態空間モデルに基づいて、粒子フィルタの一種であるモンテカルロフィルタを用いて状態推定を行う。なお、モンテカルロフィルタを用いた状態推定については既知の技術であるため、ここでは説明を省略する。

以上に示すステップＳ２及びＳ３の処理により、船体状態データ計算部２３は、履歴記憶部２２に記憶された一定期間の横揺れ角の時系列データを解析することによって、船体のＧＭをリアルタイムに計算（推定）している。この手法によれば、載荷状態での重心位置の変化や船の運動程度をリアルタイムに把握することができる。

また、ＧＭを推定するために一切の近似を利用しないため、概ね安定したＧＭの推定結果を得ることができる。例えば船舶の設計当時のＧＭが０．５２である場合に、式（１）で示す従来手法ではＧＭの算出値が０．６８程度であったのに対し、本実施形態に係る手法では０．４８〜０．５４の間の推定値が算出された。従って、本実施形態に係る手法によれば、精度良くＧＭを推定することが可能である。なお、上記式（７）においては、一期前よりも過去の変数（例えばX_n-2）が含まれる式にしても良い。

図３に戻り、ステップＳ１からステップＳ４に進んだ場合、船体状態データ計算部２３は、ステップＳ１で履歴記憶部２２に記憶された上下揺れの変位の時系列データを解析することによって、船体の喫水、排水量を計算する（Ｓ４）。例えば、上下揺れの変位の時系列データ、ＧＰＳアンテナ（図１の計測部２１に対応）の設置高さ及び船舶の長手方向の傾斜角等に基づいて喫水を計算し、さらに計算された喫水に基づいて排水量を計算する。

ステップＳ３及びＳ４からステップＳ５に進んだ場合、海象推定部２４は、ステップＳ２で計算された船体状態データに基づいて、当該船体の現在の船体応答を計算する（Ｓ５）。

ステップＳ５の処理に際しては、以下の２つの方法があり、いずれかが採用される。一つ目の方法は、予め船体状態データ（喫水、ＧＭ）と船速と入力となる海象（波高、波周期、波向き）とをパラメータとして船体応答関数を計算してデータベースを作成しておき、現在のこれらの状態に相当する船体応答関数を補間計算により求める方法である。二つ目の方法は、現在の船体の状態、海象を入力として計算式に基づいて船体応答関数を計算する方法である。いずれにしても、海象は今から求める未知数であるが、船体応答関数を計算するには必要な項であり、非線形問題として繰り返し法により最適な応答関数をリアルタイムに選択する。なお、ここでいう応答関数とは、規則的な波長の波を任意の方向から受けた場合に船体がどのように応答（運動）するかを示す波向及び波長等をパラメータとする関数である。

このステップＳ５に係る処理では、計測部２１（サテライトコンパス２）によって計測される現在の船体運動（ピッチ、ロール、ヒーブ）と船速に応じて最適な船体応答関数をリアルタイムに選択している。これにより、船速が刻々と変化する実際の海上であっても、最適な船体の応答関数を決定することができ、その結果として、後述する海象条件の推定に係る精度も向上させることができる。

また、ステップＳ１からステップＳ６に進んだ場合、船体状態データ計算部２３は、ステップＳ１で履歴記憶部２２に記録された上下揺れの変位、縦揺れ角、横揺れ角に基づいて、各船体運動（横揺れ、縦揺れ及び上下揺れ）のクロススペクトルを計算する（Ｓ６）。ステップＳ６に係る処理は既知の技術であるため、ここでは説明を省略する。

ステップＳ６において船体状態データ計算部２３は、上下揺れの変位（単位：m）、縦揺れ角（単位：rad）、横揺れ角（単位：rad）の時系列データに基づいて、横揺れオートスペクトル（単位：rad²/s）、縦揺れオートスペクトル（単位：rad²/s）上下揺れオートスペクトル（単位：m²/s）、縦揺れ−上下揺れのクロススペクトル（単位：rad・m／s）、横揺れ−上下揺れのクロススペクトル（単位：rad・m／s）、縦揺れ−横揺れのクロススペクトル（単位：rad²／s）からなる各船体運動（横揺れ、縦揺れ及び上下揺れ）のクロススペクトルを、周波数毎に求める。周波数毎に求められたクロススペクトルは、時系列データとして履歴記憶部２２に記憶される。

ステップＳ５及びＳ６からステップＳ７に進んだ場合、海象推定部２４は、ステップＳ５で計算された現在の船体応答と、ステップＳ６で計算された各船体運動（横揺れ、縦揺れ及び上下揺れ）のクロススペクトルとに基づいて、方向波スペクトルを確率統計的に計算する（Ｓ７）。以下、ステップＳ７で用いられる本実施形態に係る理論とともに、当該処理を説明する。

海洋波があらゆる方向から到来する全ての周波数を含む成分波の重ね合わせで表現できるものとすると、ある時間tにおける固定点（船舶位置）での海面変動量η(t)は、方向波スペクトルE(f,x) （単位：m²／（rad／s））を用いて次式（８）で示される。なお、式（８）においてルート記号で囲まれた部分及びε(f,x)のそれぞれは、周波数fで方向xから到来する成分波の振幅、位相である。

一方、船体動揺が波浪入力に対して線形応答であると仮定すると、ある波の出会い周波数f_eにおける方向波スペクトルE(f_e,x)と船体動揺のクロススペクトルφ_ln(f_e)との関係は、一般的に次式（９）で示される。なお、式（９）においてlとnは船体動揺のモードであり、H_l(f_e,x)、H_n*( f_e,x)は各々の動揺モードの応答関数である。また、xは波との出会い角、記号(*)は複素共役である。

この式（９）は出会い周波数をベースに示されているので、これを絶対周波数をベースとする式（次式（１０）参照）に変換する。

式（１０）において、右辺第２項から第４項は追波時の寄与、すなわちクロススペクトルに含まれる追い波中を航行する際の波の周波数成分の度合いを示している。パラメータＡ、絶対周波数と対応する３つの出会い周波数f₀₁、f₀₂、f₀₃及びヤコビアンはそれぞれ次式（１１）に示すように定義される。なお、式（１１）においてUは船速、gは重力加速度である。

ここで、出会い角xに関する積分範囲を十分に大きな数K個の微小区間に分けた場合の微小積分区間内での変動量の応答関数及び方向波スペクトルは、一定とみなすことが可能である。そのため、式（１０）を次式（１２）のように離散化することができる。なお、式（１２）においてK1(0≦K1≦K/2)は、離散的積分範囲の中で追波状態になるものの個数を表している。

ここで、縦揺れ角、横揺れ角及び上下変位の各々を任意の変動量θ、φ、ηとした場合のクロススペクトルΦ(f_e)は３×３行列となり、式（１２）は、次式（１３）のようにマトリックス表示できる。なお、式（１３）においてH(f₀₁)は3×K行列、H(f₀₂)及びH(f₀₃)は3×K1行列、E(f₀₁)はK×Kの対角行列、E(f₀₂)及びE(f₀₃)はK1×K1の対角行列である。また、記号（T）は転置行列である。

クロススペクトル行列Φ(f_e)はエルミート行列であるから、上三角行列のみを取り扱えばよい。また、式（１３）において実数部と虚数部を分けるとともに、観測に伴う誤差項Wを導入して表記する場合、式（１３）は次式（１４）で示す線形回帰モデルで表すことができる。

式（１４）において、yはクロススペクトル行列Φ(f_e)の実部と虚部で構成されるベクトルである。Aは船体動揺の応答関数の理論値で構成される係数行列である。Wは統計的な性質が平均0、分散共分散行列Σにしたがう白色雑音である。xは離散化された方向波スペクトルから構成される未知ベクトルである。

この式（１４）において、クロススペクトルが時系列的に得られるものとすると、それに対応して、方向波スペクトルを時系列的に推定することが可能になる。これは、式（１４）を時変システムとして捉えることに相当し、時刻を添字tで表して次式（１５）のように拡張できる。

式（１５）は、一般状態空間モデルにおける観測モデルと形式的に等価である。従って、方向波スペクトルが時間に関して滑らかに変化するという平滑化事前分布をシステムモデルとして導入する（次式（１６）参照）ことによって、方向波スペクトルの推定を次式（１６）に示す一般状態空間モデルの状態推定の問題に帰着できる。

式（１６）において、x_tは状態ベクトル、v_tはシステムノイズベクトル、y_tは観測ベクトル、A_tは状態遷移行列、W_tは観測ノイズベクトルである。ここで、方向波スペクトルが非負であることを考慮して状態ベクトルx_tの対数を改めてx_tと置き換えると、式（１６）は、次式（１７）に示す一般状態空間モデルに変形される。

ここで、F(x_t)は、全ての要素に対して指数をとることを意味している。また、状態ベクトルの要素は、次式（１８）のように構成される。

式（１８）において、ｍは波の絶対周波数の分割数である。式（１７）は、非線形観測モデル、すなわち非線形な状態空間モデルである。従って、状態推定には非線形フィルタリングに有効な方法を用いる必要がある。従前は粒子フィルタを用いていたが、この方法は計算負荷が非常に高い。そこで、本実施形態ではアンサンブルカルマンフィルタによる状態推定法を導入するが、非線形観測モデルである式（１７）そのままの形ではアンサンブルカルマンフィルタを適用できない。この問題を解決するために、次式（１９）に示す拡大された状態ベクトルを考える。

また、次式（２０）に示す拡大観測行列及び拡大状態遷移ベクトルを考える。

その結果、x_tについては次式（２１）となり、拡大システムモデルが得られる。またy_tについても次式（２２）に示す形式的に線形な拡大観測モデルを得ることができる。x_t及びy_tは線形観測の拡大状態空間モデルとなるので、アンサンブルカルマンフィルタによる状態推定を実現できる。なお、アンサンブルカルマンフィルタの適用については既知の技術であるため、ここでは説明を省略する。

以上に示すステップＳ７の処理により、海象推定部２４は、ステップＳ５で計算された船体応答と、ステップＳ６で計算された各船体運動のクロススペクトルとに基づいて、方向波スペクトルを確率統計的に計算する（Ｓ７）。

このステップＳ７の処理によれば、過去から現在にかけての一定期間の船体応答及び各船体運動のクロススペクトルの時系列データを確率統計的に処理することによって、リアルタイムに方向波スペクトルを計算している。そのため、精度の高い方向波スペクトルを導き出すことができる。

またこの手法によれば、アンサンブルカルマンフィルタによる状態推定に基づいて方向波スペクトルを推定しているので、従前のモンテカルロフィルタによる方法と比較して、格段に短い計算時間で精度の高い方向波スペクトルの推定が実現できる。

なお、ステップＳ７の処理が終了すると、ステップＳ８に進み、海象推定部２４は、ステップＳ７で計算された方向波スペクトルに基づいて海象条件を推定する（Ｓ８）。ステップＳ８では、ステップＳ７で計算された方向波スペクトルに基づいて、船舶が航行中の局所的な海域における波向き、波周期、有義波高等の海象条件を推定することができる。

以上、本発明の一実施形態について説明したが、上記実施形態は本発明の適用例の一つを示したものであり、本発明の技術的範囲を上記実施形態の具体的構成に限定する趣旨ではない。

１海象推定システム
２サテライトコンパス
３情報処理装置（横メタセンタ高さ推定装置）
４ディスプレイ
２１計測部
２２履歴記憶部
２３船体状態データ計算部（横メタセンタ高さ推定部）
２４海象推定部
２５出力部

Claims

不規則に変動する波浪上を航行する船体の横揺れ角の時系列データを記憶する履歴記憶手段と、
前記履歴記憶手段によって記憶された船体の横揺れ角の時系列データに基づいて、当該船体の横メタセンタ高さを推定する横メタセンタ高さ推定手段と、を有し、
前記横メタセンタ高さ推定手段は、前記船体の横揺れ角の時系列データに基づいて横揺れ固有周波数を計算し、計算された横揺れ固有周波数を観測モデルとし、横メタセンタ高さと当該船体の慣動半径とを状態変数とする一般状態空間モデルに基づく状態推定により、前記横メタセンタ高さを推定することを特徴とする横メタセンタ高さ推定装置。
前記横メタセンタ高さ推定手段は、前記横揺れ固有周波数と前記横メタセンタ高さとの関係式（１）
［数１］
（但し、Ｔは横揺れ固有周期、ｆは横揺れ固有周波数、ｋは慣動半径、ｇは重力加速度、ＧＭは横メタセンタ高さ）において、時刻ｎにおける横メタセンタ高さＧＭ_ｎ及び慣動半径ｋ_ｎを状態変数とし、横揺れ固有周期Ｔ_ｎを観測データとする式（２）
［数２］
（但し、ｖ _ｎはシステムノイズ、ｗ _ｎは観測ノイズ）で表される一般状態空間モデルに基づく状態推定を行うことによって、前記横メタセンタ高さＧＭを推定することを特徴とする請求項１に記載の横メタセンタ高さ推定装置。
不規則に変動する波浪上を航行する船体の横揺れ角の時系列データを履歴記憶部に記憶する履歴記憶工程と、
前記履歴記憶部に記憶された船体の横揺れ角の時系列データに基づいて、当該船体の横メタセンタ高さを推定する横メタセンタ高さ推定工程と、を含み、
前記横メタセンタ高さ推定工程では、まず前記船体の横揺れ角の時系列データに基づいて横揺れ固有周波数を計算し、計算された横揺れ固有周波数を観測モデルとし、横メタセンタ高さと当該船体の慣動半径とを状態変数とする一般状態空間モデルに基づく状態推定により、前記横メタセンタ高さを推定することを特徴とする横メタセンタ高さ推定方法。