(57)【要約】【課題】内部波のパラメータを推定する際の演算処理を簡略化する内部波パラメータ推定装置および内部波パラメータ推定方法を提供すること。【解決手段】内部波パラメータ推定装置は、初期観測ベクトルおよび初期観測行列と、第１初期値および第２初期値ともとに観測ベクトルを生成する観測ベクトル生成部と、初期観測行列と、第１初期値および第２初期値とをもとに観測行列を生成する観測行列生成部と、を有している。また、内部波パラメータ推定装置は、観測ベクトルと観測行列とをもとに内部波パラメータの差分ベクトルを求める連立方程式求解部と、差分ベクトルと、第１初期値および第２初期値とをもとに、モード帯域および総エネルギー量を求める内部波パラメータ計算部と、を有している。【選択図】図２
2018136671.tif 000002
【特許請求の範囲】
【請求項１】
海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成部と、
前記内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成部と、
前記初期観測ベクトル生成部において生成された前記初期観測ベクトルと、前記初期観測行列生成部において生成された前記初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成部と、
前記初期観測行列生成部において生成された前記初期観測行列と、前記第１初期値および前記第２初期値とをもとに観測行列を生成する観測行列生成部と、
前記観測ベクトル生成部において生成された前記観測ベクトルと、前記観測行列生成部において生成された前記観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解部と、
前記連立方程式求解部において求められた前記差分ベクトルと、前記第１初期値および前記第２初期値とをもとに、前記内部波パラメータとして前記モード帯域および前記総エネルギー量を求める内部波パラメータ計算部と、を有する内部波パラメータ推定装置。
【請求項２】
前記内部波パラメータ計算部は、
求めた前記モード帯域および前記総エネルギー量によって、前記第１初期値および前記第２初期値を更新するものである請求項１に記載の内部波パラメータ推定装置。
【請求項３】
前記初期観測ベクトル生成部は、
海中における複数の深度それぞれの前記内部波変位をもとに前記初期観測ベクトルを生成するものであり、
前記初期観測行列生成部は、
前記複数の深度を用いて前記初期観測行列を生成するものである請求項１又は２に記載の内部波パラメータ推定装置。
【請求項４】
海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成工程と、
前記内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成工程と、
前記初期観測ベクトル生成工程で生成した前記初期観測ベクトルと、前記初期観測行列生成工程で生成した前記初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成工程と、
前記初期観測行列生成工程で生成した前記初期観測行列と、前記第１初期値および前記第２初期値とをもとに観測行列を生成する観測行列生成工程と、
前記観測ベクトル生成工程で生成した前記観測ベクトルと、前記観測行列生成工程で生成した前記観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解工程と、
前記連立方程式求解工程で求めた前記差分ベクトルと、前記第１初期値および前記第２初期値とをもとに、前記内部波パラメータとして前記モード帯域および前記総エネルギー量を求める内部波パラメータ計算工程と、を有する内部波パラメータ推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、海洋内部波による海中の等密度面の変位のシミュレーションにおいて、海洋内部波のパラメータを推定する内部波パラメータ推定装置および内部波パラメータ推定方法に関する。
【背景技術】
【０００２】
海洋における海水は、海面から海底に向かって密度が徐々に大きくなるように成層している。ここで、海面上の風又は潮汐などの外力が働くと、海面が波立つのと同じように、海中の等密度面にも海洋内部波（以下、内部波ともいう。）と呼ばれる波が生じる。内部波は、海中の水温又は音速の分布に時間的および空間的な変動をもたらす。そして、海中の水温又は音速の分布の変動は、例えば、海中における音波伝搬のシミュレーション結果に影響を与えることが知られており、その影響の大きさを見積もるためには、内部波による等密度面の変位（以下、内部波変位ともいう。）を予測する必要がある。なお、海中の等密度面とは、海中において密度が等しい面のことである。
【０００３】
従来から、内部波変位を予測する際に使用される内部波のスペクトル特性は、内部波の周波数およびモード番号に基づいて求められる（例えば非特許文献１参照）。内部波のスペクトル特性は、海中での浮力周波数を用いて求めることができ、海中での浮力周波数は、現場において水温又は塩分濃度を測定することにより求められる（例えば非特許文献２参照）。また、非特許文献３には、モードスペクトルを推定する方法が開示されており、非特許文献４には、連立方程式を解く手法としてＱＲ法が開示されている。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】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)
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、非特許文献１〜４に記載された手法をもとにモードスペクトルを推定する場合、近似の誤差を小さくするためにモード数を大きくすると、モードスペクトルのベクトルを求めるための演算量が大きくなる。そして、モードスペクトルのベクトルからモード帯域および総エネルギー量といった内部波のパラメータを推定するためには、最小自乗法などによる演算処理が必要となる。こうした課題に鑑みて、従来から、内部波のパラメータを推定する際の演算処理を簡略化する手法が望まれている。
【課題を解決するための手段】
【０００６】
本発明に係る内部波パラメータ推定装置は、海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成部と、内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成部と、初期観測ベクトル生成部において生成された初期観測ベクトルと、初期観測行列生成部において生成された初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成部と、初期観測行列生成部において生成された初期観測行列と、第１初期値および第２初期値とをもとに観測行列を生成する観測行列生成部と、観測ベクトル生成部において生成された観測ベクトルと、観測行列生成部において生成された観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解部と、連立方程式求解部において求められた差分ベクトルと、第１初期値および第２初期値とをもとに、内部波パラメータとしてモード帯域および総エネルギー量を求める内部波パラメータ計算部と、を有するものである。
【０００７】
また、本発明に係る内部波パラメータ推定方法は、海洋内部波による等密度面の変位である内部波変位をもとに初期観測ベクトルを生成する初期観測ベクトル生成工程と、内部波変位の測定に係る周波数とコリオリパラメータと浮力周波数とをもとに初期観測行列を生成する初期観測行列生成工程と、初期観測ベクトル生成工程で生成した初期観測ベクトルと、初期観測行列生成工程で生成した初期観測行列と、モード帯域の初期値である第１初期値および総エネルギー量の初期値である第２初期値ともとに観測ベクトルを生成する観測ベクトル生成工程と、初期観測行列生成工程で生成した初期観測行列と、第１初期値および第２初期値とをもとに観測行列を生成する観測行列生成工程と、観測ベクトル生成工程で生成した観測ベクトルと、観測行列生成工程で生成した観測行列とをもとに、内部波パラメータの差分ベクトルを求める連立方程式求解工程と、連立方程式求解工程で求めた差分ベクトルと、第１初期値および第２初期値とをもとに、内部波パラメータとしてモード帯域および総エネルギー量を求める内部波パラメータ計算工程と、を有するものである。
【発明の効果】
【０００８】
本発明は、初期観測ベクトルおよび初期観測行列を用いて生成した観測ベクトルおよび観測行列をもとにモード帯域および総エネルギー量を求めることから、モード数を大きくしても演算量が変化せず、かつ最小自乗法などによる演算処理が不要となるため、内部波のパラメータを推定する際の演算処理を簡略化することができる。
【図面の簡単な説明】
【０００９】
【図１】本発明の実施の形態に係る内部波変位シミュレーションシステムの構成を示すブロック図である。
【図２】本発明の実施の形態に係る内部波パラメータ推定装置の機能的構成を示すブロック図である。
【図３】図１の内部波パラメータ推定装置の動作を示すフローチャートである。
【図４】従来技術に係る内部波パラメータ推定装置の構成を示すブロック図である。
【図５】図２の内部波パラメータ推定装置による演算回数と、図４の内部波パラメータ推定装置による演算回数とを示すグラフである。
【発明を実施するための形態】
【００１０】
実施の形態．
図１は、本発明の実施の形態に係る内部波変位シミュレーションシステムの構成を示すブロック図である。内部波変位シミュレーションシステム２００は、内部波変位を予測し、海中の音波伝搬のシミュレーションなどを行うものである。ここで、内部波変位シミュレーションシステム２００の構成内容を説明する前に、内部波のパラメータを推定し、内部波変位を予測する原理について説明する。
【００１１】
内部波変位の予測において使用される内部波のスペクトル特性Ｆ（ω，ｊ）は、周波数ωとモード番号ｊに対して与えられた下記の式（１）によって求めることができる。式（１）のＢ（ω）は下記式（２）によって定まり、式（１）のＨ（ｊ，ｊ＊）は下記式（３）によって定まる。式（３）のＨ０（ｊ＊）は下記式（４）によって定まる。モード番号ｊは、予め設定されたモード数Ｊを上限とする任意の自然数である（ｊ＝１、２、…、Ｊ）。
【００１２】
【数１】
2018136671.tif 000003
【００１３】
【数２】
2018136671.tif 000004
【００１４】
【数３】
2018136671.tif 000005
【００１５】
【数４】
2018136671.tif 000006
【００１６】
ここで、ｊ＊はモード帯域を表すパラメータであり、Ｅｎは内部波のエネルギーを表すパラメータである。ｆはコリオリパラメータ（慣性周波数）であり、ｆ＝２Ωｓｉｎ（φ）によって計算される。Ωは地球の自転角速度であり、φは緯度である。ｂおよびＮ０は、浮力周波数Ｎ（ｚ）が下記式（５）で表されるときの定数である。Ｎ０は海面での浮力周波数であり、ｂはＮ（ｚ）の深度に対する低下度合いを決めるパラメータである。なお、深度とは海面からの深さのことである。
【００１７】
【数５】
2018136671.tif 000007
【００１８】
任意のＮ（ｚ）については、式（５）の関係から、水深Ｚが十分大きいと仮定すると、下記式（６）の関係が得られるため、本実施の形態では、式（１）を下記式（７）のように変形して使用する。
【００１９】
【数６】
2018136671.tif 000008
【００２０】
【数７】
2018136671.tif 000009
【００２１】
式（７）の右辺の各要素のうち、Ｎ（ｚ）は現場において水温及び塩分濃度を測定することにより求められる。また、周波数ω、コリオリパラメータｆ、およびモード番号ｊは既知の量である。
【００２２】
従って、スペクトル特性Ｆ（ω，ｊ）を求めるためには、Ｎ（ｚ）の深度に対する低下度合いを決めるパラメータｂと内部波のエネルギーＥｎとを乗算したｂＥｎと、モード帯域ｊ＊とを知る必要がある。以下では、ｂＥｎを総エネルギー量Ｅ０と称する（Ｅ０＝ｂＥｎ）と共に、内部波のパラメータであるモード帯域ｊ＊および総エネルギー量Ｅ０を内部波パラメータともいう。
【００２３】
ここで、式（７）の右辺の既知量の項をＫ（ω）とし、式（７）の未知量を含む項をＸ（ｊ，ｊ＊，Ｅ０）として、式（７）を下記式（８）のように書き換える。以下では、Ｘ（ｊ，ｊ＊，Ｅ０）をモードスペクトルという。式（８）のＫ（ω）は下記式（９）によって定まり、式（８）のＸ（ｊ，ｊ＊，Ｅ０）は下記式（１０）によって定まる。
【００２４】
【数８】
2018136671.tif 000010
【００２５】
【数９】
2018136671.tif 000011
【００２６】
【数１０】
2018136671.tif 000012
【００２７】
ここで、モードスペクトルＸ（ｊ，ｊ＊，Ｅ０）は、以下のように推定することができる。内部波変位の周波数ωの成分が、深度ｚ１，ｚ２，・・・，ｚＬでそれぞれＰ個測定されたとし、各深度に対応するＰ個の成分を、それぞれ、ζｐ（ｚ１），ζｐ（ｚ２），・・・，ζｐ（ｚＬ）（ｐ＝１，２，…，Ｐ）とする。すると、異なる深度に対する内部波変位の積の期待値は、下記式（１１）のように近似することができる。測定深度数Ｌは、内部波変位を測定する深度の数である。変位数Ｐは、各深度に対応する内部波変位の個数である。
【００２８】
【数１１】
2018136671.tif 000013
【００２９】
ここで、Ｊは内部波変位に影響を与えているモード数であり、十分大きな数を与えることにより、式（１１）による近似の誤差を小さくすることができる。また、Ｗｊ（ｚ）は、下記式（１２）および式（１３）の微分方程式を満足する内部波モードであり、αｊはｊ番目の内部波モードＷｊ（ｚ）に対応する固有波数である。
【００３０】
【数１２】
2018136671.tif 000014
【００３１】
【数１３】
2018136671.tif 000015
【００３２】
ｍ＝１，２，・・・，Ｌ、およびｎ＝１，２，・・・，Ｌに対して式（１１）を連立すると、初期観測ベクトルをｙ、初期観測行列をＥ、未知のモードスペクトルのベクトルをｘとした下記式（１４）の連立方程式が得られる。式（１４）の初期観測ベクトルｙは下記式（１５）によって定まり、式（１４）の初期観測行列Ｅは下記式（１６）によって定まり、式（１４）のモードスペクトルのベクトルｘは下記式（１７）によって定まる。
【００３３】
【数１４】
2018136671.tif 000016
【００３４】
【数１５】
2018136671.tif 000017
【００３５】
【数１６】
2018136671.tif 000018
【００３６】
【数１７】
2018136671.tif 000019
【００３７】
式（１４）を一般的な線形逆問題の手法によって解くことにより、モードスペクトルのベクトルｘの値を得ることができる。
【００３８】
また、式（１０）で表されるモードスペクトルＸ（ｊ，ｊ＊，Ｅ０）をモード帯域ｊ＊および総エネルギー量Ｅ０の関数とみて、モード帯域の初期値ｊ＊（０）および総エネルギー量の初期値Ｅ０（０）の周りでテーラー展開し、１次の項までをとると、「ｊ＊＝ｊ＊（０）＋Δｊ＊」および「Ｅ０＝Ｅ０（０）＋ΔＥ０」におけるＸ（ｊ，ｊ＊，Ｅ０）の値は、下記式（１８）のように近似することができる。以降では、モード帯域ｊ＊の初期値ｊ＊（０）を単に「初期値ｊ＊（０）」ともいい、総エネルギー量Ｅ０の初期値Ｅ０（０）を単に「初期値Ｅ０（０）」ともいう。なお、初期値ｊ＊（０）は本発明の「第１初期値」に相当し、初期値Ｅ０（０）は本発明の「第２初期値」に相当する。
【００３９】
【数１８】
2018136671.tif 000020
【００４０】
ここで、ｊ＊による偏微分は下記式（１９）のようになる。式（１９）のＨ’０（ｊ＊）は下記式（２０）で表され、式（２０）の最後の無限級数は下記式（２１）の公式を使って計算することができる。一方、Ｅ０による偏微分は、下記式（２２）のようになる。
【００４１】
【数１９】
2018136671.tif 000021
【００４２】
【数２０】
2018136671.tif 000022
【００４３】
【数２１】
2018136671.tif 000023
【００４４】
【数２２】
2018136671.tif 000024
【００４５】
ここで、式（１８）を式（１７）に代入すると、式（１４）の連立方程式は、下記式（２３）のように変形することができる。式（２３）のｘ０は下記式（２４）によって定まる内部波パラメータの初期ベクトルであり、式（２３）のＤは下記式（２５）によって定まる微分ベクトルである。
【００４６】
【数２３】
2018136671.tif 000025
【００４７】
【数２４】
2018136671.tif 000026
【００４８】
【数２５】
2018136671.tif 000027
【００４９】
ここで、式（２３）を整理すると、下記式（２６）のように、観測ベクトルをΔｙ、観測行列をＭ、未知の内部波パラメータの差分ベクトルをΔｘとする新たな連立方程式が得られる。観測ベクトルΔｙは下記式（２７）によって定まり、観測行列Ｍは下記式（２８）によって定まり、差分ベクトルΔｘは下記式（２９）によって定まる。
【００５０】
【数２６】
2018136671.tif 000028
【００５１】
【数２７】
2018136671.tif 000029
【００５２】
【数２８】
2018136671.tif 000030
【００５３】
【数２９】
2018136671.tif 000031
【００５４】
式（２６）を一般的な連立方程式の解法によって解き、内部波パラメータの差分ベクトルΔｘ＝［Δｊ＊  ΔＥ０］Ｔ（Ｔは行列の転置を表す記号である）を得ることができたら、モード帯域ｊ＊は下記式（３０）によって求めることができ、総エネルギー量Ｅ０は下記式（３１）によって求めることができる。
【００５５】
【数３０】
2018136671.tif 000032
【００５６】
【数３１】
2018136671.tif 000033
【００５７】
なお、初期値ｊ＊（０）および初期値Ｅ０（０）と、真のモード帯域および真の総エネルギー量との差が大きいほど、式（１８）の線形近似による誤差が大きくなる。初期値ｊ＊（０）および初期値Ｅ０（０）は、例えば、下記式（３２）および式（３３）に示す典型値に設定することができる。式（３２）および式（３３）に示す初期値は、真値と大差ないため、線形近似の誤差を小さくすることができる。
【００５８】
【数３２】
2018136671.tif 000034
【００５９】
【数３３】
2018136671.tif 000035
【００６０】
ところで、式（３０）および式（３１）により求めたモード帯域ｊ＊および総エネルギー量Ｅ０は、式（３２）および式（３３）に示す初期値よりも真値に近づいていると考えられる。そのため、式（３０）により求めたモード帯域ｊ＊を新たな初期値ｊ＊（０）とし、式（３１）により求めた総エネルギー量Ｅ０を新たな初期値Ｅ０（０）として、式（２６）の連立方程式を構成するようにしてもよい。このように再構築した式（２６）の連立方程式を解くようにすれば、線形近似の誤差をさらに小さくすることができる。
【００６１】
続いて、図１を参照して、内部波変位シミュレーションシステム２００の構成について説明する。図１に示すように、内部波変位シミュレーションシステム２００は、内部波変位測定装置１１０と、内部波パラメータ推定装置１００と、スペクトル特性演算装置１２０と、内部波変位予測処理装置１３０と、を有している。
【００６２】
内部波変位測定装置１１０は、複数の深度ｚ１，ｚ２，・・・，ｚＬのそれぞれについて、Ｐ個の内部波変位であるζｐ（ｚ１），ζｐ（ｚ２），・・・，ζｐ（ｚＬ）（ｐ＝１，２，…，Ｐ）を測定するものである。本実施の形態において、測定深度数Ｌは２以上の任意の自然数であり、変位数Ｐは２以上の任意の自然数であるものとする。すなわち、内部波変位測定装置１１０は、深度の異なる複数の地点において、複数の内部波変位を求めるものである。内部波変位測定装置１１０が求める内部波変位は、内部波変位の周波数ωの成分であり、周波数ωは予め設定されている。内部波変位測定装置１１０は、内部波変位の測定深度である深度ｚ１，ｚ２，・・・，ｚＬと、測定した内部波変位ζｐ（ｚ１），ζｐ（ｚ２），・・・，ζｐ（ｚＬ）と、内部波変位の測定に係る周波数ωとを、内部波パラメータ推定装置１００へ出力するものである。
【００６３】
内部波パラメータ推定装置１００は、内部波変位測定装置１１０から入力された内部波変位の測定深度、複数の内部波変位、および周波数ωなどをもとに、内部波のパラメータとして、モード帯域ｊ＊および総エネルギー量Ｅ０を求めるものである。また、内部波パラメータ推定装置１００は、求めたモード帯域ｊ＊および総エネルギー量Ｅ０をスペクトル特性演算装置１２０へ出力するものである。
【００６４】
スペクトル特性演算装置１２０は、内部波パラメータ推定装置１００から出力されたモード帯域ｊ＊および総エネルギー量Ｅ０を用いて、式（８）により、スペクトル特性Ｆ（ω，ｊ）を求めるものである。また、スペクトル特性演算装置１２０は、求めたスペクトル特性Ｆ（ω，ｊ）を内部波変位予測処理装置１３０へ出力するものである。
【００６５】
内部波変位予測処理装置１３０は、スペクトル特性演算装置１２０から出力されたスペクトル特性Ｆ（ω，ｊ）をもとに、内部波変位をシミュレーションして予測するものである。また、内部波変位予測処理装置１３０は、内部波変位の予測の結果をもとに、海中の音波伝搬のシミュレーションなどを行うものである。
【００６６】
図２は、本発明の実施の形態に係る内部波パラメータ推定装置の機能的構成を示すブロック図である。図２に示すように、内部波パラメータ推定装置１００は、初期観測ベクトル生成部１０と、初期観測行列生成部２０と、観測ベクトル生成部３０と、観測行列生成部４０と、連立方程式求解部５０と、内部波パラメータ計算部６０と、記憶部７０と、を有している。
【００６７】
記憶部７０には、浮力周波数Ｎ（ｚ）およびコリオリパラメータｆが予め記憶されている。また、記憶部７０には、モード帯域の初期値ｊ＊（０）と総エネルギー量の初期値Ｅ０（０）が予め格納されている。
【００６８】
初期観測ベクトル生成部１０には、内部波変位測定装置１１０から、深度ｚ１，ｚ２，・・・，ｚＬのそれぞれにおいて測定された内部波変位ζｐ（ｚ１），ζｐ（ｚ２），・・・，ζｐ（ｚＬ）（ｐ＝１，２，…，Ｐ）が入力される。
【００６９】
初期観測ベクトル生成部１０は、内部波変位ζｐ（ｚ１），ζｐ（ｚ２），・・・，ζｐ（ｚＬ）を用いて、式（１１）および式（１５）により、初期観測ベクトルｙを求めるものである。また、初期観測ベクトル生成部１０は、求めた初期観測ベクトルｙを観測ベクトル生成部３０へ出力するものである。
【００７０】
初期観測行列生成部２０には、内部波変位測定装置１１０から、内部波変位の測定深度である深度ｚ１，ｚ２，・・・，ｚＬと、周波数ωとが入力される。また、初期観測行列生成部２０は、コリオリパラメータｆ、および浮力周波数Ｎ（ｚ）を記憶部７０から取得するものである。初期観測行列生成部２０は、入力された周波数ωと、記憶部７０から取得したコリオリパラメータｆおよび浮力周波数Ｎ（ｚ）とを用いて、式（２）および式（９）により、式（７）の右辺の既知量の項Ｋ（ω）を計算するものである。
【００７１】
また、初期観測行列生成部２０は、内部波変位の測定深度、周波数ω、コリオリパラメータｆ、および浮力周波数Ｎ（ｚ）を用いて、式（１２）および式（１３）を満足する内部波モードＷｊ（ｚｍ）（ｊ＝１，２，・・・，Ｊ、ｍ＝１，２，・・・，Ｌ）を求めるものである。初期観測行列生成部２０は、一般的な微分方程式の解法により、内部波モードＷｊ（ｚｍ）を求めるように構成されている。
【００７２】
さらに、初期観測行列生成部２０は、既知量の項Ｋ（ω）と内部波モードＷｊ（ｚｍ）（ｊ＝１，２，・・・，Ｊ、ｍ＝１，２，・・・，Ｌ）とを用いて、式（１６）により、初期観測行列Ｅを求めるものである。そして、初期観測行列生成部２０は、求めた初期観測行列Ｅを、観測ベクトル生成部３０および観測行列生成部４０に出力するものである。
【００７３】
観測ベクトル生成部３０には、初期観測ベクトル生成部１０から初期観測ベクトルｙが入力され、初期観測行列生成部２０から初期観測行列Ｅが入力される。また、観測ベクトル生成部３０は、モード帯域の初期値ｊ＊（０）および総エネルギー量の初期値Ｅ０（０）を記憶部７０から取得するものである。
【００７４】
観測ベクトル生成部３０は、記憶部７０から取得したモード帯域の初期値ｊ＊（０）および総エネルギー量の初期値Ｅ０（０）を用いて、式（３）、式（４）、式（１０）、および式（２４）により、内部波パラメータの初期ベクトルｘ０を求めるものである。また、観測ベクトル生成部３０は、求めた初期ベクトルｘ０と、入力された初期観測ベクトルｙおよび初期観測行列Ｅとを用いて、式（２７）により、観測ベクトルΔｙを求めるものである。そして、観測ベクトル生成部３０は、求めた観測ベクトルΔｙを連立方程式求解部５０に出力するものである。
【００７５】
観測行列生成部４０には、初期観測行列生成部２０から初期観測行列Ｅが入力される。また、観測行列生成部４０は、モード帯域の初期値ｊ＊（０）、および総エネルギー量の初期値Ｅ０（０）を記憶部７０から取得するものである。
【００７６】
観測行列生成部４０は、記憶部７０から取得したモード帯域の初期値ｊ＊（０）、および総エネルギー量の初期値Ｅ０（０）を用いて、式（４）、式（１９）、式（２０）、式（２１）、式（２２）、および式（２５）により、微分ベクトルＤを求めるものである。また、観測行列生成部４０は、求めた微分ベクトルＤと、入力された初期観測行列Ｅとを用いて、式（２８）により、観測行列Ｍを求めるものである。そして、観測行列生成部４０は、求めた観測行列Ｍを連立方程式求解部５０に出力するものである。
【００７７】
連立方程式求解部５０には、観測ベクトル生成部３０から観測ベクトルΔｙが入力され、観測行列生成部４０から観測行列Ｍが入力される。連立方程式求解部５０は、入力された観測ベクトルΔｙおよび観測行列Ｍを用いて、式（２６）により、内部波パラメータの差分ベクトルΔｘを求めるものである。連立方程式求解部５０は、ＱＲ法などの一般的な連立方程式の解法によって式（２６）を解き、内部波パラメータの差分ベクトルΔｘを求めるように構成されている。また、連立方程式求解部５０は、求めた内部波パラメータの差分ベクトルΔｘを、内部波パラメータ計算部６０に出力するものである。
【００７８】
内部波パラメータ計算部６０には、連立方程式求解部５０から内部波パラメータの差分ベクトルΔｘが入力される。また、内部波パラメータ計算部６０は、モード帯域の初期値ｊ＊（０）および総エネルギー量の初期値Ｅ０（０）を記憶部７０から取得するものである。内部波パラメータ計算部６０は、入力された内部波パラメータの差分ベクトルΔｘと、取得したモード帯域の初期値ｊ＊（０）および総エネルギー量の初期値Ｅ０（０）とを用いて、式（３０）および式（３１）により、モード帯域ｊ＊および総エネルギー量Ｅ０を内部波パラメータとして求めるものである。そして、内部波パラメータ計算部６０は、求めた内部波パラメータをスペクトル特性演算装置１２０に出力する。
【００７９】
ここで、内部波パラメータ計算部６０は、求めた内部波パラメータによって、記憶部７０内のモード帯域の初期値ｊ＊（０）および総エネルギー量の初期値Ｅ０（０）を更新するようにしてもよい。すなわち、内部波パラメータ計算部６０は、式（３０）により推定したモード帯域ｊ＊を初期値ｊ＊（０）に設定すると共に、式（３１）により推定した総エネルギー量Ｅ０を初期値Ｅ０（０）に設定するようにしてもよい。このように、予め設定された初期値よりも真値に近いと考えられるモード帯域ｊ＊および総エネルギー量Ｅ０によって初期値ｊ＊（０）および初期値Ｅ０（０）を更新すれば、線形近似の誤差をさらに小さくすることができる。
【００８０】
ところで、上述した原理の説明では、初期値ｊ＊（０）および初期値Ｅ０（０）として、式（３２）および（３３）の値を与える場合を例示したが、これに限らず、他に真値に近いと考えられる値があれば、その値を初期値として用いてもよい。
【００８１】
ここで、初期観測ベクトル生成部１０、初期観測行列生成部２０、観測ベクトル生成部３０、観測行列生成部４０、連立方程式求解部５０、および内部波パラメータ計算部６０は、それぞれの各機能を実現する回路デバイスなどのハードウェアで実現することもできるし、例えば、マイコン、ＤＳＰ（Ｄｉｇｉｔａｌ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｏｒ）、又はＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）等の演算装置上で実行されるソフトウェアとして実現することもできる。後者の場合、記憶部７０には、初期観測ベクトル生成部１０、初期観測行列生成部２０、観測ベクトル生成部３０、観測行列生成部４０、連立方程式求解部５０、および内部波パラメータ計算部６０の動作プログラムが格納されていてもよい。また、記憶部７０は、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）、フラッシュメモリ等のＰＲＯＭ（Ｐｒｏｇｒａｍｍａｂｌｅ  ＲＯＭ）、又はＨＤＤ（Ｈａｒｄ  Ｄｉｓｋ  Ｄｒｉｖｅ）等により構成することができる。
【００８２】
図３は、内部波パラメータ推定装置の動作を示すフローチャートである。図３に基づき、内部波パラメータ推定装置１００による内部波パラメータ推定方法について説明する。
【００８３】
初期観測ベクトル生成部１０は、海中における複数の深度それぞれの内部波変位をもとに初期観測ベクトルｙを生成し、生成した初期観測ベクトルｙを観測ベクトル生成部３０へ出力する（ステップＳ１０１／初期観測ベクトル生成工程）。
【００８４】
また、初期観測行列生成部２０は、内部波変位の測定に係る周波数ωと、コリオリパラメータｆと、浮力周波数Ｎ（ｚ）とを用いて、式（７）の右辺の既知量の項Ｋ（ω）を計算する。次いで、初期観測行列生成部２０は、内部波変位の測定深度、周波数ω、コリオリパラメータｆ、および浮力周波数Ｎ（ｚ）を用いて、内部波モードＷｊ（ｚｍ）を求める。そして、初期観測行列生成部２０は、既知量の項Ｋ（ω）と内部波モードＷｊ（ｚｍ）とをもとに初期観測行列Ｅを生成し、生成した初期観測行列Ｅを観測ベクトル生成部３０および観測行列生成部４０に出力する（ステップＳ１０２／初期観測行列生成工程）。
【００８５】
次に、観測ベクトル生成部３０は、初期値ｊ＊（０）および初期値Ｅ０（０）を用いて内部波パラメータの初期ベクトルｘ０を求める。そして、観測ベクトル生成部３０は、求めた初期ベクトルｘ０と、初期観測ベクトルｙおよび初期観測行列Ｅとをもとに観測ベクトルΔｙを求めて連立方程式求解部５０に出力する（ステップＳ１０３／観測ベクトル生成工程）。
【００８６】
また、観測行列生成部４０は、初期値ｊ＊（０）および初期値Ｅ０（０）を用いて微分ベクトルＤを求める。そして、観測行列生成部４０は、求めた微分ベクトルＤと、初期観測行列Ｅとをもとに観測行列Ｍを求めて連立方程式求解部５０に出力する（ステップＳ１０４／観測行列生成工程）。
【００８７】
連立方程式求解部５０は、観測ベクトルΔｙおよび観測行列Ｍを用いて内部波パラメータの差分ベクトルΔｘを求め、求めた差分ベクトルΔｘを内部波パラメータ計算部６０に出力する（ステップＳ１０５／連立方程式求解工程）。
【００８８】
続いて、内部波パラメータ計算部６０は、内部波パラメータの差分ベクトルΔｘと、初期値ｊ＊（０）および初期値Ｅ０（０）とを用いてモード帯域ｊ＊および総エネルギー量Ｅ０を求めるものである。そして、内部波パラメータ計算部６０は、求めたモード帯域ｊ＊および総エネルギー量Ｅ０を内部波パラメータとしてスペクトル特性演算装置１２０に出力する（ステップＳ１０６／内部波パラメータ計算工程）。
【００８９】
また、内部波パラメータ計算部６０は、内部波パラメータ計算工程で求めたモード帯域ｊ＊および総エネルギー量Ｅ０によって、記憶部７０内の初期値ｊ＊（０）および初期値Ｅ０（０）を更新する（ステップＳ１０７／初期値更新工程）。そして、内部波パラメータ推定装置１００は、上記ステップＳ１０３〜Ｓ１０７の一連の処理を繰り返し実行する。
【００９０】
ここで、上記の動作説明は、図３に付したステップ番号の順に説明したが、これに限定されるものではない。例えば、ステップＳ１０１の初期観測ベクトル生成工程とステップＳ１０２の初期観測行列生成工程とは並行して行うことができ、ステップＳ１０３の観測ベクトル生成工程とステップＳ１０４の観測行列生成工程とも並行して行うことができる。
【００９１】
以上のように、本実施の形態における内部波パラメータ推定装置１００は、初期観測ベクトルｙおよび初期観測行列Ｅを用いて生成した観測ベクトルΔｙおよび観測行列Ｍをもとにモード帯域ｊ＊および総エネルギー量Ｅ０を求める。そのため、モード数Ｊを大きくしても演算量が変化せず、かつ最小自乗法などによる演算処理が不要となることから、内部波のパラメータを推定する際の演算処理を簡略化することができる。したがって、内部波変位シミュレーションシステム２００によれば、少ない演算量で内部波変位を精度よくシミュレーションすることができる。
【００９２】
また、内部波パラメータ計算部６０には、求めたモード帯域ｊ＊および総エネルギー量Ｅ０によって、初期値ｊ＊（０）および初期値Ｅ０（０）を更新する機能をもたせることができる。このようにすれば、設定された初期値をさらに真値に近づけることができるため、線形近似の誤差をより小さくすることができる。
【００９３】
さらに、初期観測ベクトル生成部１０は、海中における複数の深度それぞれの内部波変位をもとに初期観測ベクトルｙを生成するものであり、初期観測行列生成部２０は、複数の深度を用いて初期観測行列Ｅを生成するものである。すなわち、本実施の形態では、測定深度数Ｌが２以上の自然数に設定されている。したがって、より広域の内部波変位を用いることができることから、内部波パラメータの推定精度が向上するため、内部波変位のシミュレーションの精度向上を図ることができる。
【００９４】
ここで、図４および図５を参照して、内部波パラメータ推定装置１００によって得られる効果を詳細に説明する。図４は、従来技術に係る内部波パラメータ推定装置の構成を示すブロック図である。図５は、図２の内部波パラメータ推定装置による演算回数と、図４の内部波パラメータ推定装置による演算回数とを示すグラフである。図５では、横軸にモード数Ｊをとり、縦軸に浮動小数点演算回数をとっている。なお、図５には、内部波変位の測定深度数Ｌが３０の場合を例示している。
【００９５】
図４に示すように、従来の構成である内部波パラメータ推定装置１０００は、初期観測ベクトルｙを生成する初期観測ベクトル生成部１１００と、初期観測行列Ｅを生成する初期観測行列生成部１２００と、初期観測ベクトルｙと初期観測行列Ｅとからモードスペクトルのベクトルｘを求める連立方程式求解部１５００と、を有している。
【００９６】
すなわち、内部波パラメータ推定装置１０００の連立方程式求解部１５００には、初期観測ベクトルｙおよび初期観測行列Ｅが入力されるようになっている。そして、連立方程式求解部１５００は、初期観測ベクトルｙおよび初期観測行列Ｅをもとに、ＱＲ法などの一般的な連立方程式の解法によって式（１４）を解き、モードスペクトルのベクトルｘを出力するようになっている。そのため、モード帯域ｊ＊および総エネルギー量Ｅ０を推定するためには、モードスペクトルのベクトルｘが得られた後に、最小自乗法などの演算処理を行う必要がある。
【００９７】
ここで、例えば、下記式（３４）のように誤差関数を定義する。ｘｊは、求めたモードスペクトルのベクトルｘのｊ番目の要素である。
【００９８】
【数３４】
2018136671.tif 000036
【００９９】
つまり、従来の構成では、式（３４）の誤差関数が最小となるように、モード帯域ｊ＊および総エネルギー量Ｅ０を推定する必要がある。
【０１００】
また、式（１４）および式（１６）の初期観測行列Ｅは、行数がＬ（Ｌ＋１）／２であり、列数がＪである。したがって、式（１１）における近似の誤差を小さくするためにモード数Ｊを大きくすると、モードスペクトルのベクトルｘを求めるための演算量が大きくなる。
【０１０１】
ここで、例えば、連立方程式を解く手法として広く使われているＱＲ法を用いる場合を考える。ＱＲ法は、初期観測行列Ｅが「Ｅ＝ＱＲ」のように対角行列Ｑと上三角行列Ｒとの積で表せることを利用するものであり、浮動小数点演算が約２Ｊ２［Ｌ（Ｌ＋１）／２−Ｊ／３］回必要であることが知られている。
【０１０２】
図５のように、測定深度数Ｌが３０の場合、従来の構成では、浮動小数点演算が約２Ｊ２［４６５−Ｊ／３］回必要となる。すなわち、従来の構成では、図５の破線グラフＧｙのように、モード数Ｊを大きくするにつれて演算量が増加していくことがわかる。
【０１０３】
一方、本実施の形態では、式（２６）の連立方程式において、観測行列Ｍの列数は常に２である。よって、式（２６）の連立方程式をＱＲ法で解く場合には、浮動小数点演算の回数は約８×［Ｌ（Ｌ＋１）／２−２／３］回となる。
【０１０４】
ここで、線形近似の誤差は、浮動小数点演算を繰り返し行う毎に小さくなり、数回の繰り返しによって十分に小さくなる。したがって、本実施の形態では、浮動小数点演算を繰り返す回数である繰り返し回数Ｓが２以上に設定されている。この場合、浮動小数点演算の回数は約８×Ｓ×［Ｌ（Ｌ＋１）／２−２／３］回となる。図５では、繰り返し回数Ｓが１０に設定されている場合を例示しており、この場合、浮動小数点演算の回数は約８０×［Ｌ（Ｌ＋１）／２−２／３］回となる。
【０１０５】
すなわち、従来の構成では、図５の破線グラフＧｙのように、モード数Ｊを大きくするほど、浮動小数点演算の回数が増加していたが、内部波パラメータ推定装置１００では、図５の実線グラフＧｘのように、浮動小数点演算の回数がモード数Ｊに依存しない。つまり、内部波パラメータ推定装置１００においては、モード数Ｊを大きくしても、浮動小数点演算の回数が多くならない。
【０１０６】
したがって、式（１１）における近似の誤差を小さくするために、モード数Ｊを大きくすることを考えるとき、例えば、繰り返し回数Ｓを１０に設定した場合は、モード数を７以上にすれば、従来の構成よりも少ない演算回数で内部波パラメータを推定することができる。なお、非特許文献３では、モード数Ｊを１５としており、かかる設定の場合、図５に演算回数差ΔＣとして示すように、内部波パラメータ推定装置１００と従来の構成との間に、演算回数の大きな違いが生じる。つまり、内部波パラメータ推定装置１００によれば、モード数Ｊを大きくすればするほど、連立方程式を解くための演算回数が従来よりも削減され、優位性が向上する。
【０１０７】
また、内部波パラメータ推定装置１００は、式（２６）の連立方程式「Δｙ＝ＭΔｘ」を解くことにより、内部波パラメータを得ることができるため、従来の構成のように、最小自乗法などによる内部波パラメータの推定処理を行う必要がない。すなわち、内部波パラメータ推定装置１００によれば、内部波のパラメータを推定する際の演算処理を簡略化することができる。
【０１０８】
なお、上述した実施の形態は、内部波パラメータ推定装置および内部波パラメータ推定方法における好適な具体例であり、本発明の技術的範囲は、これらの態様に限定されるものではない。例えば、上記実施の形態では、測定深度数Ｌが２以上の自然数である場合を例示したが、これに限らず、測定深度数Ｌは１であってもよい。ただし、測定深度数Ｌは、大きな値に設定した方が、より広域の内部波変位を用いることができるため、内部波パラメータの推定精度が向上することから、内部波変位のシミュレーションの精度向上を図ることができる。また、上記実施の形態では、変位数Ｐが２以上の自然数である場合を例示したが、これに限らず、変位数Ｐは１であってもよい。ただし、変位数Ｐについても、大きな値に設定した方が、より広域の内部波変位を用いることができるため、内部波パラメータの推定精度が向上することから、内部波変位のシミュレーションの精度向上を図ることができる。
【符号の説明】
【０１０９】
１０  初期観測ベクトル生成部、２０  初期観測行列生成部、３０  観測ベクトル生成部、４０  観測行列生成部、５０  連立方程式求解部、６０  内部波パラメータ計算部、７０  記憶部、１００  内部波パラメータ推定装置、１１０  内部波変位測定装置、１２０  スペクトル特性演算装置、１３０  内部波変位予測処理装置、２００  内部波変位シミュレーションシステム。
【図１】
2018136671.tif 000037
【図２】
2018136671.tif 000038
【図３】
2018136671.tif 000039
【図４】
2018136671.tif 000040
【図５】
2018136671.tif 000041