特許 6558760 海象推定装置及び海象推定方法

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

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

(11)【特許番号】6558760

(24)【登録日】2019年7月26日

(45)【発行日】2019年8月14日

(54)【発明の名称】海象推定装置及び海象推定方法

(51)【国際特許分類】

   G01W 1/00 20060101AFI20190805BHJP

   G08G 3/00 20060101ALI20190805BHJP

【ＦＩ】

   G01W1/00 Z

   G08G3/00 A

【請求項の数】4

【全頁数】21

(21)【出願番号】特願2017-551970(P2017-551970)

(86)(22)【出願日】2016年11月20日

(86)【国際出願番号】JP2016084373

(87)【国際公開番号】WO2017086482

(87)【国際公開日】20170526

【審査請求日】2018年3月29日

(31)【優先権主張番号】特願2015-228228(P2015-228228)

(32)【優先日】2015年11月20日

(33)【優先権主張国】JP

(73)【特許権者】

【識別番号】502298192

【氏名又は名称】流体テクノ株式会社

(73)【特許権者】

【識別番号】501168814

【氏名又は名称】国立研究開発法人水産研究・教育機構

(73)【特許権者】

【識別番号】000166247

【氏名又は名称】古野電気株式会社

(74)【代理人】

【識別番号】100120086

【弁理士】

【氏名又は名称】▲高▼津 一也

(72)【発明者】

【氏名】玉島 正裕

(72)【発明者】

【氏名】寺田 大介

(72)【発明者】

【氏名】松田 秋彦

(72)【発明者】

【氏名】前野 仁

(72)【発明者】

【氏名】林 孝二

【審査官】 伊藤 昭治

(56)【参考文献】

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

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

【文献】 寺田大介 ほか，逐次型Ｂａｙｅｓ法による方向波スペクトルのオンライン推定，日本航海学会論文集，２００３年，１０８巻，９７−１０４頁

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

Ｇ０１Ｗ １／００

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

【請求項1】

船体の横揺れ角、縦揺れ角及び上下揺れ変位に係る船体運動データの履歴を記憶する履歴記憶手段と、
前記履歴記憶手段によって記憶された前記船体運動データの履歴に基づいて、当該船体の現在の船体状態データを計算する船体状態データ計算手段と、
前記履歴記憶手段によって記憶された前記船体運動データの履歴に基づいて、当該船体の横揺れ、縦揺れ及び上下揺れの各船体運動のクロススペク卜ルを計算するクロススペク卜ル計算手段と、
前記船体状態データ計算手段によって計算された船体状態データに基づいて、船体応答関数を計算する船体応答関数計算手段と、
前記クロススペク卜ル計算手段によって計算された各船体運動のクロススペク卜ルと前記船体応答関数計算手段によって計算された船体応答関数とに基づいて、当該船体が航行中の海域における海象条件を推定する海象推定手段と、
を備えた海象推定装置であって、
前記船体状態データは、当該船体の喫水、排水量及び横メタセンタ高さを含む
ことを特徴とする海象推定装置。

【請求項2】

船体の横揺れ角、縦揺れ角及び上下揺れ変位に係る船体運動データの履歴を記憶する履歴記憶手段と、
前記履歴記憶手段によって記憶された前記船体運動データの履歴に基づいて、当該船体の現在の船体状態データを計算する船体状態データ計算手段と、
前記履歴記憶手段によって記憶された前記船体運動データの履歴に基づいて、当該船体の横揺れ、縦揺れ及び上下揺れの各船体運動のクロススペク卜ルを計算するクロススペク卜ル計算手段と、
前記船体状態データ計算手段によって計算された船体状態データに基づいて、船体応答関数を計算する船体応答関数計算手段と、
前記クロススペク卜ル計算手段によって計算された各船体運動のクロススペク卜ルと前記船体応答関数計算手段によって計算された船体応答関数とに基づいて、当該船体が航行中の海域における海象条件を推定する海象推定手段と、
を備えた海象推定装置であって、
サテライトコンパスと、
当該船体の規則波における応答関数と、前記船体状態データ及び当該船体の船速とを対応付けて記憶する応答関数記憶手段と、をさらに備え、
前記船体応答関数計算手段は、前記現在の船体状態データと前記サテライトコンパスによって計測される当該船体の船速とに対応付けられた応答関数を選択し、選択された前記応答関数に基づいて当該船体の船体応答関数を計算する
ことを特徴とする海象推定装置。

【請求項3】

前記海象推定手段は、前記各船体運動のクロススペク卜ルと前記船体応答関数とに基づいて方向波スペク卜ルを確率統計的に計算し、計算された方向波スペク卜ルに基づいて海象条件を推定することを特徴とする請求項１又は２に記載の海象推定装置。

【請求項4】

船体の横揺れ角、縦揺れ角及び上下揺れ変位に係る船体運動データの履歴を記憶する履歴記憶部を備えた海象推定装置における海象推定方法であって、
船体状態データ計算部が、前記履歴記憶部において記憶された前記船体運動データの履歴に基づいて、当該船体の現在の船体状態データを計算する船体状態データ計算工程と、
クロススペク卜ル計算部が、前記履歴記憶部において記憶された前記船体運動データの履歴に基づいて、当該船体の横揺れ、縦揺れ及び上下揺れの各船体運動のクロススペク卜ルを計算するクロススペク卜ル計算工程と、
船体応答関数計算部が、前記船体状態データ計算工程において計算された船体状態データに基づいて、船体応答関数を計算する船体応答関数計算工程と、
海象推定部が、前記クロススペクトル計算工程において計算された各船体運動のクロススペク卜ルと前記船体応答関数計算工程において計算された船体応答関数とに基づいて、当該船体が航行中の海域における海象条件を推定する海象推定工程と、
を含み、
前記船体状態データは、当該船体の喫水、排水量及び横メタセンタ高さを含む
ことを特徴とする海象推定方法。

【発明の詳細な説明】

【参照による取り込み】

【0001】

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

【技術分野】

【0002】

本発明は、海象推定装置及び海象推定方法に関する。

【背景技術】

【0003】

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

【0004】

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

【0005】

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

【0006】

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

【0007】

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

【先行技術文献】

【非特許文献】

【0008】

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

【発明の概要】

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

【0009】

ところで、上記非特許文献１に開示された技術では、海象条件を推定する際に船体応答関数を使用する。船体応答関数とは、規則的な波長の波を任意の方向から受けた場合に船体がどのように応答（運動）するかを示す波向及び波長等をパラメータとする関数である。非特許文献１に開示された技術では、この船体応答関数は、固定値の船速を基に設定されていた。

【0010】

しかしながら、実際の海上では船速は刻々と変化するものであるため、固定値の船速を基に決定される上記船体応答関数の精度は低かった。またその結果として、海象条件の推定に係る精度も低いものとなっていた。

【0011】

本発明は、上記のような事情に鑑みてなされたものであって、船速の変化に応じて精度良く海象条件の推定が可能な海象推定装置及び海象推定方法を提供することを目的とする。

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

【0012】

上記の目的を達成するために、本発明に係る海象推定装置は、船体の横揺れ角、縦揺れ角及び上下揺れ変位に係る船体運動データの履歴を記憶する履歴記憶手段と、前記履歴記憶手段によって記憶された前記船体運動データの履歴に基づいて、当該船体の現在の船体状態データを計算する船体状態データ計算手段と、前記履歴記憶手段によって記憶された前記船体運動データの履歴に基づいて、当該船体の横揺れ、縦揺れ及び上下揺れの各船体運動のクロススペクトルを計算するクロススペクトル計算手段と、前記船体状態データ計算手段によって計算された船体状態データに基づいて、船体応答関数を計算する船体応答関数計算手段と、前記クロススペクトル計算手段によって計算された各船体運動のクロススペクトルと前記船体応答関数計算手段によって計算された船体応答関数とに基づいて、当該船体が航行中の海域における海象条件を推定する海象推定手段と、を備えたことを特徴とする。

【0013】

また上記の目的を達成するために、本発明に係る海象推定方法は、船体の横揺れ角、縦揺れ角及び上下揺れ変位に係る船体運動データの履歴を記憶する履歴記憶部を備えた海象推定装置における海象推定方法であって、船体状態データ計算部が、前記履歴記憶部において記憶された前記船体運動データの履歴に基づいて、当該船体の現在の船体状態データを計算する船体状態データ計算工程と、クロススペクトル計算部が、前記履歴記憶部において記憶された前記船体運動データの履歴に基づいて、当該船体の横揺れ、縦揺れ及び上下揺れの各船体運動のクロススペクトルを計算するクロススペクトル計算工程と、船体応答関数計算部が、前記船体状態データ計算工程において計算された船体状態データに基づいて、船体応答関数を計算する船体応答関数計算工程と、海象推定部が、前記クロススペクトル計算工程において計算された各船体運動のクロススペクトルと前記船体応答関数計算工程において計算された船体応答関数とに基づいて、当該船体が航行中の海域における海象条件を推定する海象推定工程と、を含むことを特徴とする。

【発明の効果】

【0014】

本発明によれば、船速の変化に応じて精度良く海象条件を推定することができる。

【図面の簡単な説明】

【0015】

【図1】本実施形態に係る海象推定装置を含む海象推定システムのハードウェア構成例を示す図である。

【図2】本実施形態に係る海象推定装置を含む海象推定システムの機能構成例を示す図である。

【図3】本実施形態に係る海象推定装置を含む海象推定システムの制御ロジックの一例を示すフローチャートである。

【図4】本実施形態に係る船体応答関数計算部が使用するデータベースのデータ構造の一例を示す図である。

【図5A】本実施形態に係る船体応答関数を説明するための図である。

【図5B】本実施形態に係る船体応答関数を説明するための図である。

【図5C】本実施形態に係る船体応答関数を説明するための図である。

【図6A】本実施形態に係る船体運動データの一例を示す図である。

【図6B】本実施形態に係る船体運動データの一例を示す図である。

【図6C】本実施形態に係る船体運動データの一例を示す図である。

【図7A】本実施形態に係るクロススペクトルの一例を示す図である。

【図7B】本実施形態に係るクロススペクトルの一例を示す図である。

【図7C】本実施形態に係るクロススペクトルの一例を示す図である。

【図7D】本実施形態に係るクロススペクトルの一例を示す図である。

【図7E】本実施形態に係るクロススペクトルの一例を示す図である。

【図7F】本実施形態に係るクロススペクトルの一例を示す図である。

【図8A】本実施形態に係る方向波スペクトルの一例を示す図である。

【図8B】本実施形態に係る方向波スペクトルの一例を示す図である。

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

【0016】

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

【0017】

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

【0018】

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

【0019】

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

【0020】

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

【0021】

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

【0022】

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

【0023】

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

【0024】

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

【0025】

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

【0026】

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

【0027】

図２に示す海象推定システム１は、計測部２１、履歴記憶部２２、クロススペクトル計算部２３、船体状態データ計算部２４、船体応答関数計算部２５、海象推定部２６、応答関数記憶部２７、出力部２８を有する。この海象推定システム１は船舶に搭載される。なお、履歴記憶部２２、クロススペクトル計算部２３、船体状態データ計算部２４、船体応答関数計算部２５、海象推定部２６及び応答関数記憶部２７の各機能は、情報処理装置３によって実現される。

【0028】

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

【0029】

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

【0030】

クロススペクトル計算部２３は、履歴記憶部２２によって記憶された一定期間の船体運動データの時系列データに基づいて、各船体運動（横揺れ、縦揺れ及び上下揺れ）のクロススペクトルを計算する。このクロススペクトル計算部２３は図１の演算処理装置３２等によって実現される。

【0031】

船体状態データ計算部２４は、履歴記憶部２２によって記憶された一定期間の船体運動データの履歴である時系列データに基づいて船体状態データを計算する船体状態データ計算手段である。ここでいう船体状態データとは、船体の喫水、排水量、ＧＭ等の船体の状態に係るデータである。この船体状態データ計算部２４は図１の演算処理装置３２等によって実現される。

【0032】

船体応答関数計算部２５は、船体状態データ計算部２４によって計算された船体状態データと、応答関数記憶部２７に記憶された船体応答関数に係るデータとに基づいて、当該船体の船体応答関数を計算する船体応答関数計算手段である。この船体応答関数計算部２５は図１の演算処理装置３２等によって実現される。

【0033】

海象推定部２６は、クロススペクトル計算部２３によって計算されたクロススペクトルと船体応答関数計算部２５によって計算された船体応答関数とに基づいて、当該海象推定システム１が搭載された船舶が航行中の海域における局所的な海象条件を推定する海象推定手段である。ここでいう海象条件とは、船舶の航行領域における波浪の波高、波周期、波向き等の海象に係る情報である。

【0034】

この海象推定部２６は、まずクロススペクトルと船体応答関数とに基づいて、方向波スペクトルを確率統計的に計算する。その後、計算された方向波スペクトルに基づいて海象条件を推定する。詳細処理については後述する。この海象推定部２６は図１の演算処理装置３２等によって実現される。

【0035】

応答関数記憶部２７は、船体応答関数計算部２５が船体応答関数を計算する際に用いる船体応答関数に係るデータを記憶する応答関数記憶手段である。船体応答関数とは、規則的な波長の波を任意の方向から受けた場合に船体がどのように応答（運動）するかを示す波向及び波長等をパラメータとする関数である。また、船体応答関数に係るデータとは、船体の規則波における船体応答関数と、船体状態データ（喫水やＧＭ）、船速や減衰係数とを対応付けたデータベースである。この応答関数記憶部２７は図１のメモリ装置３１等によって実現される。

【0036】

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

【0037】

以上に示す構成により、本実施形態に係る海象推定システム１では、計測部２１によって計測された船体運動データに基づいて、海象推定部２６が海象条件を推定する。出力部２８は、推定された海象条件を出力する。

【0038】

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

【0039】

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

【0040】

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

【0041】

なお、図３に示すステップＳ１の処理を繰り返し行うことにより、過去から現在までの一定期間の横揺れ角、縦揺れ角、上下揺れの変位の時系列データ、すなわち時間毎の角度等のデータが、履歴記憶部２２に記憶される。例えば、１秒〜１０秒間隔で船体運動データは計測並びに記憶され、各データは小数点以下２桁以上の値で表される。なお、船体運動データは、海象状況による運動以外に針路変更や増減速の影響を含んだ非定常な時系列データを含む。ステップＳ１では、このようなデータを分離し以降のステップで使用しない。この処理には既往の技術である赤池情報量基準（ＡＩＣ：Akaike’s Information Criterion）法が用いられる。その後、ステップＳ２、Ｓ４及びＳ６の各処理に進む。

【0042】

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

【0043】

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

【数1】

【0044】

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

【数2】

【0045】

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

【数3】

【0046】

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

【数4】

【0047】

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

【0048】

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

【0049】

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

【0050】

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

【0051】

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

【数5】

【0052】

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

【数6】

【0053】

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

【0054】

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

【0055】

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

【0056】

図３に戻り、ステップＳ１からステップＳ４に進んだ場合、船体状態データ計算部２４は、ステップＳ１で履歴記憶部２２に記憶された上下揺れの変位の時系列データを解析することによって、現在の船体の喫水、排水量を計算する（Ｓ４）。

【0057】

ステップＳ４では、船体状態データ計算部２４は、例えばＧＰＳアンテナ（図１の計測部２１に対応）の設置高さ、船舶の長手方向の傾斜角、及び海上保安庁が公表している日本沿岸潮汐調和定数表に基づいて計算される潮汐データを基に、喫水を計算（推定）する。また、上下揺れの変位の時系列データに基づいて作成される上下揺れの運動方程式を、ステップＳ１の横揺れの場合と同様に離散化し、上下揺れの固有周波数を計算（推定）する。その後、推定された固有周波数、船体線図に基づいて予め計算される喫水と船舶の長手方向の傾斜角の組合せに対応したＴＰＣ（Tons Per Centimeter）及び排水量の関係を利用して、排水量を計算（推定）する。

【0058】

また例えばステップＳ４では、履歴記憶部２２に記憶された縦揺れ角、上下揺れの変位の時系列データから統計解析手法を用いて排水量とトリムを計算する。これにより、船体のハイドロ計算表から船首尾喫水を確定させる。本実施形態に係る海象推定システム１が搭載された船体が喫水計を備える場合には、喫水計で計測されたデータから排水量を計算する。喫水計がない場合、ステップＳ１〜Ｓ２の場合と同様、自己組織型状態空間モデル解析法を用いて排水量を計算する。

【0059】

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

【0060】

ステップＳ５の処理に際しては、以下の２つの方法があり、いずれかが採用される。一つ目の方法は、予め船体状態データ（喫水、ＧＭ）と船速と入力となる海象（波高、波周期、波向き）とをパラメータとして船体応答関数を計算してデータベース（例えば図４）を作成しておき、現在のこれらの状態に相当する船体応答関数を、1次補間、2次補間等の補間計算により求める方法である。二つ目の方法は、現在の船体の状態、海象を入力として計算式に基づいて船体応答関数を直接計算する方法である。いずれにしても、海象は今から求める未知数であるが、船体応答関数を計算するには必要な項であり、非線形問題として繰り返し法により最適な応答関数をリアルタイムに選択する。なお、ここでいう船体応答関数とは、図５に示すように、任意の波高において船体の各方向（0度：船尾方向〜１８０度：船首方向）から来る波による船体応答を計算したもの、すなわち船体応答の特性を関数表示したものであり、船体状態により異なる。また船体応答とは、波浪などの外力によって船体の受ける影響を表し、船体運動、構造変化（構造応答）として表れる。

【0061】

図４は、本実施形態に係る船体応答関数計算部が使用するデータベースのデータ構造の一例を示す図である。

【0062】

図４に示すデータベース４１は、ＧＭ（１≦ｋ≦ｏ）、排水量（１≦ｊ≦ｎ）、船速（１≦ｇ≦ｍ）毎の船体応答関数ｆ_ｘｙｚ（１≦ｘ≦ｏ、１≦ｙ≦ｎ、１≦ｚ≦ｍ）が格納されている。なお、船体応答関数ｆ_ｘｙｚの添字ｈ、ｐ、ｒはそれぞれ上下揺れ、縦揺れ、横揺れを表す。このようなデータベースが応答関数記憶部２７（図２参照）に記憶される。ステップＳ５において船体応答関数計算部２５は、ステップＳ３及びＳ４で計算された船体状態データに対応する船体応答関数を、データベースから選択する又はデータベースから選択して１次補間、２次補間等の補間計算により求める。

【0063】

図５Ａ〜図５Ｃは、本実施形態に係る船体応答関数を説明するための図である。

【0064】

図５Ａの縦軸は上下揺れの運動振幅（Ａｍｐ．）／（波高／２）（Ａｗ）であり、横軸は波長／垂線間長さ（Ｌｐｐ）である。図５Ｂの縦軸は縦揺れの運動振幅（Ａｍｐ．）／（波高／２）（Ａｗ）であり、横軸は波長／垂線間長さ（Ｌｐｐ）である。図５Ｃの縦軸は横揺れの運動振幅（Ａｍｐ．）／（波高／２）（Ａｗ）であり、横軸は波長／垂線間長さ（ＬＰＰ）である。図５Ａ〜図５Ｃでは、船体の上下揺れ、縦揺れ、横揺れそれぞれの応答関数ｆ_ｘｙｚについて、海象条件、波高、波長（λ）、波向き（ω）をパラメータとし運動振幅の計算例を示している。

【0065】

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

【0066】

また、ステップＳ１からステップＳ６に進んだ場合、クロススペクトル計算部２３が、ステップＳ１で履歴記憶部２２に記録された上下揺れの変位、縦揺れ角、横揺れ角の時系列データに基づいて、各船体運動（横揺れ、縦揺れ及び上下揺れ）のクロススペクトルを計算する（Ｓ６）。

【0067】

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

【0068】

なお、ステップＳ６においてクロススペクトル計算部２３は、自己回帰モデル解析方を用いてクロススペクトルを計算する。すなわち、上下揺れの変位、縦揺れ角、横揺れ角の時系列データに自己回帰モデルを当てはめ、推定された自己回帰係数と分散共分散行列を用いてクロススペクトルを計算する。この際、モデル次数の決定が重要であるが、ＡＩＣ法を用いて自動的にモデル次数を決定する。これら一連の手順は、ＭＡＩＣＥ（Minimum AIC Estimation）法と呼ばれ、１９８０年代に確立された方法である。

【0069】

ステップＳ５及びＳ６からステップＳ７に進んだ場合、海象推定部２６は、ステップＳ５で計算された現在の船体応答関数と、ステップＳ６で計算された各船体運動（横揺れ、縦揺れ及び上下揺れ）のクロススペクトルとに基づいて、方向波スペクトルを確率統計的に計算する（Ｓ７）。

【0070】

図６Ａ〜図６Ｃは、本実施形態に係る船体運動データの一例を示す図である。

【0071】

図６Ａ〜図６Ｃでは、長さ７０ｍ程度の船についてステップＳ１で計測される波浪中の船体運動データを示している。すなわち、図６Ａ〜図６Ｃは、それぞれ縦揺れ角、横揺れ角、上下揺れの変位の時系列データを示している。

【0072】

図７Ａ〜図７Ｆは、本実施形態に係るクロススペクトルの一例を示す図である。

【0073】

図７Ａ〜図７Ｆでは、ステップＳ６において計算されるクロススペクトルの解析結果の一例を示している。すなわち、図７Ａ〜図７Ｆでは、それぞれ縦揺れオートスペクトル、縦揺れ−横揺れクロススペクトル、縦揺れ−上下揺れのクロススペクトル、横揺れオートスペクトル、横揺れ−上下揺れのクロススペクトル、上下揺れオートスペクトルを示している。

【0074】

図８Ａ及び図８Ｂは、本実施形態に係る方向波スペクトルの一例を示す図である。

【0075】

図８Ａ及び図８Ｂでは、ステップＳ７で逆計算により得られた方向波スペクトルの1次元スペクトル、方向分布をそれぞれ示している。この1次元スペクトル、方向分布に基づいて、海象条件、波高、波周期及び波向きが求められる。以下、ステップＳ７で用いられる本実施形態に係る理論とともに、当該処理を説明する。

【0076】

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

【数7】

【0077】

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

【数8】

【0078】

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

【数9】

【0079】

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

【数10】

【0080】

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

【数11】

【0081】

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

【数12】

【0082】

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

【数13】

【0083】

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

【0084】

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

【数14】

【0085】

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

【数15】

【0086】

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

【数16】

【0087】

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

【数17】

【0088】

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

【数18】

【0089】

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

【数19】

【0090】

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

【数20】

【数21】

【0091】

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

【0092】

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

【0093】

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

【0094】

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

【0095】

具体的には、方向波スペクトルを周波数に関して積分することにより、波向きを推定し、方向波スペクトルを方向角に関して積分することにより、波周期（平均、ゼロクロス、極値間）、有義波高を推定する。なお、安定的に正常な解析結果を得るために、正常な解析結果と異常な解析結果のスペクトルのパターンを認識し、異常な解析結果をスクリーニングする方法を用いる。前述の通り、ステップＳ６に係るクロススペクトルの計算は、自己回帰モデルを当てはめることによって行う。この際、モデル次数の決定にＡＩＣを用いているため、ＡＩＣの値を比較することによって、時系列の定常性を判定できる（この方法は、局所定常自己回帰モデル解析法と呼ばれる）。異常な解析結果のスクリーニングは、次の手順で行う。すなわち、まず第１の手順で、ＡＩＣによる時系列データの定常性を判定する。定常である場合に次に進んで、第２の手順で、クロススペクトルを計算する。その後第３の手順で、方向波スペクトルを計算する。その後第４の手順で、推定した方向波スペクトルと動揺の応答関数から応答スペクトルを計算する。その後第５の手順で、第４の手順で計算した応答スペクトルと第２の手順で計算されたクロススペクトルのオート成分を比較する。その後第６の手順で、第５の手順で得られた差分を元に、推定結果の信頼性を機械学習に基づいて評価する。

【0096】

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

【符号の説明】

【0097】

１ 海象推定システム
２ サテライトコンパス
３ 情報処理装置（海象推定装置）
４ ディスプレイ
２１ 計測部
２２ 履歴記憶部
２３ クロススペクトル計算部
２４ 船体状態データ計算部
２５ 船体応答関数計算部
２６ 海象推定部
２７ 応答関数記憶部
２８ 出力部

【図1】

【図2】

【図3】

【図4】

【図5A】

【図5B】

【図5C】

【図6A】

【図6B】

【図6C】

【図7A】

【図7B】

【図7C】

【図7D】

【図7E】

【図7F】

【図8A】

【図8B】