JP2018136671A

JP2018136671A - 内部波パラメータ推定装置および内部波パラメータ推定方法

Info

Publication number: JP2018136671A
Application number: JP2017029861A
Authority: JP
Inventors: 到森下; Itaru Morishita
Original assignee: Oki Electric Industry Co Ltd
Current assignee: Oki Electric Industry Co Ltd
Priority date: 2017-02-21
Filing date: 2017-02-21
Publication date: 2018-08-30

Abstract

【課題】内部波のパラメータを推定する際の演算処理を簡略化する内部波パラメータ推定装置および内部波パラメータ推定方法を提供すること。【解決手段】内部波パラメータ推定装置は、初期観測ベクトルおよび初期観測行列と、第１初期値および第２初期値ともとに観測ベクトルを生成する観測ベクトル生成部と、初期観測行列と、第１初期値および第２初期値とをもとに観測行列を生成する観測行列生成部と、を有している。また、内部波パラメータ推定装置は、観測ベクトルと観測行列とをもとに内部波パラメータの差分ベクトルを求める連立方程式求解部と、差分ベクトルと、第１初期値および第２初期値とをもとに、モード帯域および総エネルギー量を求める内部波パラメータ計算部と、を有している。【選択図】図２

Description

本発明は、海洋内部波による海中の等密度面の変位のシミュレーションにおいて、海洋内部波のパラメータを推定する内部波パラメータ推定装置および内部波パラメータ推定方法に関する。

海洋における海水は、海面から海底に向かって密度が徐々に大きくなるように成層している。ここで、海面上の風又は潮汐などの外力が働くと、海面が波立つのと同じように、海中の等密度面にも海洋内部波（以下、内部波ともいう。）と呼ばれる波が生じる。内部波は、海中の水温又は音速の分布に時間的および空間的な変動をもたらす。そして、海中の水温又は音速の分布の変動は、例えば、海中における音波伝搬のシミュレーション結果に影響を与えることが知られており、その影響の大きさを見積もるためには、内部波による等密度面の変位（以下、内部波変位ともいう。）を予測する必要がある。なお、海中の等密度面とは、海中において密度が等しい面のことである。

従来から、内部波変位を予測する際に使用される内部波のスペクトル特性は、内部波の周波数およびモード番号に基づいて求められる（例えば非特許文献１参照）。内部波のスペクトル特性は、海中での浮力周波数を用いて求めることができ、海中での浮力周波数は、現場において水温又は塩分濃度を測定することにより求められる（例えば非特許文献２参照）。また、非特許文献３には、モードスペクトルを推定する方法が開示されており、非特許文献４には、連立方程式を解く手法としてＱＲ法が開示されている。

Munk, "9.9 An inconclusive discussion," in The evolution of physical oceanography (MIT, 1981) Millard, Owens, Foronoff, "On the calculation of the Brunt-Vaisala frequency," Deep-Sea Research 37(1), 1990. Colosi, Duda, Lin, Lynch, Newhall, Cornuelle, " Observations of sound-speed fluctuations on the New Jersey continental shelf in the summer of 2006," J. Acoust. Soc. Am. 131(2), 2012. Gene H. Golub and Charles F. Van Loan, "5.3 The Full- Rank Least Squares Problem," in Matrix Computations 4th edition (The Johns Hopkins University Press, 2013)

しかしながら、非特許文献１〜４に記載された手法をもとにモードスペクトルを推定する場合、近似の誤差を小さくするためにモード数を大きくすると、モードスペクトルのベクトルを求めるための演算量が大きくなる。そして、モードスペクトルのベクトルからモード帯域および総エネルギー量といった内部波のパラメータを推定するためには、最小自乗法などによる演算処理が必要となる。こうした課題に鑑みて、従来から、内部波のパラメータを推定する際の演算処理を簡略化する手法が望まれている。

本発明に係る内部波パラメータ推定装置は、海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成部と、内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成部と、初期観測ベクトル生成部において生成された初期観測ベクトルと、初期観測行列生成部において生成された初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成部と、初期観測行列生成部において生成された初期観測行列と、第１初期値および第２初期値とをもとに観測行列を生成する観測行列生成部と、観測ベクトル生成部において生成された観測ベクトルと、観測行列生成部において生成された観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解部と、連立方程式求解部において求められた差分ベクトルと、第１初期値および第２初期値とをもとに、内部波パラメータとしてモード帯域および総エネルギー量を求める内部波パラメータ計算部と、を有するものである。

また、本発明に係る内部波パラメータ推定方法は、海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成工程と、内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成工程と、初期観測ベクトル生成工程で生成した初期観測ベクトルと、初期観測行列生成工程で生成した初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成工程と、初期観測行列生成工程で生成した初期観測行列と、第１初期値および第２初期値とをもとに観測行列を生成する観測行列生成工程と、観測ベクトル生成工程で生成した観測ベクトルと、観測行列生成工程で生成した観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解工程と、連立方程式求解工程で求めた差分ベクトルと、第１初期値および第２初期値とをもとに、内部波パラメータとしてモード帯域および総エネルギー量を求める内部波パラメータ計算工程と、を有するものである。

本発明は、初期観測ベクトルおよび初期観測行列を用いて生成した観測ベクトルおよび観測行列をもとにモード帯域および総エネルギー量を求めることから、モード数を大きくしても演算量が変化せず、かつ最小自乗法などによる演算処理が不要となるため、内部波のパラメータを推定する際の演算処理を簡略化することができる。

本発明の実施の形態に係る内部波変位シミュレーションシステムの構成を示すブロック図である。本発明の実施の形態に係る内部波パラメータ推定装置の機能的構成を示すブロック図である。図１の内部波パラメータ推定装置の動作を示すフローチャートである。従来技術に係る内部波パラメータ推定装置の構成を示すブロック図である。図２の内部波パラメータ推定装置による演算回数と、図４の内部波パラメータ推定装置による演算回数とを示すグラフである。

実施の形態．
図１は、本発明の実施の形態に係る内部波変位シミュレーションシステムの構成を示すブロック図である。内部波変位シミュレーションシステム２００は、内部波変位を予測し、海中の音波伝搬のシミュレーションなどを行うものである。ここで、内部波変位シミュレーションシステム２００の構成内容を説明する前に、内部波のパラメータを推定し、内部波変位を予測する原理について説明する。

内部波変位の予測において使用される内部波のスペクトル特性Ｆ（ω，ｊ）は、周波数ωとモード番号ｊに対して与えられた下記の式（１）によって求めることができる。式（１）のＢ（ω）は下記式（２）によって定まり、式（１）のＨ（ｊ，ｊ_＊）は下記式（３）によって定まる。式（３）のＨ_０（ｊ_＊）は下記式（４）によって定まる。モード番号ｊは、予め設定されたモード数Ｊを上限とする任意の自然数である（ｊ＝１、２、…、Ｊ）。

ここで、ｊ_＊はモード帯域を表すパラメータであり、Ｅｎは内部波のエネルギーを表すパラメータである。ｆはコリオリパラメータ（慣性周波数）であり、ｆ＝２Ωｓｉｎ（φ）によって計算される。Ωは地球の自転角速度であり、φは緯度である。ｂおよびＮ_０は、浮力周波数Ｎ（ｚ）が下記式（５）で表されるときの定数である。Ｎ_０は海面での浮力周波数であり、ｂはＮ（ｚ）の深度に対する低下度合いを決めるパラメータである。なお、深度とは海面からの深さのことである。

任意のＮ（ｚ）については、式（５）の関係から、水深Ｚが十分大きいと仮定すると、下記式（６）の関係が得られるため、本実施の形態では、式（１）を下記式（７）のように変形して使用する。

式（７）の右辺の各要素のうち、Ｎ（ｚ）は現場において水温及び塩分濃度を測定することにより求められる。また、周波数ω、コリオリパラメータｆ、およびモード番号ｊは既知の量である。

従って、スペクトル特性Ｆ（ω，ｊ）を求めるためには、Ｎ（ｚ）の深度に対する低下度合いを決めるパラメータｂと内部波のエネルギーＥｎとを乗算したｂＥｎと、モード帯域ｊ_＊とを知る必要がある。以下では、ｂＥｎを総エネルギー量Ｅ_０と称する（Ｅ_０＝ｂＥｎ）と共に、内部波のパラメータであるモード帯域ｊ_＊および総エネルギー量Ｅ_０を内部波パラメータともいう。

ここで、式（７）の右辺の既知量の項をＫ（ω）とし、式（７）の未知量を含む項をＸ（ｊ，ｊ_＊，Ｅ_０）として、式（７）を下記式（８）のように書き換える。以下では、Ｘ（ｊ，ｊ_＊，Ｅ_０）をモードスペクトルという。式（８）のＫ（ω）は下記式（９）によって定まり、式（８）のＸ（ｊ，ｊ_＊，Ｅ_０）は下記式（１０）によって定まる。

ここで、モードスペクトルＸ（ｊ，ｊ_＊，Ｅ_０）は、以下のように推定することができる。内部波変位の周波数ωの成分が、深度ｚ_１，ｚ_２，・・・，ｚ_ＬでそれぞれＰ個測定されたとし、各深度に対応するＰ個の成分を、それぞれ、ζ_ｐ（ｚ_１），ζ_ｐ（ｚ_２），・・・，ζ_ｐ（ｚ_Ｌ）（ｐ＝１，２，…，Ｐ）とする。すると、異なる深度に対する内部波変位の積の期待値は、下記式（１１）のように近似することができる。測定深度数Ｌは、内部波変位を測定する深度の数である。変位数Ｐは、各深度に対応する内部波変位の個数である。

ここで、Ｊは内部波変位に影響を与えているモード数であり、十分大きな数を与えることにより、式（１１）による近似の誤差を小さくすることができる。また、Ｗ_ｊ（ｚ）は、下記式（１２）および式（１３）の微分方程式を満足する内部波モードであり、α_ｊはｊ番目の内部波モードＷ_ｊ（ｚ）に対応する固有波数である。

ｍ＝１，２，・・・，Ｌ、およびｎ＝１，２，・・・，Ｌに対して式（１１）を連立すると、初期観測ベクトルをｙ、初期観測行列をＥ、未知のモードスペクトルのベクトルをｘとした下記式（１４）の連立方程式が得られる。式（１４）の初期観測ベクトルｙは下記式（１５）によって定まり、式（１４）の初期観測行列Ｅは下記式（１６）によって定まり、式（１４）のモードスペクトルのベクトルｘは下記式（１７）によって定まる。

式（１４）を一般的な線形逆問題の手法によって解くことにより、モードスペクトルのベクトルｘの値を得ることができる。

また、式（１０）で表されるモードスペクトルＸ（ｊ，ｊ_＊，Ｅ_０）をモード帯域ｊ_＊および総エネルギー量Ｅ_０の関数とみて、モード帯域の初期値ｊ_＊ ^（０）および総エネルギー量の初期値Ｅ_０ ^（０）の周りでテーラー展開し、１次の項までをとると、「ｊ_＊＝ｊ_＊ ^（０）＋Δｊ_＊」および「Ｅ_０＝Ｅ_０ ^（０）＋ΔＥ_０」におけるＸ（ｊ，ｊ_＊，Ｅ_０）の値は、下記式（１８）のように近似することができる。以降では、モード帯域ｊ_＊の初期値ｊ_＊ ^（０）を単に「初期値ｊ_＊ ^（０）」ともいい、総エネルギー量Ｅ_０の初期値Ｅ_０ ^（０）を単に「初期値Ｅ_０ ^（０）」ともいう。なお、初期値ｊ_＊ ^（０）は本発明の「第１初期値」に相当し、初期値Ｅ_０ ^（０）は本発明の「第２初期値」に相当する。

ここで、ｊ_＊による偏微分は下記式（１９）のようになる。式（１９）のＨ’_０（ｊ_＊）は下記式（２０）で表され、式（２０）の最後の無限級数は下記式（２１）の公式を使って計算することができる。一方、Ｅ_０による偏微分は、下記式（２２）のようになる。

ここで、式（１８）を式（１７）に代入すると、式（１４）の連立方程式は、下記式（２３）のように変形することができる。式（２３）のｘ_０は下記式（２４）によって定まる内部波パラメータの初期ベクトルであり、式（２３）のＤは下記式（２５）によって定まる微分ベクトルである。

ここで、式（２３）を整理すると、下記式（２６）のように、観測ベクトルをΔｙ、観測行列をＭ、未知の内部波パラメータの差分ベクトルをΔｘとする新たな連立方程式が得られる。観測ベクトルΔｙは下記式（２７）によって定まり、観測行列Ｍは下記式（２８）によって定まり、差分ベクトルΔｘは下記式（２９）によって定まる。

式（２６）を一般的な連立方程式の解法によって解き、内部波パラメータの差分ベクトルΔｘ＝［Δｊ_＊ ΔＥ_０］^Ｔ（Ｔは行列の転置を表す記号である）を得ることができたら、モード帯域ｊ_＊は下記式（３０）によって求めることができ、総エネルギー量Ｅ_０は下記式（３１）によって求めることができる。

なお、初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）と、真のモード帯域および真の総エネルギー量との差が大きいほど、式（１８）の線形近似による誤差が大きくなる。初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）は、例えば、下記式（３２）および式（３３）に示す典型値に設定することができる。式（３２）および式（３３）に示す初期値は、真値と大差ないため、線形近似の誤差を小さくすることができる。

ところで、式（３０）および式（３１）により求めたモード帯域ｊ_＊および総エネルギー量Ｅ_０は、式（３２）および式（３３）に示す初期値よりも真値に近づいていると考えられる。そのため、式（３０）により求めたモード帯域ｊ_＊を新たな初期値ｊ_＊ ^（０）とし、式（３１）により求めた総エネルギー量Ｅ_０を新たな初期値Ｅ_０ ^（０）として、式（２６）の連立方程式を構成するようにしてもよい。このように再構築した式（２６）の連立方程式を解くようにすれば、線形近似の誤差をさらに小さくすることができる。

続いて、図１を参照して、内部波変位シミュレーションシステム２００の構成について説明する。図１に示すように、内部波変位シミュレーションシステム２００は、内部波変位測定装置１１０と、内部波パラメータ推定装置１００と、スペクトル特性演算装置１２０と、内部波変位予測処理装置１３０と、を有している。

内部波変位測定装置１１０は、複数の深度ｚ_１，ｚ_２，・・・，ｚ_Ｌのそれぞれについて、Ｐ個の内部波変位であるζ_ｐ（ｚ_１），ζ_ｐ（ｚ_２），・・・，ζ_ｐ（ｚ_Ｌ）（ｐ＝１，２，…，Ｐ）を測定するものである。本実施の形態において、測定深度数Ｌは２以上の任意の自然数であり、変位数Ｐは２以上の任意の自然数であるものとする。すなわち、内部波変位測定装置１１０は、深度の異なる複数の地点において、複数の内部波変位を求めるものである。内部波変位測定装置１１０が求める内部波変位は、内部波変位の周波数ωの成分であり、周波数ωは予め設定されている。内部波変位測定装置１１０は、内部波変位の測定深度である深度ｚ_１，ｚ_２，・・・，ｚ_Ｌと、測定した内部波変位ζ_ｐ（ｚ_１），ζ_ｐ（ｚ_２），・・・，ζ_ｐ（ｚ_Ｌ）と、内部波変位の測定に係る周波数ωとを、内部波パラメータ推定装置１００へ出力するものである。

内部波パラメータ推定装置１００は、内部波変位測定装置１１０から入力された内部波変位の測定深度、複数の内部波変位、および周波数ωなどをもとに、内部波のパラメータとして、モード帯域ｊ_＊および総エネルギー量Ｅ_０を求めるものである。また、内部波パラメータ推定装置１００は、求めたモード帯域ｊ_＊および総エネルギー量Ｅ_０をスペクトル特性演算装置１２０へ出力するものである。

スペクトル特性演算装置１２０は、内部波パラメータ推定装置１００から出力されたモード帯域ｊ_＊および総エネルギー量Ｅ_０を用いて、式（８）により、スペクトル特性Ｆ（ω，ｊ）を求めるものである。また、スペクトル特性演算装置１２０は、求めたスペクトル特性Ｆ（ω，ｊ）を内部波変位予測処理装置１３０へ出力するものである。

内部波変位予測処理装置１３０は、スペクトル特性演算装置１２０から出力されたスペクトル特性Ｆ（ω，ｊ）をもとに、内部波変位をシミュレーションして予測するものである。また、内部波変位予測処理装置１３０は、内部波変位の予測の結果をもとに、海中の音波伝搬のシミュレーションなどを行うものである。

図２は、本発明の実施の形態に係る内部波パラメータ推定装置の機能的構成を示すブロック図である。図２に示すように、内部波パラメータ推定装置１００は、初期観測ベクトル生成部１０と、初期観測行列生成部２０と、観測ベクトル生成部３０と、観測行列生成部４０と、連立方程式求解部５０と、内部波パラメータ計算部６０と、記憶部７０と、を有している。

記憶部７０には、浮力周波数Ｎ（ｚ）およびコリオリパラメータｆが予め記憶されている。また、記憶部７０には、モード帯域の初期値ｊ_＊ ^（０）と総エネルギー量の初期値Ｅ_０ ^（０）が予め格納されている。

初期観測ベクトル生成部１０には、内部波変位測定装置１１０から、深度ｚ_１，ｚ_２，・・・，ｚ_Ｌのそれぞれにおいて測定された内部波変位ζ_ｐ（ｚ_１），ζ_ｐ（ｚ_２），・・・，ζ_ｐ（ｚ_Ｌ）（ｐ＝１，２，…，Ｐ）が入力される。

初期観測ベクトル生成部１０は、内部波変位ζ_ｐ（ｚ_１），ζ_ｐ（ｚ_２），・・・，ζ_ｐ（ｚ_Ｌ）を用いて、式（１１）および式（１５）により、初期観測ベクトルｙを求めるものである。また、初期観測ベクトル生成部１０は、求めた初期観測ベクトルｙを観測ベクトル生成部３０へ出力するものである。

初期観測行列生成部２０には、内部波変位測定装置１１０から、内部波変位の測定深度である深度ｚ_１，ｚ_２，・・・，ｚ_Ｌと、周波数ωとが入力される。また、初期観測行列生成部２０は、コリオリパラメータｆ、および浮力周波数Ｎ（ｚ）を記憶部７０から取得するものである。初期観測行列生成部２０は、入力された周波数ωと、記憶部７０から取得したコリオリパラメータｆおよび浮力周波数Ｎ（ｚ）とを用いて、式（２）および式（９）により、式（７）の右辺の既知量の項Ｋ（ω）を計算するものである。

また、初期観測行列生成部２０は、内部波変位の測定深度、周波数ω、コリオリパラメータｆ、および浮力周波数Ｎ（ｚ）を用いて、式（１２）および式（１３）を満足する内部波モードＷ_ｊ（ｚ_ｍ）（ｊ＝１，２，・・・，Ｊ、ｍ＝１，２，・・・，Ｌ）を求めるものである。初期観測行列生成部２０は、一般的な微分方程式の解法により、内部波モードＷ_ｊ（ｚ_ｍ）を求めるように構成されている。

さらに、初期観測行列生成部２０は、既知量の項Ｋ（ω）と内部波モードＷ_ｊ（ｚ_ｍ）（ｊ＝１，２，・・・，Ｊ、ｍ＝１，２，・・・，Ｌ）とを用いて、式（１６）により、初期観測行列Ｅを求めるものである。そして、初期観測行列生成部２０は、求めた初期観測行列Ｅを、観測ベクトル生成部３０および観測行列生成部４０に出力するものである。

観測ベクトル生成部３０には、初期観測ベクトル生成部１０から初期観測ベクトルｙが入力され、初期観測行列生成部２０から初期観測行列Ｅが入力される。また、観測ベクトル生成部３０は、モード帯域の初期値ｊ_＊ ^（０）および総エネルギー量の初期値Ｅ_０ ^（０）を記憶部７０から取得するものである。

観測ベクトル生成部３０は、記憶部７０から取得したモード帯域の初期値ｊ_＊ ^（０）および総エネルギー量の初期値Ｅ_０ ^（０）を用いて、式（３）、式（４）、式（１０）、および式（２４）により、内部波パラメータの初期ベクトルｘ_０を求めるものである。また、観測ベクトル生成部３０は、求めた初期ベクトルｘ_０と、入力された初期観測ベクトルｙおよび初期観測行列Ｅとを用いて、式（２７）により、観測ベクトルΔｙを求めるものである。そして、観測ベクトル生成部３０は、求めた観測ベクトルΔｙを連立方程式求解部５０に出力するものである。

観測行列生成部４０には、初期観測行列生成部２０から初期観測行列Ｅが入力される。また、観測行列生成部４０は、モード帯域の初期値ｊ_＊ ^（０）、および総エネルギー量の初期値Ｅ_０ ^（０）を記憶部７０から取得するものである。

観測行列生成部４０は、記憶部７０から取得したモード帯域の初期値ｊ_＊ ^（０）、および総エネルギー量の初期値Ｅ_０ ^（０）を用いて、式（４）、式（１９）、式（２０）、式（２１）、式（２２）、および式（２５）により、微分ベクトルＤを求めるものである。また、観測行列生成部４０は、求めた微分ベクトルＤと、入力された初期観測行列Ｅとを用いて、式（２８）により、観測行列Ｍを求めるものである。そして、観測行列生成部４０は、求めた観測行列Ｍを連立方程式求解部５０に出力するものである。

連立方程式求解部５０には、観測ベクトル生成部３０から観測ベクトルΔｙが入力され、観測行列生成部４０から観測行列Ｍが入力される。連立方程式求解部５０は、入力された観測ベクトルΔｙおよび観測行列Ｍを用いて、式（２６）により、内部波パラメータの差分ベクトルΔｘを求めるものである。連立方程式求解部５０は、ＱＲ法などの一般的な連立方程式の解法によって式（２６）を解き、内部波パラメータの差分ベクトルΔｘを求めるように構成されている。また、連立方程式求解部５０は、求めた内部波パラメータの差分ベクトルΔｘを、内部波パラメータ計算部６０に出力するものである。

内部波パラメータ計算部６０には、連立方程式求解部５０から内部波パラメータの差分ベクトルΔｘが入力される。また、内部波パラメータ計算部６０は、モード帯域の初期値ｊ_＊ ^（０）および総エネルギー量の初期値Ｅ_０ ^（０）を記憶部７０から取得するものである。内部波パラメータ計算部６０は、入力された内部波パラメータの差分ベクトルΔｘと、取得したモード帯域の初期値ｊ_＊ ^（０）および総エネルギー量の初期値Ｅ_０ ^（０）とを用いて、式（３０）および式（３１）により、モード帯域ｊ_＊および総エネルギー量Ｅ_０を内部波パラメータとして求めるものである。そして、内部波パラメータ計算部６０は、求めた内部波パラメータをスペクトル特性演算装置１２０に出力する。

ここで、内部波パラメータ計算部６０は、求めた内部波パラメータによって、記憶部７０内のモード帯域の初期値ｊ_＊ ^（０）および総エネルギー量の初期値Ｅ_０ ^（０）を更新するようにしてもよい。すなわち、内部波パラメータ計算部６０は、式（３０）により推定したモード帯域ｊ_＊を初期値ｊ_＊ ^（０）に設定すると共に、式（３１）により推定した総エネルギー量Ｅ_０を初期値Ｅ_０ ^（０）に設定するようにしてもよい。このように、予め設定された初期値よりも真値に近いと考えられるモード帯域ｊ_＊および総エネルギー量Ｅ_０によって初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）を更新すれば、線形近似の誤差をさらに小さくすることができる。

ところで、上述した原理の説明では、初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）として、式（３２）および（３３）の値を与える場合を例示したが、これに限らず、他に真値に近いと考えられる値があれば、その値を初期値として用いてもよい。

ここで、初期観測ベクトル生成部１０、初期観測行列生成部２０、観測ベクトル生成部３０、観測行列生成部４０、連立方程式求解部５０、および内部波パラメータ計算部６０は、それぞれの各機能を実現する回路デバイスなどのハードウェアで実現することもできるし、例えば、マイコン、ＤＳＰ（ＤｉｇｉｔａｌＳｉｇｎａｌＰｒｏｃｅｓｓｏｒ）、又はＣＰＵ（ＣｅｎｔｒａｌＰｒｏｃｅｓｓｉｎｇＵｎｉｔ）等の演算装置上で実行されるソフトウェアとして実現することもできる。後者の場合、記憶部７０には、初期観測ベクトル生成部１０、初期観測行列生成部２０、観測ベクトル生成部３０、観測行列生成部４０、連立方程式求解部５０、および内部波パラメータ計算部６０の動作プログラムが格納されていてもよい。また、記憶部７０は、ＲＡＭ（ＲａｎｄｏｍＡｃｃｅｓｓＭｅｍｏｒｙ）、フラッシュメモリ等のＰＲＯＭ（ＰｒｏｇｒａｍｍａｂｌｅＲＯＭ）、又はＨＤＤ（ＨａｒｄＤｉｓｋＤｒｉｖｅ）等により構成することができる。

図３は、内部波パラメータ推定装置の動作を示すフローチャートである。図３に基づき、内部波パラメータ推定装置１００による内部波パラメータ推定方法について説明する。

初期観測ベクトル生成部１０は、海中における複数の深度それぞれの内部波変位をもとに初期観測ベクトルｙを生成し、生成した初期観測ベクトルｙを観測ベクトル生成部３０へ出力する（ステップＳ１０１／初期観測ベクトル生成工程）。

また、初期観測行列生成部２０は、内部波変位の測定に係る周波数ωと、コリオリパラメータｆと、浮力周波数Ｎ（ｚ）とを用いて、式（７）の右辺の既知量の項Ｋ（ω）を計算する。次いで、初期観測行列生成部２０は、内部波変位の測定深度、周波数ω、コリオリパラメータｆ、および浮力周波数Ｎ（ｚ）を用いて、内部波モードＷ_ｊ（ｚ_ｍ）を求める。そして、初期観測行列生成部２０は、既知量の項Ｋ（ω）と内部波モードＷ_ｊ（ｚ_ｍ）とをもとに初期観測行列Ｅを生成し、生成した初期観測行列Ｅを観測ベクトル生成部３０および観測行列生成部４０に出力する（ステップＳ１０２／初期観測行列生成工程）。

次に、観測ベクトル生成部３０は、初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）を用いて内部波パラメータの初期ベクトルｘ_０を求める。そして、観測ベクトル生成部３０は、求めた初期ベクトルｘ_０と、初期観測ベクトルｙおよび初期観測行列Ｅとをもとに観測ベクトルΔｙを求めて連立方程式求解部５０に出力する（ステップＳ１０３／観測ベクトル生成工程）。

また、観測行列生成部４０は、初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）を用いて微分ベクトルＤを求める。そして、観測行列生成部４０は、求めた微分ベクトルＤと、初期観測行列Ｅとをもとに観測行列Ｍを求めて連立方程式求解部５０に出力する（ステップＳ１０４／観測行列生成工程）。

連立方程式求解部５０は、観測ベクトルΔｙおよび観測行列Ｍを用いて内部波パラメータの差分ベクトルΔｘを求め、求めた差分ベクトルΔｘを内部波パラメータ計算部６０に出力する（ステップＳ１０５／連立方程式求解工程）。

続いて、内部波パラメータ計算部６０は、内部波パラメータの差分ベクトルΔｘと、初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）とを用いてモード帯域ｊ_＊および総エネルギー量Ｅ_０を求めるものである。そして、内部波パラメータ計算部６０は、求めたモード帯域ｊ_＊および総エネルギー量Ｅ_０を内部波パラメータとしてスペクトル特性演算装置１２０に出力する（ステップＳ１０６／内部波パラメータ計算工程）。

また、内部波パラメータ計算部６０は、内部波パラメータ計算工程で求めたモード帯域ｊ_＊および総エネルギー量Ｅ_０によって、記憶部７０内の初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）を更新する（ステップＳ１０７／初期値更新工程）。そして、内部波パラメータ推定装置１００は、上記ステップＳ１０３〜Ｓ１０７の一連の処理を繰り返し実行する。

ここで、上記の動作説明は、図３に付したステップ番号の順に説明したが、これに限定されるものではない。例えば、ステップＳ１０１の初期観測ベクトル生成工程とステップＳ１０２の初期観測行列生成工程とは並行して行うことができ、ステップＳ１０３の観測ベクトル生成工程とステップＳ１０４の観測行列生成工程とも並行して行うことができる。

以上のように、本実施の形態における内部波パラメータ推定装置１００は、初期観測ベクトルｙおよび初期観測行列Ｅを用いて生成した観測ベクトルΔｙおよび観測行列Ｍをもとにモード帯域ｊ_＊および総エネルギー量Ｅ_０を求める。そのため、モード数Ｊを大きくしても演算量が変化せず、かつ最小自乗法などによる演算処理が不要となることから、内部波のパラメータを推定する際の演算処理を簡略化することができる。したがって、内部波変位シミュレーションシステム２００によれば、少ない演算量で内部波変位を精度よくシミュレーションすることができる。

また、内部波パラメータ計算部６０には、求めたモード帯域ｊ_＊および総エネルギー量Ｅ_０によって、初期値ｊ_＊ ^（０）および初期値Ｅ_０ ^（０）を更新する機能をもたせることができる。このようにすれば、設定された初期値をさらに真値に近づけることができるため、線形近似の誤差をより小さくすることができる。

さらに、初期観測ベクトル生成部１０は、海中における複数の深度それぞれの内部波変位をもとに初期観測ベクトルｙを生成するものであり、初期観測行列生成部２０は、複数の深度を用いて初期観測行列Ｅを生成するものである。すなわち、本実施の形態では、測定深度数Ｌが２以上の自然数に設定されている。したがって、より広域の内部波変位を用いることができることから、内部波パラメータの推定精度が向上するため、内部波変位のシミュレーションの精度向上を図ることができる。

ここで、図４および図５を参照して、内部波パラメータ推定装置１００によって得られる効果を詳細に説明する。図４は、従来技術に係る内部波パラメータ推定装置の構成を示すブロック図である。図５は、図２の内部波パラメータ推定装置による演算回数と、図４の内部波パラメータ推定装置による演算回数とを示すグラフである。図５では、横軸にモード数Ｊをとり、縦軸に浮動小数点演算回数をとっている。なお、図５には、内部波変位の測定深度数Ｌが３０の場合を例示している。

図４に示すように、従来の構成である内部波パラメータ推定装置１０００は、初期観測ベクトルｙを生成する初期観測ベクトル生成部１１００と、初期観測行列Ｅを生成する初期観測行列生成部１２００と、初期観測ベクトルｙと初期観測行列Ｅとからモードスペクトルのベクトルｘを求める連立方程式求解部１５００と、を有している。

すなわち、内部波パラメータ推定装置１０００の連立方程式求解部１５００には、初期観測ベクトルｙおよび初期観測行列Ｅが入力されるようになっている。そして、連立方程式求解部１５００は、初期観測ベクトルｙおよび初期観測行列Ｅをもとに、ＱＲ法などの一般的な連立方程式の解法によって式（１４）を解き、モードスペクトルのベクトルｘを出力するようになっている。そのため、モード帯域ｊ_＊および総エネルギー量Ｅ_０を推定するためには、モードスペクトルのベクトルｘが得られた後に、最小自乗法などの演算処理を行う必要がある。

ここで、例えば、下記式（３４）のように誤差関数を定義する。ｘ_ｊは、求めたモードスペクトルのベクトルｘのｊ番目の要素である。

つまり、従来の構成では、式（３４）の誤差関数が最小となるように、モード帯域ｊ_＊および総エネルギー量Ｅ_０を推定する必要がある。

また、式（１４）および式（１６）の初期観測行列Ｅは、行数がＬ（Ｌ＋１）／２であり、列数がＪである。したがって、式（１１）における近似の誤差を小さくするためにモード数Ｊを大きくすると、モードスペクトルのベクトルｘを求めるための演算量が大きくなる。

ここで、例えば、連立方程式を解く手法として広く使われているＱＲ法を用いる場合を考える。ＱＲ法は、初期観測行列Ｅが「Ｅ＝ＱＲ」のように対角行列Ｑと上三角行列Ｒとの積で表せることを利用するものであり、浮動小数点演算が約２Ｊ^２［Ｌ（Ｌ＋１）／２−Ｊ／３］回必要であることが知られている。

図５のように、測定深度数Ｌが３０の場合、従来の構成では、浮動小数点演算が約２Ｊ^２［４６５−Ｊ／３］回必要となる。すなわち、従来の構成では、図５の破線グラフＧｙのように、モード数Ｊを大きくするにつれて演算量が増加していくことがわかる。

一方、本実施の形態では、式（２６）の連立方程式において、観測行列Ｍの列数は常に２である。よって、式（２６）の連立方程式をＱＲ法で解く場合には、浮動小数点演算の回数は約８×［Ｌ（Ｌ＋１）／２−２／３］回となる。

ここで、線形近似の誤差は、浮動小数点演算を繰り返し行う毎に小さくなり、数回の繰り返しによって十分に小さくなる。したがって、本実施の形態では、浮動小数点演算を繰り返す回数である繰り返し回数Ｓが２以上に設定されている。この場合、浮動小数点演算の回数は約８×Ｓ×［Ｌ（Ｌ＋１）／２−２／３］回となる。図５では、繰り返し回数Ｓが１０に設定されている場合を例示しており、この場合、浮動小数点演算の回数は約８０×［Ｌ（Ｌ＋１）／２−２／３］回となる。

すなわち、従来の構成では、図５の破線グラフＧｙのように、モード数Ｊを大きくするほど、浮動小数点演算の回数が増加していたが、内部波パラメータ推定装置１００では、図５の実線グラフＧｘのように、浮動小数点演算の回数がモード数Ｊに依存しない。つまり、内部波パラメータ推定装置１００においては、モード数Ｊを大きくしても、浮動小数点演算の回数が多くならない。

したがって、式（１１）における近似の誤差を小さくするために、モード数Ｊを大きくすることを考えるとき、例えば、繰り返し回数Ｓを１０に設定した場合は、モード数を７以上にすれば、従来の構成よりも少ない演算回数で内部波パラメータを推定することができる。なお、非特許文献３では、モード数Ｊを１５としており、かかる設定の場合、図５に演算回数差ΔＣとして示すように、内部波パラメータ推定装置１００と従来の構成との間に、演算回数の大きな違いが生じる。つまり、内部波パラメータ推定装置１００によれば、モード数Ｊを大きくすればするほど、連立方程式を解くための演算回数が従来よりも削減され、優位性が向上する。

また、内部波パラメータ推定装置１００は、式（２６）の連立方程式「Δｙ＝ＭΔｘ」を解くことにより、内部波パラメータを得ることができるため、従来の構成のように、最小自乗法などによる内部波パラメータの推定処理を行う必要がない。すなわち、内部波パラメータ推定装置１００によれば、内部波のパラメータを推定する際の演算処理を簡略化することができる。

なお、上述した実施の形態は、内部波パラメータ推定装置および内部波パラメータ推定方法における好適な具体例であり、本発明の技術的範囲は、これらの態様に限定されるものではない。例えば、上記実施の形態では、測定深度数Ｌが２以上の自然数である場合を例示したが、これに限らず、測定深度数Ｌは１であってもよい。ただし、測定深度数Ｌは、大きな値に設定した方が、より広域の内部波変位を用いることができるため、内部波パラメータの推定精度が向上することから、内部波変位のシミュレーションの精度向上を図ることができる。また、上記実施の形態では、変位数Ｐが２以上の自然数である場合を例示したが、これに限らず、変位数Ｐは１であってもよい。ただし、変位数Ｐについても、大きな値に設定した方が、より広域の内部波変位を用いることができるため、内部波パラメータの推定精度が向上することから、内部波変位のシミュレーションの精度向上を図ることができる。

１０初期観測ベクトル生成部、２０初期観測行列生成部、３０観測ベクトル生成部、４０観測行列生成部、５０連立方程式求解部、６０内部波パラメータ計算部、７０記憶部、１００内部波パラメータ推定装置、１１０内部波変位測定装置、１２０スペクトル特性演算装置、１３０内部波変位予測処理装置、２００内部波変位シミュレーションシステム。

Claims

海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成部と、
前記内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成部と、
前記初期観測ベクトル生成部において生成された前記初期観測ベクトルと、前記初期観測行列生成部において生成された前記初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成部と、
前記初期観測行列生成部において生成された前記初期観測行列と、前記第１初期値および前記第２初期値とをもとに観測行列を生成する観測行列生成部と、
前記観測ベクトル生成部において生成された前記観測ベクトルと、前記観測行列生成部において生成された前記観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解部と、
前記連立方程式求解部において求められた前記差分ベクトルと、前記第１初期値および前記第２初期値とをもとに、前記内部波パラメータとして前記モード帯域および前記総エネルギー量を求める内部波パラメータ計算部と、を有する内部波パラメータ推定装置。
前記内部波パラメータ計算部は、
求めた前記モード帯域および前記総エネルギー量によって、前記第１初期値および前記第２初期値を更新するものである請求項１に記載の内部波パラメータ推定装置。
前記初期観測ベクトル生成部は、
海中における複数の深度それぞれの前記内部波変位をもとに前記初期観測ベクトルを生成するものであり、
前記初期観測行列生成部は、
前記複数の深度を用いて前記初期観測行列を生成するものである請求項１又は２に記載の内部波パラメータ推定装置。
海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成工程と、
前記内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成工程と、
前記初期観測ベクトル生成工程で生成した前記初期観測ベクトルと、前記初期観測行列生成工程で生成した前記初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成工程と、
前記初期観測行列生成工程で生成した前記初期観測行列と、前記第１初期値および前記第２初期値とをもとに観測行列を生成する観測行列生成工程と、
前記観測ベクトル生成工程で生成した前記観測ベクトルと、前記観測行列生成工程で生成した前記観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解工程と、
前記連立方程式求解工程で求めた前記差分ベクトルと、前記第１初期値および前記第２初期値とをもとに、前記内部波パラメータとして前記モード帯域および前記総エネルギー量を求める内部波パラメータ計算工程と、を有する内部波パラメータ推定方法。