海象推定装置及び海象推定方法

参照による取り込み

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

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

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

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

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

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

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

平山 次清、井関 俊夫、石田 茂資、“遭遇波浪のリアルタイム把握法と最新の成果”、日本造船学会、平成15年12月、p.74−96

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

図1に示す海象推定システム1は、サテライトコンパス2、情報処理装置(海象推定装置)3、ディスプレイ4を備える。この海象推定システム1は船舶の船体に搭載される。

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

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

メモリ装置31は、情報処理装置3の起動時に補助記憶装置35に記憶されたプログラム(図2のクロススペクトル計算部23〜海象推定部26の機能を実現するプログラム)等を読み出して記憶するRAM(Random Access Memory)等の記憶装置である。このメモリ装置31は、プログラムの実行に必要なファイル、データ等も記憶する。

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

補助記憶装置35は、プログラムやファイル、データ等を記憶するHDD(Hard Disk Drive)等の記憶装置である。この補助記憶装置35には、図2のクロススペクトル計算部23〜海象推定部26の機能を実現するプログラム等が記憶される。

ドライブ装置36は、記憶媒体37に記憶されたプログラム(例えば、図2のクロススペクトル計算部23〜海象推定部26の機能を実現するプログラム)を読み出す装置である。ドライブ装置36によって読み出されたプログラムは、補助記憶装置35にインストールされる。記憶媒体37は、上記のプログラム等を記録したUSB(Universal Serial Bus)メモリ、SDメモリカード等の記憶媒体である。

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

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

図2に示す海象推定システム1は、計測部21、履歴記憶部22、クロススペクトル計算部23、船体状態データ計算部24、船体応答関数計算部25、海象推定部26、応答関数記憶部27、出力部28を有する。この海象推定システム1は船舶に搭載される。なお、履歴記憶部22、クロススペクトル計算部23、船体状態データ計算部24、船体応答関数計算部25、海象推定部26及び応答関数記憶部27の各機能は、情報処理装置3によって実現される。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式(3)をベクトル形式で表示すると、次式(4)のように表される。

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

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

その後ステップS3において、船体状態データ計算部24は、ステップS2で算出された横揺れ固有周波数に基づき、GMを計算(推定)する(S3)。

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

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

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

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

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

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

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

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

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

ステップS3及びS4からステップS5に進んだ場合、船体応答関数計算部25は、ステップS3及びS4で計算された船体状態データに基づいて、当該船体の現在の船体応答関数を計算する(S5)。

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

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

図4に示すデータベース41は、GM(1≦k≦o)、排水量(1≦j≦n)、船速(1≦g≦m)毎の船体応答関数f_xyz(1≦x≦o、1≦y≦n、1≦z≦m)が格納されている。なお、船体応答関数f_xyzの添字h、p、rはそれぞれ上下揺れ、縦揺れ、横揺れを表す。このようなデータベースが応答関数記憶部27(図2参照)に記憶される。ステップS5において船体応答関数計算部25は、ステップS3及びS4で計算された船体状態データに対応する船体応答関数を、データベースから選択する又はデータベースから選択して1次補間、2次補間等の補間計算により求める。

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

図5Aの縦軸は上下揺れの運動振幅(Amp.)/(波高/2)(Aw)であり、横軸は波長/垂線間長さ(Lpp)である。図5Bの縦軸は縦揺れの運動振幅(Amp.)/(波高/2)(Aw)であり、横軸は波長/垂線間長さ(Lpp)である。図5Cの縦軸は横揺れの運動振幅(Amp.)/(波高/2)(Aw)であり、横軸は波長/垂線間長さ(LPP)である。図5A〜図5Cでは、船体の上下揺れ、縦揺れ、横揺れそれぞれの応答関数f_xyzについて、海象条件、波高、波長(λ)、波向き(ω)をパラメータとし運動振幅の計算例を示している。

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

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

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

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

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

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

図6A〜図6Cでは、長さ70m程度の船についてステップS1で計測される波浪中の船体運動データを示している。すなわち、図6A〜図6Cは、それぞれ縦揺れ角、横揺れ角、上下揺れの変位の時系列データを示している。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1 海象推定システム 2 サテライトコンパス 3 情報処理装置(海象推定装置) 4 ディスプレイ 21 計測部 22 履歴記憶部 23 クロススペクトル計算部 24 船体状態データ計算部 25 船体応答関数計算部 26 海象推定部 27 応答関数記憶部 28 出力部