JP2009075017A - 表層流速推定方法、装置並びにプログラム - Google Patents

Abstract

【課題】海洋レーダによって得られるデータを用いた沿岸表層流動の観測において、外因性ノイズの影響を排除してドップラースペクトルの正当なピーク周波数を短時間で検出することができるようにする。
【解決手段】ドップラースペクトルの正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出すると共にドップラースペクトルの負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出し、部分領域毎の正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値との合算値が最大である部分領域と連接し且つ部分領域毎の合算値が所定の閾値よりも大きい部分領域からなる連接領域の中からピーク周波数を検出するようにした。
【選択図】図１

Description

本発明は、表層流速推定方法、装置並びにプログラムに関する。さらに詳述すると、本発明は、海洋レーダによって得られるデータを用いた沿岸表層流動の観測に用いて好適な表層流速推定方法、装置並びにプログラムに関する。

海洋レーダは、陸上からのリモートセンシングにより広域の流速、波高、風向といった表層流動を天候に関係なく長期間連続観測できる装置であり、沿岸環境モニタリングや港湾管理、沿岸防災、漁業情報の提供など多分野への活用が進められている。

海洋レーダによる観測は、送信アンテナから観測領域に向けてレーダ波を照射し、海面でブラッグ散乱されて戻ってくる電波である海洋波からのエコーを受信アンテナで受信し、エコー強度の時系列からドップラースペクトルを求めて表層流動を抽出する。

受信アンテナで受信されるエコーはブラッグ共鳴条件を満たす波浪成分波による後方散乱エコーであり、このエコーのドップラースペクトルは第一次散乱及び第二次散乱のスペクトルによって構成される。表層流動の流速は第一次散乱スペクトルを解析することによって得られ、波浪情報は第一次散乱スペクトル及び第二次散乱スペクトルを解析することによって得られる。

後方散乱エコーのドップラースペクトルから流速を求める際の第一次散乱スペクトル周波数の検出のアルゴリズムは、通常、ドップラースペクトル分布の正の周波数領域と負の周波数領域とにおいて最も大きなスペクトル値を持つピークを検出してそれらのピーク周波数を得るようにしている（特許文献１）。

特開平１１−８３９９２号

しかしながら、レーダ波照射中に例えばレーダ波をほぼ直角に横切る大きな船舶が存在すると大きなエコーを受信してしまう等のように、受信アンテナで受信されるエコーは外因性のノイズ及び受信機のノイズの影響を受けることは避けられない。そして、このような外因性ノイズが混入するとドップラースペクトルのピークが不明瞭になったり外因性ノイズによるピークが顕著になったりするので、流速に対応する正当な第一次散乱スペクトル周波数を検出することができなくなり、表層流動の観測精度が低下してしまうという問題がある。

また、突発的な事象を対象とした沿岸環境モニタリングや港湾管理、沿岸防災においては特にリアルタイム若しくはほぼリアルタイムでの情報提供が必要とされるが、従来の方法はドップラースペクトルのピーク周波数の検出に多大な時間を必要とするので真に必要とされるタイミングで情報を提供することができるとは言い難い。

そこで、本発明は、外因性ノイズの影響を排除してドップラースペクトルの正当なピーク周波数を短時間で検出することができる表層流速推定方法、装置並びにプログラムを提供することを目的とする。

かかる目的を達成するため、請求項１記載の表層流速推定方法は、ビーム方向距離別・周波数別ドップラースペクトルのピーク周波数を検出して表層流速を推定する方法において、ドップラースペクトルの正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出すると共にドップラースペクトルの負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出し、部分領域毎の正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値との合算値が最大である部分領域と連接し且つ部分領域毎の合算値が所定の閾値よりも大きい部分領域からなる連接領域の中からピーク周波数を検出するようにしている。

したがって、この表層流速推定方法によると、誤差情報の検定処理を行って元のドップラースペクトルデータから一定の条件を満たす領域を抽出し、当該検定処理が施された領域の中からピーク周波数を検出するようにしているので、自然の雑音や人工雑音等の外因性ノイズが含まれる実際の観測データからノイズを適切に除去して適確な流動情報が検出される。

また、請求項２記載の表層流速推定方法は、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行うステップと、ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離するステップと、正の領域のスペクトル及び負の領域のスペクトルについて周波数方向の平均値を算出するステップと、正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出するステップと、正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に負の領域のスペクトルの偏差について部分領域の平均値を算出するステップと、正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値とを部分領域毎に合算して部分領域毎の合算値を算出するステップと、部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成するステップと、二値化行列について隣接行列を作成するステップと、隣接行列を定常状態になるまで掛け合わせて定常解行列を作成するステップと、定常解行列から部分領域毎の合算値が最大である部分領域と連接し且つ部分領域毎の合算値が二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して二値化行列からフィルタ基礎行列を作成するステップと、フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成するステップと、正の領域のスペクトルの偏差と負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出するステップと、フィルタ行列をスペクトル偏差合算値に適用してフィルタ適用スペクトルを作成するステップと、フィルタ適用スペクトルからピーク周波数を検出するステップと、ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出するステップとを有するようにしている。

また、請求項３記載の表層流速推定装置は、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行う手段と、ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離する手段と、正の領域のスペクトル及び負の領域のスペクトルについて周波数方向の平均値を算出する手段と、正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出する手段と、正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に負の領域のスペクトルの偏差について部分領域の平均値を算出する手段と、正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値とを部分領域毎に合算して部分領域毎の合算値を算出する手段と、部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成する手段と、二値化行列について隣接行列を作成する手段と、隣接行列を定常状態になるまで掛け合わせて定常解行列を作成する手段と、定常解行列から部分領域毎の合算値が最大である部分領域と連接し且つ部分領域毎の合算値が二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して二値化行列からフィルタ基礎行列を作成する手段と、フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成する手段と、正の領域のスペクトルの偏差と負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出する手段と、フィルタ行列をスペクトル偏差合算値に適用してフィルタ適用スペクトルを作成する手段と、フィルタ適用スペクトルからピーク周波数を検出する手段と、ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出する手段とを有するようにしている。

また、請求項４記載の表層流速推定プログラムは、ビーム方向距離別・周波数別ドップラースペクトルのピーク周波数を検出して表層流速の推定を行う際に、少なくとも、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行う手段、ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離する手段、正の領域のスペクトル及び負の領域のスペクトルについて周波数方向の平均値を算出する手段、正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出する手段、正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に負の領域のスペクトルの偏差について部分領域の平均値を算出する手段、正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値とを部分領域毎に合算して部分領域毎の合算値を算出する手段、部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成する手段、二値化行列について隣接行列を作成する手段、隣接行列を定常状態になるまで掛け合わせて定常解行列を作成する手段、定常解行列から部分領域毎の合算値が最大である部分領域と連接し且つ部分領域毎の合算値が二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して二値化行列からフィルタ基礎行列を作成する手段、フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成する手段、正の領域のスペクトルの偏差と負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出する手段、フィルタ行列をスペクトル偏差合算値に適用してフィルタ適用スペクトルを作成する手段、フィルタ適用スペクトルからピーク周波数を検出する手段、ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出する手段としてコンピュータを機能させるようにしている。

したがって、この表層流速推定方法、装置並びにプログラムによると、誤差情報の検定処理を行って元のドップラースペクトルデータから一定の条件を満たす領域を抽出し、当該検定処理が施された領域の中から検出されたピーク周波数を用いて表層流速を推定するようにしているので、自然の雑音や人工雑音等の外因性ノイズが含まれる実際の観測データからノイズを適切に除去して検出された適確な流動情報に基づく表層流速が推定される。

なお、本発明においては、ドップラースペクトルデータとは、流速情報を含む周波数方向のドップラースペクトルの値のデータという意味で用いる。

本発明の表層流速推定方法、装置並びにプログラムによれば、自然の雑音や人工雑音等の外因性ノイズが含まれる実際の観測データからノイズを適切に除去して適確な流動情報具体的にはドップラースペクトルの正当なピーク周波数を検出することができると共に、当該適確な流動情報に基づいて表層流速が推定されるので、推定精度の向上が可能になる。また、適確な流動情報を検出することができるので、例えばピーク周波数の検出結果のチェック等の煩雑な後処理が必要とされず推定に係る手間を省くことが可能になると共に推定を短時間で行うことが可能になる。さらにこれにより、本発明によれば、リアルタイム海域流動モニタリングに資する流動情報を提供することができるので、突発的な事象を対象とした沿岸環境モニタリングや港湾管理、沿岸防災に対して有用な情報提供を行うことが可能になる。

以下、本発明の構成を図面に示す最良の形態に基づいて詳細に説明する。

図１から図１２に、本発明の表層流速推定方法、装置並びにプログラムの実施形態の一例を示す。この表層流速推定方法は、ビーム方向距離別・周波数別ドップラースペクトルのピーク周波数を検出して表層流速を推定する方法において、ドップラースペクトルの正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出すると共にドップラースペクトルの負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出し、部分領域毎の正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値との合算値が最大である部分領域と連接し且つ部分領域毎の合算値が所定の閾値よりも大きい部分領域からなる連接領域の中からピーク周波数を検出するようにしている。

そして、上記表層流速推定方法は、図１に示すように、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行うステップ（Ｓ１）と、ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離するステップ（Ｓ２）と、正の領域のスペクトル及び負の領域のスペクトルについて周波数方向の平均値を算出するステップ（Ｓ３）と、正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出するステップ（Ｓ４）と、正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に負の領域のスペクトルの偏差について部分領域の平均値を算出するステップ（Ｓ５）と、正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値とを部分領域毎に合算して部分領域毎の合算値を算出するステップ（Ｓ６）と、部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成するステップ（Ｓ７）と、二値化行列について隣接行列を作成するステップ（Ｓ８）と、隣接行列を定常状態になるまで掛け合わせて定常解行列を作成するステップ（Ｓ９）と、定常解行列から部分領域毎の合算値が最大である部分領域と連接し且つ部分領域毎の合算値が二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して二値化行列からフィルタ基礎行列を作成するステップ（Ｓ１０）と、フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成するステップ（Ｓ１１）と、正の領域のスペクトルの偏差と負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出するステップ（Ｓ１２）と、フィルタ行列をスペクトル偏差合算値に適用してフィルタ適用スペクトルを作成するステップ（Ｓ１３）と、フィルタ適用スペクトルからピーク周波数を検出するステップ（Ｓ１４）と、ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出するステップ（Ｓ１５）とからなる処理構成によって実現される。

また、上記表層流速推定方法は、本発明の表層流速推定装置として実現される。本発明の表層流速推定装置は、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行う手段と、ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離する手段と、正の領域のスペクトル及び負の領域のスペクトルについて周波数方向の平均値を算出する手段と、正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出する手段と、正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に負の領域のスペクトルの偏差について部分領域の平均値を算出する手段と、正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値とを部分領域毎に合算して部分領域毎の合算値を算出する手段と、部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成する手段と、二値化行列について隣接行列を作成する手段と、隣接行列を定常状態になるまで掛け合わせて定常解行列を作成する手段と、定常解行列から部分領域毎の合算値が最大である部分領域と連接し且つ部分領域毎の合算値が二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して二値化行列からフィルタ基礎行列を作成する手段と、フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成する手段と、正の領域のスペクトルの偏差と負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出する手段と、フィルタ行列をスペクトル偏差合算値に適用してフィルタ適用スペクトルを作成する手段と、フィルタ適用スペクトルからピーク周波数を検出する手段と、ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出する手段とを有する。

上述の表層流速推定方法並びに表層流速推定装置は、本発明の表層流速推定プログラムをコンピュータ上で実行することによっても実現される。本実施形態では、表層流速推定プログラムをコンピュータ上で実行する場合を例に挙げて説明する。

表層流速推定プログラム１７を実行するための本実施形態の表層流速推定装置１０の全体構成を図２に示す。この表層流速推定装置１０は、制御部１１、記憶部１２、入力部１３、表示部１４及びメモリ１５を備え相互にバス等の信号回線により接続されている。また、表層流速推定装置１０にはデータサーバ１６が通信回線等により接続されており、その通信回線等を介して相互にデータや制御指令等の信号の送受信（出入力）が行われる。

制御部１１は記憶部１２に記憶されている表層流速推定プログラム１７によって表層流速推定装置１０全体の制御並びに表層流速の推定に係る演算を行うものであり、例えばＣＰＵ（中央演算処理装置）である。記憶部１２は少なくともデータやプログラムを記憶可能な装置であり、例えばハードディスクである。メモリ１５は制御部１１が各種制御や演算を実行する際の作業領域であるメモリ空間となる。

入力部１３は少なくとも作業者の命令を制御部１１に与えるためのインターフェイスであり、例えばキーボードである。

表示部１４は制御部１１の制御により文字や図形等の表示を行うものであり、例えばディスプレイである。

また、データサーバ１６は少なくともデータを記憶可能なサーバである。

表層流速推定装置１０の制御部１１には、表層流速推定プログラム１７を実行することにより、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行うデータ読込部１１ａと、ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離するスペクトル分離部１１ｂと、正の領域のスペクトル及び負の領域のスペクトルについて周波数方向の平均値を算出するスペクトル平均算出部１１ｃと、正の領域のスペクトルから正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に負の領域のスペクトルから負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出するスペクトル偏差算出部１１ｄと、正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に負の領域のスペクトルの偏差について部分領域の平均値を算出する部分領域偏差平均値算出部１１ｅと、正の領域のスペクトルの偏差の部分領域の平均値と負の領域のスペクトルの偏差の部分領域の平均値とを部分領域毎に合算して部分領域毎の合算値を算出する部分領域偏差平均値合算部１１ｆと、部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成する二値化行列作成部１１ｇと、二値化行列について隣接行列を作成する隣接行列作成部１１ｈと、隣接行列を定常状態になるまで掛け合わせて定常解行列を作成する定常解行列作成部１１ｉと、定常解行列から部分領域毎の合算値が最大である部分領域と連接し且つ部分領域毎の合算値が二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して二値化行列からフィルタ基礎行列を作成するフィルタ基礎行列作成部１１ｊと、フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成するフィルタ行列作成部１１ｋと、正の領域のスペクトルの偏差と負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出するスペクトル偏差合算部１１ｌと、フィルタ行列をスペクトル偏差合算値に適用してフィルタ適用スペクトルを作成するフィルタ適用部１１ｍと、フィルタ適用スペクトルからピーク周波数を検出するピーク周波数検出部１１ｎと、ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出する流速算出部１１ｏとが構成される。

はじめに、第一段階として元のドップラースペクトルデータを用いて部分領域スペクトルデータを作成する処理（Ｓ１〜Ｓ５）について説明する。

本発明の表層流速推定方法の実行にあたっては、まず、制御部１１のデータ読込部１１ａは、処理対象領域のドップラースペクトルデータの読み込みを行う（Ｓ１）。

本発明では、海洋レーダ受信アンテナで受信した電波（即ちエコー）であって海洋波からの後方散乱波のドップラースペクトルのデータが用いられる。

本発明で用いるデータを収集するためのレーダは、距離方向・周波数方向で得られるドップラースペクトルを求めることができればどのような方式のものであっても良い。具体的には例えば、レーダの送受信方式が例えばＦＭＣＷ（Ｆrequency Ｍodulated Ｃontinuous Ｗave の略）方式やＦＭＩＣＷ（Ｆrequency Ｍodulated and Ｉnterrupted Ｃontinuous Ｗave の略）方式であるレーダ、空中線方式が例えば機械回転式やＤＢＦ（Ｄigital Ｂeam Ｆorming の略）方式やフェイズドアレイ方式であるレーダが用いられ得る。また、中心周波数が例えば２４．５１５ＭＨｚ若しくは４１．９ＭＨｚ、周波数掃引幅が例えば３００ｋＨｚ程度、送信出力が例えば５０Ｗ〜１００Ｗ程度、ビーム幅が例えば１３度〜２０度程度、距離分解能が例えば０．５ｋｍ程度のレーダが用いられ得る。なお、エコーからドップラースペクトルを求める方法自体は周知の技術であるのでここでは詳細については省略する（例えば、野崎憲朗・梅原俊彦：２４ＭＨｚ帯を使った短波海洋レーダによる海流、波浪の観察（前編），ＨＡＭ Ｊournal，６３，pp.77-85，１９８９年； 野崎憲朗・梅原俊彦：２４ＭＨｚ帯を使った短波海洋レーダによる海流、波浪の観察（後編），ＨＡＭ Ｊournal，６４，pp.103-108，１９８９年； 土木学会 海岸工学委員会編 研究現況レビュー小委員会：陸上設置型レーダによる沿岸海洋観測，２００１年； 藤井智史：海洋レーダの技術と歴史，沿岸海洋研究，４１，pp.97-108，２００４年）。

本実施形態では、同一の領域を異なる方向から観測することができるように観測対象領域の海洋沿岸に離間して設置された二つのレーダ局によって実際に受信されたエコーのドップラースペクトルを用いて説明する。本実施形態で用いるドップラースペクトルを図３に示す。なお、図３において、横軸は周波数（単位はＨｚ）、縦軸はビーム方向の距離（単位はｋｍ）、コンターの値はスペクトルの強度（単位はｄＢ×１００）をそれぞれ表す。なお、レーダ局は送信アンテナと受信アンテナと両アンテナを制御等する制御装置とを少なくとも有する。

本実施形態では、各レーダ局１回の観測で、８．２度間隔の８方位について、ビーム方向距離５００ｍ幅で３２ｋｍまで（即ち６４レンジ）のドップラースペクトルデータが収集される。そして、図３に示されるように、本実施形態で用いるドップラースペクトルは、ビーム方向距離８ｋｍ周辺にノイズがあると共に、ビーム方向距離１３ｋｍ周辺に見られるようにスペクトルが大きな位置においても複数のピークを有するものである。

本実施形態では、ドップラースペクトルデータは、スペクトルデータベース１８としてデータサーバ１６に予め保存される。なお、ドップラースペクトルデータは、観測領域別、観測日時別、レーダ別、ビーム方位別、ビーム方向距離別、周波数別にスペクトルデータベース１８に保存される。

ここで、各処理の説明において明示しないが、以降の各種の演算やメモリ１５への記憶やメモリ１５からの読み込みなどの全ての処理はレーダ毎・ビーム方位毎に行われる。

データ読込部１１ａは、スペクトルデータベース１８から処理対象領域の処理対象日時におけるビーム方向距離別・周波数別のドップラースペクトルデータを読み込み、メモリ１５に記憶させる。図３に示される本実施形態で用いるドップラースペクトルは、或る処理対象領域・日時における特定のレーダ局によって受信された特定のビーム方位のエコーのドップラースペクトルを周波数別・ビーム方向距離別に表示したものである。

次に、制御部１１のスペクトル分離部１１ｂは、Ｓ１の処理において読み込まれたドップラースペクトルデータを正の周波数領域のスペクトルデータと負の周波数領域のスペクトルデータとに分離する（Ｓ２）。

具体的には、スペクトル分離部１１ｂは、Ｓ１の処理においてメモリ１５に記憶された処理対象領域の処理対象日時におけるビーム方向距離別・周波数別のドップラースペクトルデータをメモリ１５から読み込み、周波数が正の領域のスペクトルＡi,jと負の領域のスペクトルＢi,jとに分離する。ここに、Ａ，Ｂはスペクトルの値（単位はｄＢ）、ｉはビーム方向距離の番号、ｊは周波数の番号をそれぞれ表す。

本実施形態では、ビーム方向距離の幅は５００ｍであるので、距離０．５ｋｍに番号１、距離１ｋｍに番号２、以降同様にして５００ｍ毎に順に連続番号を付与し、距離３２ｋｍに番号６４を付与する（例えば、坂井伸一 他：広域流動観測のための高性能沿岸海洋レーダの開発，研究報告Ｕ０２０５６，電力中央研究所，２００３年）。したがって、本実施形態では、ビーム方向距離番号ｉは１から６４までの自然数である。

また、本実施形態では、中心周波数が４１．９ＭＨｚであり、この中心周波数を挟んでプラスマイナス１．７５５３７５Ｈｚの周波数帯を受信している。ここで、本実施形態の説明では、断りがない限り中心周波数を０Ｈｚとして記述する。また、周波数のピッチは０．００５９６Ｈｚである。さらにまた、本実施形態では、−０．０９５３６Ｈｚ〜０．０９５３６Ｈｚの間を除いて処理している。そして、本実施形態では、正の領域に対しては周波数０．１０１３２Ｈｚに番号１、周波数０．１０７２８Ｈｚに番号２、以降同様にして０．００５９６Ｈｚ毎に順に連続番号を付与し、周波数１．２２７７６Ｈｚに番号１９１を付与する。また、負の領域に対しては周波数−１．２２７７６Ｈｚに番号１、周波数−１．２２１８０Ｈｚに番号２、以降同様にして０．００５９６Ｈｚ毎に順に連続番号を付与し、周波数−０．１０１３２Ｈｚに番号１９１を付与する。よって、本実施形態では、周波数番号ｊは１から１９１までの自然数である。

なお、ドップラースペクトルデータにおいて、周波数が正の領域のデータはレーダ局に近付いてくる海面波からの散乱応答のスペクトルデータであり、周波数が負の領域のデータはレーダ局から遠ざかる海面波からの散乱応答のスペクトルデータである。

そして、スペクトル分離部１１ｂは、正の領域スペクトルＡi,j及び負の領域スペクトルＢi,j（ただし、ｉ＝１，２，３，…，６４、並びに、ｊ＝１，２，３，…，１９１）をメモリ１５に記憶させる。

次に、制御部１１のスペクトル平均算出部１１ｃは、Ｓ２の処理において設定された正の領域スペクトル及び負の領域スペクトルについて周波数方向の平均値を算出する（Ｓ３）。

具体的には、スペクトル平均算出部１１ｃは、Ｓ２の処理においてメモリ１５に記憶された正の領域スペクトルＡi,j及び負の領域スペクトルＢi,jをメモリ１５から読み込み、正負の領域別に、ビーム方向距離番号ｉ毎のスペクトルの平均値（すなわち、周波数方向の平均値）である正の領域スペクトルＡi,jの平均値Ａ^-i,・及び負の領域スペクトルＢi,jの平均値Ｂ^-i,・を算出する。

そして、スペクトル平均算出部１１ｃは、正の領域スペクトル平均値Ａ^-i,・及び負の領域スペクトル平均値Ｂ^-i,・（ただし、ｉ＝１，２，３…，６４）をメモリ１５に記憶させる。

次に、制御部１１のスペクトル偏差算出部１１ｄは、Ｓ２の処理において設定された正の領域スペクトル及び負の領域スペクトルとＳ３の処理において算出された正の領域スペクトル平均値及び負の領域スペクトル平均値とを用いて平均値からの偏差を算出する（Ｓ４）。

具体的には、スペクトル偏差算出部１１ｄは、Ｓ２の処理においてメモリ１５に記憶された正の領域スペクトルＡi,j及び負の領域スペクトルＢi,jをメモリ１５から読み込むと共に、Ｓ３の処理においてメモリ１５に記憶された正の領域スペクトル平均値Ａ^-i,・及び負の領域スペクトル平均値Ｂ^-i,・をメモリ１５から読み込み、数式１ａによって正の領域スペクトル偏差Ａ¹i,jを算出すると共に数式１ｂによって負の領域スペクトル偏差Ｂ¹i,jを算出する。

（数１ａ） Ａ¹i,j＝Ａi,j−Ａ^-i,・
（数１ｂ） Ｂ¹i,j＝Ｂi,j−Ｂ^-i,・

このＳ４の処理は、ドップラースペクトルのビーム方向距離番号ｉ別・周波数番号ｊ別の値とビーム方向距離番号ｉ毎の平均値との差を算出するものであり、ドップラースペクトルの値は全体的に負であることが多いので、この後の処理の容易さを考慮してピーク周波数及びその周辺領域のスペクトルの値を正の値にすることを目的としている。

図４は、図３に示される本実施形態におけるドップラースペクトルデータに対してＳ２からＳ４までの処理を行って算出される正の領域スペクトル偏差Ａ¹i,j及び負の領域スペクトル偏差Ｂ¹i,jである。図４（Ａ）は周波数が負の領域のデータであり、レーダ局から遠ざかる海面波からの散乱応答のスペクトルデータである。また、図４（Ｂ）は周波数が正の領域のデータであり、レーダ局に近付いてくる海面波からの散乱応答のスペクトルデータである。

図４において、縦軸は図３と同じくビーム方向距離（単位はｋｍ）であり、横軸は周波数と変換係数とを掛け合わせて算出されるビーム方向流速（単位はｃｍ／ｓ）である。なお、周波数と変換係数とを掛け合わせてビーム方向流速を算出する方法自体は周知の技術であるのでここでは詳細については省略する（例えば、坂井伸一 他：広域流動観測のための高性能沿岸海洋レーダの開発，研究報告Ｕ０２０５６，電力中央研究所，２００３年）。

そして、スペクトル偏差算出部１１ｄは、正の領域スペクトル偏差Ａ¹i,j及び負の領域スペクトル偏差Ｂ¹i,j（ただし、ｉ＝１，２，３，…，６４、並びに、ｊ＝１，２，３，…，１９１）をメモリ１５に記憶させる。

次に、制御部１１の部分領域偏差平均値算出部１１ｅは、Ｓ４の処理において算出された正の領域スペクトル偏差及び負の領域スペクトル偏差を用いて正の領域スペクトル偏差の部分領域の平均値及び負の領域スペクトル偏差の部分領域の平均値を算出する（Ｓ５）。

具体的には、部分領域偏差平均値算出部１１ｅは、Ｓ４の処理においてメモリ１５に記憶された正の領域スペクトル偏差Ａ¹i,j及び負の領域スペクトル偏差Ｂ¹i,jをメモリ１５から読み込む。

そして、数式２及び数式３によって、元のドップラースペクトルデータのビーム方向距離を小さい方から順にいくつかずつ括って元の距離幅よりも広い距離幅であると共に元のドップラースペクトルデータの周波数を小さい方から順にいくつかずつ括って元の周波数ピッチよりも広い周波数幅である部分領域毎に、正の領域スペクトル偏差Ａ¹i,jの部分領域平均値Ａ²k,l及び負の領域スペクトル偏差Ｂ¹i,jの部分領域平均値Ｂ²k,lを算出する。この際、元のデータのビーム方向距離のいくつかを括って作られる部分領域の距離レンジに対して距離が小さい方から順に部分領域の距離レンジ番号ｋ（＝１，２，３，…）を付与すると共に、元のデータの周波数のいくつかを括って作られる部分領域の周波数レンジに対して元の周波数番号が小さい方から順に部分領域の周波数レンジ番号ｌ（＝１，２，３，…）を付与する。

数式２及び数式３において、ist＝im×(ｋ−１)＋１及びjst＝jm×(ｌ−１)＋１である。imはドップラースペクトルデータの部分領域を形成する際の元のドップラースペクトルデータのビーム方向距離の括りの数であり、jmは元のドップラースペクトルデータの周波数の括りの数である。

正の領域及び負の領域スペクトル偏差の部分領域平均値Ａ²k,l及びＢ²k,lの算出は、元のドップラースペクトルデータの解像度を下げてスペクトルの値の変化を滑らかにするために行う処理である。したがって、ビーム方向距離の括りの数im及び周波数の括りの数jmは原則として２以上の自然数として設定される。具体的には、括りの数im，jmは、これらによって決定される部分領域のビーム方向距離の幅及び周波数（即ち流速）の幅がピーク周波数を捕捉する際の解像度となるので、元のドップラースペクトルデータのビーム方向距離ピッチ及び周波数ピッチと共に流速推定の目的に基づいて必要とされる推定精度等を考慮して適宜設定される。例えば、部分領域のビーム方向距離の幅としては１〜４ｋｍ程度とすることが考えられ、周波数を変換した流速の幅としては１０〜３０ｃｍ／ｓ程度とすることが考えられる。ただし、部分領域のビーム方向距離の幅及び流速の幅は、これらに限定されるものではなく、どちらについてもこれより小さい値であっても良いしこれより大きい値であっても良い。なお、例えば高い推定精度が必要とされる場合や元のドップラースペクトルデータのビーム方向距離のピッチが大きい場合などにはビーム方向距離を括らないことにしてビーム方向距離の括りの数im＝１としても良い。

本実施形態では、元データのビーム方向距離ピッチが５００ｍであるところ距離の括りの数im＝４に設定され、また、元データの周波数ピッチが０．００５９６Ｈｚであるところ周波数の括りの数jm＝７に設定される。これにより、本実施形態では、部分領域のビーム方向距離の幅は２ｋｍとなり、周波数を変換した流速の幅は約１４ｃｍ／ｓとなる。

そして、本実施形態では、ビーム方向距離番号ｉ＝１，２，３，…，６４であるので距離の括りの数im＝４とすると部分領域の距離レンジ番号ｋは１から１６までとなる。また、本実施形態では、数式２及び数式３のjstはjst＝jm×(ｌ−１)＋２とし、周波数番号ｊ＝１，２，３，…，１９１であるので周波数の括りの数jm＝７とすると部分領域の周波数レンジ番号ｌは１から２７までとなる。

図５は、図４に示される正の領域スペクトル偏差Ａ¹i,j及び負の領域スペクトル偏差Ｂ¹i,jに対してＳ５の処理を行って算出される正の領域スペクトル偏差の部分領域平均値Ａ²k,l（図５（Ｂ））及び負の領域スペクトル偏差の部分領域平均値Ｂ²k,l（図５（Ａ））である。図５における横軸及び縦軸は図４と同様である。

図５に示すように、Ｓ５の処理を行うことによって元のスペクトルデータと比べて値の変化が滑らかなデータになり、スペクトルのノイズが平準化される。

ビーム方向距離の括りの数im及び周波数の括りの数jmの値は、表層流速推定プログラム１７上に予め規定しておくと共にＳ５の処理を実行する段階でこれらの値を部分領域偏差平均値算出部１１ｅがプログラムから読み込むようにしても良いし、ビーム方向距離の括りの数im及び周波数の括りの数jmの値を作業者が指定すると共に指定された値を部分領域偏差平均値算出部１１ｅが読み込むようにしても良い。

作業者がビーム方向距離の括りの数im及び周波数の括りの数jmの値を指定する場合には、部分領域偏差平均値算出部１１ｅが、Ｓ５の処理を実行する段階でビーム方向距離の括りの数im及び周波数の括りの数jmの値の指定を要求する内容のメッセージを表示部１４に表示し、入力部１３を介して入力された作業者の指定の値を読み込むようにする。

そして、部分領域偏差平均値算出部１１ｅは、正の領域スペクトル偏差の部分領域平均値Ａ²k,l及び負の領域スペクトル偏差の部分領域平均値Ｂ²k,l（ただし、ｋ＝１，２，３，…，１６、並びに、ｌ＝１，２，３，…，２７）をメモリ１５に記憶させる。

次に、第二段階として部分領域スペクトルデータを用いてフィルタ行列を作成する処理（Ｓ６〜Ｓ１１）について説明する。

周波数とビーム方向距離とについて展開されるドップラースペクトルデータからピーク周波数を検出するためのフィルタとして当該ドップラースペクトルデータにおけるピークを含むと共に当該ピークと連接し且つ所定の閾値よりも大きいスペクトルからなる領域（連接領域と呼ぶ）を用いることが本発明の特徴であり、第二段階の処理はこのフィルタを作成するための処理である。

本実施形態では、ピークを含み所定の閾値よりも大きいスペクトルからなる連接領域をドップラースペクトルデータから抽出する方法として隣接行列を用いるようにしている。

制御部１１の部分領域偏差平均値合算部１１ｆは、Ｓ５の処理において算出された正の領域スペクトル偏差の部分領域平均値と負の領域スペクトル偏差の部分領域平均値とを部分領域毎に合算する（Ｓ６）。

具体的には、部分領域偏差平均値合算部１１ｆは、Ｓ５の処理においてメモリ１５に記憶された正の領域スペクトル偏差の部分領域平均値Ａ²k,l及び負の領域スペクトル偏差の部分領域平均値Ｂ²k,lをメモリ１５から読み込み、数式４によって部分領域毎に正の領域スペクトル偏差の部分領域平均値Ａ²k,lと負の領域スペクトル偏差の部分領域平均値Ｂ²k,lとを足し合わせて正の領域と負の領域とのスペクトル偏差の部分領域平均値の合算値Ｃ²k,l（以下適宜、部分領域偏差平均合算値Ｃ²k,lと表記する）を算出する。

（数４）Ｃ²k,l＝Ａ²k,l＋Ｂ²k,l

図６は、図５に示される正の領域スペクトル偏差の部分領域平均値Ａ²k,lと負の領域スペクトル偏差の部分領域平均値Ｂ²k,lとを足し合わせた部分領域偏差平均合算値Ｃ²k,lである。図６における横軸及び縦軸は図４と同様である。

図６に示すように、Ｓ６の処理を行うことによって他の領域に対するピーク周波数周辺の突出の程度がより顕著になる。

そして、部分領域偏差平均値合算部１１ｆは、部分領域偏差平均合算値Ｃ²k,l（ただし、ｋ＝１，２，３，…，１６、並びに、ｌ＝１，２，３，…，２７）をメモリ１５に記憶させる。

次に、制御部１１の二値化行列作成部１１ｇは、Ｓ６の処理において算出された部分領域偏差平均合算値を用いて部分領域毎に当該部分領域に対応する成分の値を二値化して二値化行列を作成する（Ｓ７）。

具体的には、二値化行列作成部１１ｇは、Ｓ６の処理においてメモリ１５に記憶された部分領域偏差平均合算値Ｃ²k,lをメモリ１５から読み込み、部分領域毎に部分領域偏差平均合算値Ｃ²k,lが二値化閾値Ｔbよりも大きい場合には当該部分領域に対応する成分ｆk,lの値を１とすると共に二値化閾値Ｔb以下の場合には当該部分領域に対応する成分ｆk,lの値を０として二値化行列Ｆ（成分ｆk,l）を作成する。

部分領域毎の二値化は、レーダ局を介して得られるドップラースペクトルデータにおける流速情報を含むスペクトル情報を抜き出すために行う処理である。そのため、二値化閾値Ｔbは元のドップラースペクトルデータにおける流速情報を含まないノイズの影響部分を排除することができる値に設定される。具体的には、ビーム方向距離の限界近辺の距離帯では電磁波が減衰して流速情報よりもノイズの影響が卓越していると考えられるので、当該距離帯におけるスペクトルの最大値を二値化閾値Ｔbとすることが考えられる。

本実施形態では、距離方向にスペクトルが減衰していると考えられるビーム方向距離３０ｋｍの位置における全周波数方向での部分領域偏差平均合算値Ｃ²k,lの最大値を二値化閾値Ｔbとして設定する。

また、二値化行列作成部１１ｇは、さらに、部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域を判別し、当該部分領域の距離レンジ番号ｋを最大部分領域距離レンジ番号ｋ’とすると共に当該部分領域の周波数レンジ番号ｌを最大部分領域周波数レンジ番号ｌ’とする。

図７は、図６に示される部分領域偏差平均合算値Ｃ²k,lに対してＳ７の処理を行って作成される部分領域に対応する成分ｆk,lからなる二値化行列Ｆである。図７において、横軸は周波数レンジ番号ｌであり、縦軸は距離レンジ番号ｋである。また、図中の枡目は周波数レンジ番号ｌと距離レンジ番号ｋとの組み合わせに対応する部分領域を表す。

そして、灰色の枡目は対応する部分領域の部分領域偏差平均合算値Ｃ²k,lが二値化閾値Ｔb以下であって二値化行列Ｆの成分ｆk,l＝０であることを表し、白色の枡目は対応する部分領域の合算値Ｃ²k,lが二値化閾値Ｔbよりも大きく成分ｆk,l＝１であることを表す。さらに、黒色の枡目は合算値Ｃ²k,lが最大となる部分領域であることを表す。したがって、黒色の枡目に対応する部分領域の距離レンジ番号ｋが最大部分領域距離レンジ番号ｋ’になると共に周波数レンジ番号ｌが最大部分領域周波数レンジ番号ｌ’になる。

本実施形態では、図７に示すように、最大部分領域距離レンジ番号ｋ’＝３であると共に最大部分領域周波数レンジ番号ｌ’＝２２である。

そして、二値化行列作成部１１ｇは、二値化行列Ｆの成分ｆk,l（ただし、ｋ＝１，２，３，…，１６、並びに、ｌ＝１，２，３，…，２７）の値並びに最大部分領域距離レンジ番号ｋ’及び最大部分領域周波数レンジ番号ｌ’をメモリ１５に記憶させる。

次に、制御部１１の隣接行列作成部１１ｈは、Ｓ７の処理において作成された二値化行列を用いて隣接行列を作成する（Ｓ８）。

具体的には、隣接行列作成部１１ｈは、まず、距離レンジ番号ｋの順で周波数レンジ番号ｌの順に部分領域番号ｍを部分領域毎に付与する。例えば、距離レンジ番号ｋ＝１で周波数レンジ番号ｌ＝１である部分領域には番号ｍ＝１を付与し、距離レンジ番号ｋ＝１で周波数レンジ番号ｌ＝２である部分領域には番号ｍ＝２を付与する。

そして、本実施形態では、周波数レンジ番号ｌの最大値は２７であるので、距離レンジ番号ｋ＝２で周波数レンジ番号ｌ＝１である部分領域には番号ｍ＝２８が付与され、距離レンジ番号ｋ＝１６（本実施形態の距離レンジ番号ｋの最大値である）で周波数レンジ番号ｌ＝１である部分領域には番号ｍ＝４０６が付与される。

隣接行列作成部１１ｈは、また、Ｓ７の処理においてメモリ１５に記憶された最大部分領域距離レンジ番号ｋ’及び最大部分領域周波数レンジ番号ｌ’をメモリ１５から読み込み、距離レンジ番号ｋが最大部分領域距離レンジ番号ｋ’であると共に周波数レンジ番号ｌが最大部分領域周波数レンジ番号ｌ’である部分領域の番号ｍを最大部分領域番号ｍ’とする。

そして、隣接行列作成部１１ｈは、距離レンジ番号ｋと周波数レンジ番号ｌとの組み合わせに対応する部分領域番号ｍ（ただし、ｋ＝１，２，３，…，１６、ｌ＝１，２，３，…，２７、並びに、ｍ＝１，２，３，…，４３２）及び最大部分領域番号ｍ’をメモリ１５に記憶させる。

次に、隣接行列作成部１１ｈは、Ｓ７の処理においてメモリ１５に記憶された二値化行列Ｆ、並びに、前述の処理においてメモリ１５に記憶された距離レンジ番号ｋ及び周波数レンジ番号ｌの組み合わせと部分領域番号ｍとの間の対応関係をメモリ１５から読み込む。

隣接行列作成部１１ｈは、まず、隣接行列Ｇの成分ｇm,mを全てゼロにする。そして、隣接行列作成部１１ｈは、部分領域の二値化行列Ｆの成分ｆk,lの値が１の場合は当該部分領域の距離レンジ番号ｋと周波数レンジ番号ｌとの組み合わせに対応する部分領域番号ｍについての隣接行列Ｇの成分ｇm,m＝１とする。

また、隣接行列作成部１１ｈは、二値化行列Ｆの成分ｆk,lの値が１である部分領域に対して距離レンジ番号がｋ±１若しくは周波数レンジ番号がｌ±１である上下左右の位置の部分領域の二値化行列Ｆの成分ｆk±1,l及びｆk,l±1の値を確認する。

そして、隣接行列作成部１１ｈは、二値化行列Ｆの成分ｆk±1,l若しくはｆk,l±1の値が１の場合には、部分領域番号ｍ、並びに、距離レンジ番号ｋ±１と周波数レンジ番号ｍとの組み合わせ若しくは距離レンジ番号ｍと周波数レンジ番号ｌ±１との組み合わせに対応する番号ｎについての隣接行列Ｇの成分ｇm,n＝１とする。

以上により、隣接行列作成部１１ｈは、隣接行列Ｇの成分ｇm,nの値を設定する。なお、番号ｎの設定方法は部分領域番号ｍの設定方法と同じである。

図８は、図７に示される二値化行列Ｆに対してＳ８の処理を行って作成される隣接行列Ｇである。図７と同様に図８における横軸は周波数レンジ番号ｌであり、縦軸は距離レンジ番号ｋである。上下左右で部分領域が隣接している場合即ち隣接行列Ｇの成分ｇm,n＝１の場合に枡目内の番号の間に棒線を表示している。

そして、隣接行列作成部１１ｈは、隣接行列Ｇの成分ｇm,n（ただし、ｍ＝１，２，３，…，４３２、並びに、ｎ＝１，２，３，…，４３２）の値をメモリ１５に記憶させる。

次に、制御部１１の定常解行列作成部１１ｉは、Ｓ８の処理において作成された隣接行列を用いて定常解行列を作成する（Ｓ９）。

具体的には、定常解行列作成部１１ｉは、Ｓ８の処理においてメモリ１５に記憶された隣接行列Ｇをメモリ１５から読み込み、隣接行列Ｇをｐ回掛け合わせた行列Ｇ^pの各成分ｇ^pm,nと１とを比較し、小さい値を新たな成分ｇ^pm,nとして採用する（数式５）。

（数５）ｇ^pm,n＝ｍｉｎ(１，行列Ｇ^p＝Ｇ^p-1Ｇの成分ｇ^pm,n)

そして、定常解行列作成部１１ｉは、隣接行列Ｇをｐ乗した行列Ｇ^pが行列Ｇ^p-1と同等になるまで隣接行列Ｇを掛け合わせる処理を行い、Ｇ^p＝Ｇ^p-1である行列Ｇ^pを定常解行列Ｇsとする。

そして、定常解行列作成部１１ｉは、定常解行列Ｇsの成分ｇ^pm,n（ただし、ｍ＝１，２，３，…，４３２、並びに、ｎ＝１，２，３，…，４３２）の値をメモリ１５に記憶させる。

次に、制御部１１のフィルタ基礎行列作成部１１ｊは、Ｓ７の処理において作成された二値化行列とＳ９の処理において作成された定常解行列とを用いてフィルタ基礎行列を作成する（Ｓ１０）。

具体的には、フィルタ基礎行列作成部１１ｊは、Ｓ７の処理においてメモリ１５に記憶された二値化行列ＦとＳ９の処理においてメモリ１５に記憶された定常解行列Ｇsとをメモリ１５から読み込む。そして、フィルタ基礎行列作成部１１ｊは、定常解行列Ｇsにおいて最大部分領域番号ｍ’を行とする成分ｇ^pm',nが１である列番号ｎに対応する距離レンジ番号ｋ及び周波数レンジ番号ｌを除く二値化行列Ｆの成分ｆk,lを全てゼロにしてフィルタ基礎行列Ｈ（成分ｈk,l）を作成する。

図９は、図８に示される隣接行列Ｇに対してＳ９及びＳ１０の処理を行って作成されるフィルタ基礎行列Ｈである。図７と同様に図９において横軸は周波数レンジ番号ｌであり縦軸は距離レンジ番号ｋである、また、図中の枡目は部分領域を表す。そして、灰色の枡目は当該部分領域に対応する成分ｈk,lの値がゼロであることを表し、白色の枡目は当該部分領域に対応する成分ｈk,lの値が１である（以下、フィルタ透過成分と呼ぶ）ことを表す。さらに、黒色の枡目は部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域であることを表す。

図９に示すように、Ｓ９及びＳ１０の処理を行うことによって、図７に示される二値化行列Ｆに対し、部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域（図では黒色の枡目）と連接する部分領域（即ち連接領域。図では白色の枡目）のみが成分ｈk,lの値が１として（即ちフィルタ透過成分として）残る一方で、合算値Ｃ²k,lが最大となる部分領域と連接していない部分領域の成分ｈk,lの値がゼロになり、スペクトルのノイズが除去される。

そして、フィルタ基礎行列作成部１１ｊは、フィルタ基礎行列Ｈの成分ｈk,l（ただし、ｋ＝１，２，３，…，１６、並びに、ｌ＝１，２，３，…，２７）の値をメモリ１５に記憶させる。

次に、制御部１１のフィルタ行列作成部１１ｋは、Ｓ１０の処理において作成されたフィルタ基礎行列を用いてフィルタ行列を作成する（Ｓ１１）。

具体的には、フィルタ行列作成部１１ｋは、Ｓ１０の処理においてメモリ１５に記憶されたフィルタ基礎行列Ｈをメモリ１５から読み込む。そして、フィルタ行列作成部１１ｋは、フィルタ基礎行列Ｈにおいて値が１である成分（即ちフィルタ透過成分）ｈk,lが透過成分最大周波数幅ｌmax個よりも多く連接している行ｋについては部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域の位置から左右(ｌmax−１)／２個の成分を除く成分の値をゼロとしてフィルタ行列Ｘ（成分ｘk,l）を作成する。

透過成分最大周波数幅ｌmaxは、元のドップラースペクトルデータからピーク周波数を含む領域を抽出する際の周波数方向の幅即ちフィルタ透過成分に対応する領域の周波数方向の拡がりの程度を定義づけるものであり、部分領域の個数として定義される。透過成分最大周波数幅ｌmaxの値は、特定の値に限定されるものではなく、例えば周波数を流速に換算した場合の流速の幅を考慮して適宜設定される。

本実施形態では、流速に換算した場合にピーク周波数の両側で約５０ｃｍ／ｓ程度の拡がりを有するフィルタにすることを目的とし、流速に換算した場合に約１００ｃｍ／ｓとなるように透過成分最大周波数幅ｌmax＝７に設定される。なお、本実施形態のように、ピークを中心とする全体の幅として透過成分最大周波数幅ｌmaxを定義する場合にはｌmaxの値は奇数に設定される。

図１０は、図９に示されるフィルタ基礎行列Ｈに対してＳ１１の処理を行って作成されるフィルタ行列Ｘである。図７と同様に図１０において横軸は周波数レンジ番号ｌであり縦軸は距離レンジ番号ｋであると共に図中の枡目は部分領域を表している。そして、灰色の枡目は当該部分領域に対応する成分ｘk,lの値がゼロであることを表し、白色の枡目は当該部分領域に対応する成分ｘk,lの値が１であることを表す。さらに、黒色の枡目は部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域であることを表す。

図１０に示すように、Ｓ１１の処理を行うことによって、図９に示されるフィルタ基礎行列Ｈと比べて成分ｘk,lの値が１である部分領域（図では白色の枡目）が部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域（図では黒色の枡目）から一定の範囲に限定され、スペクトルのノイズが更に除去される。

そして、フィルタ行列作成部１１ｋは、フィルタ行列Ｘの成分ｘk,l（ただし、ｋ＝１，２，３，…，１６、並びに、ｌ＝１，２，３，…，２７）の値をメモリ１５に記憶させる。

続いて、第三段階としてフィルタ行列を用いてピーク周波数を検出する処理（Ｓ１２〜Ｓ１４）及び流速を推定する処理（Ｓ１５）について説明する。

制御部１１のスペクトル偏差合算部１１ｌは、Ｓ４の処理において算出された正の領域スペクトル偏差と負の領域スペクトル偏差とをビーム方向距離毎・周波数毎に合算する（Ｓ１２）。

具体的には、スペクトル偏差合算部１１ｌは、Ｓ４の処理においてメモリ１５に記憶された正の領域スペクトル偏差Ａ¹i,j及び負の領域スペクトル偏差Ｂ¹i,jをメモリ１５から読み込み、これら正の領域スペクトル偏差Ａ¹i,jと負の領域スペクトル偏差Ｂ¹i,jとを足し合わせてスペクトル偏差合算値Ｃ¹i,jを算出する。

図１１は、図４に示される正の領域スペクトル偏差Ａ¹i,jと負の領域スペクトル偏差Ｂ¹i,jとを足し合わせたスペクトル偏差合算値Ｃ¹i,jである。図１１の横軸及び縦軸は図４と同様である。

そして、スペクトル偏差合算部１１ｌは、スペクトル偏差合算値Ｃ¹i,jをメモリ１５に記憶させる。

次に、制御部１１のフィルタ適用部１１ｍは、Ｓ１１の処理において作成されたフィルタ行列をＳ１２の処理において算出されたスペクトル偏差合算値に適用する（Ｓ１３）。

具体的には、フィルタ適用部１１ｍは、Ｓ１１の処理においてメモリ１５に記憶されたフィルタ行列Ｘをメモリ１５から読み込むと共にＳ１２の処理においてメモリ１５に記憶されたスペクトル偏差合算値Ｃ¹i,jをメモリ１５から読み込み、スペクトル偏差合算値Ｃ¹i,jにフィルタ行列Ｘを適用する。

すなわち、フィルタ適用部１１ｍは、フィルタ行列Ｘにおいて値が１になっている成分ｘk,lの距離レンジ番号ｋ及び周波数レンジ番号ｌに対応するビーム方向距離番号ｉ及び周波数番号ｊのスペクトルのみスペクトル偏差合算値Ｃ¹i,jの値とするフィルタ適用スペクトルＤi,jを作成する。なお、本実施形態では、フィルタ適用スペクトルＤi,jのｉ＝１及び１９１については、フィルタ行列Ｘのｋ＝１及び２７の成分ｘk,lの値を用いる。

図１２は、図１０に示されるフィルタ行列Ｘを図１１に示されるスペクトル偏差合算値Ｃ¹i,jに適用した結果得られるフィルタ適用スペクトルＤi,jである。図１２の横軸及び縦軸は図４と同様である。

図１２に示すように、Ｓ１３の処理を行うことによって局所的にスペクトルの強度が非常に強くなったり強度がばらついたりという外因性ノイズが除去され、元のスペクトルデータから正当なピークを含み且つそれと連接するスペクトルを有する領域のみが残る。

そして、フィルタ適用部１１ｍは、フィルタ適用スペクトルＤi,jをメモリ１５に記憶させる。

次に、制御部１１のピーク周波数検出部１１ｎは、Ｓ１３の処理において作成されたフィルタ適用スペクトルを用いてピーク周波数を検出する（Ｓ１４）。

具体的には、ピーク周波数検出部１１ｎは、Ｓ１３の処理においてメモリ１５に記憶されたフィルタ適用スペクトルＤi,jを読み込み、フィルタ適用スペクトルＤi,jの周波数方向のモーメントＤ^mi,・及び平均値Ｄ^-i,・をビーム方向距離番号ｉ毎に算出し、これらを用いて周波数方向の重心を算定する。なお、この重心の算定自体はモーメント法として周知の技術であるのでここでは詳細については省略する（例えば、亀田弘行ら：確率・統計解析，技報堂出版，１９８１年）。

そして、ピーク周波数検出部１１ｎは、ビーム方向距離番号ｉ毎の重心をピーク周波数としてメモリ１５に記憶させる。

次に、制御部１１の流速算出部１１ｏは、Ｓ１４の処理において検出されたピーク周波数を用いてビーム方向流速を算出する（Ｓ１５）。

具体的には、流速算出部１１ｏは、Ｓ１４の処理においてメモリ１５に記憶されたビーム方向距離番号ｉ毎のピーク周波数を読み込み、このピーク周波数と変換係数とを掛け合わせてビーム方向距離番号ｉ毎のビーム方向流速を算出する。周波数と変換係数とを掛け合わせてビーム方向流速を算出する方法自体は周知の技術であるのでここでは詳細については省略する（例えば、坂井伸一 他：広域流動観測のための高性能沿岸海洋レーダの開発，研究報告Ｕ０２０５６，電力中央研究所，２００３年）。なお、変換係数は表層流速推定プログラム１７上に予め規定しておく。

そして、流速算出部１１ｏは、ビーム方向距離番号ｉ毎のビーム方向流速をメモリ１５に記憶させる。

以上の処理によってレーダ別・ビーム方位別のビーム方向距離番号毎のビーム方向流速即ち視線流速が推定される。

表層流速推定装置１０は、流速推定の目的や推定結果の活用目的に合わせて適宜、推定された流速を推定結果データとしてデータサーバ１６に蓄積したり処理対象領域の地形図と重ね合わせて表示部１４に表示したりする。この際、Ｓ１５までの処理によって推定されるのはレーダ別のビーム方向流速即ち視線流速であるので、表層流速推定装置１０はこのレーダ別の視線流速をビーム方位毎・ビーム方向距離毎にベクトル合成して地点別の水平流速を算出するようにしても良い。

以上の構成を有する本発明の表層流速推定方法、装置並びにプログラムによれば、従来の技術と比較して誤差（ノイズ）情報の検定処理が優れており、その効果により流動情報を適確に検出することができる。そのため、本発明は、自然の雑音や人口雑音等の外因性ノイズが含まれる実際の観測データに対して推定精度の低下や後処理の煩雑さを避けることが可能になる。その結果、本手法により、リアルタイム海域流動モニタリングに資する流動の情報を短時間で提供することが可能になる

なお、上述の形態は本発明の好適な形態の一例ではあるがこれに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々変形実施可能である。例えば、本実施形態では、観測対象領域の海洋沿岸に設置された二つのレーダ局によって受信されたエコーのドップラースペクトルデータを表層流速推定装置１０と通信回線等によって接続されたデータサーバ１６内にスペクトルデータベース１８として一旦蓄積させるようにしているが、これに限られず、レーダ局からのデータを例えばＰＨＳを利用した無線通信回線を介して表層流速推定装置１０に直接入力するようにしても良い。この場合にはほぼリアルタイムでの推定及び情報提供が可能になる。

また、本実施形態では、二つのレーダ局によって受信されたエコーのドップラースペクトルデータをスペクトルデータベース１８として一つに集約し当該データを用いて表層流速推定装置１０が流速の推定を行うと共に必要に応じて当該推定結果を用いて情報提供を行うようにしているが、これに限られず、表層流速推定装置１０を各レーダ局に設置して各レーダ局で流速の推定を行うと共に、流速データを用いて情報提供を行う装置を更に設けて当該装置に推定された流速データを有線・無線の通信回線を介して転送し集約して当該装置が情報提供を行うようにしても良い。

また、本実施形態では、ピーク周波数を検出するためのフィルタ行列を作成する際に隣接行列を用いるようにしている（Ｓ８〜Ｓ１１）が、周波数とビーム方向距離とについて展開されるドップラースペクトルデータにおけるピークを含むと共に当該ピークと連接し且つ所定の閾値よりも大きいスペクトルからなる領域を抽出することができる方法であれば、フィルタ行列を作成する方法は本実施形態の方法に限られるものではない。具体的には、本実施形態のＳ８の処理において、図７に示される二値化行列Ｆに対して隣接行列Ｇを作成するようにしているが、これに限られず、例えば画像処理のラベリング処理を用いるようにしても良い。

この場合には、まず、図７に示される二値化行列Ｆの白い枡目即ち部分領域偏差平均合算値Ｃ²k,lが二値化閾値Ｔbよりも大きい部分領域が連接する連接領域を画像処理によって抽出する。続いて、抽出された連接領域毎に異なるラベルを当該連接領域に含まれる部分領域に付与する。そして、部分領域偏差平均合算値Ｃ²k,lが最大となる部分領域（即ち本実施形態では最大部分領域距離レンジ番号ｋ’且つ最大部分領域周波数レンジ番号ｌ’の部分領域）に付与されたラベルと同じラベルの部分領域に対応する成分を１にすると共に異なるラベルの部分領域に対応する成分を０にしてフィルタ行列Ｈを作成する。

本発明の表層流速推定方法、装置並びにプログラムを大阪湾明石海峡付近での実際の海域観測の結果のデータを用いた流速の推定に適用した実施例を図１３から図１６を用いて説明する。

本実施例では、大阪湾明石海峡領域において二台のデジタル・ビーム・フォーミング方式レーダを用いて得られた海域観測結果を用いた。二台のレーダの設置点及びビーム方位は図１３に示す通りであり、図中のＳt.Ａ及びＳt.Ｂが付された記号●がレーダの設置点、当該記号●のそれぞれから放射状に出ている直線がビーム方位である。海域観測は２００３年１月１８日から同年２月２０日までの約一ヶ月に亘って行われた。

さらに、本発明の表層流速推定方法、装置並びにプログラムによる流速の推定結果と比較するため、上記観測期間中の１月３１日及び２月１７日に、船舶に取り付けたＡＤＣＰ（Ａcoustic Ｄoppler Ｃurrent Ｐrofiler：超音波ドップラー流速計）を用いて海域の流動を直接観測した。観測に使用したＡＤＣＰは６００ｋＨｚのワークホース及びブロードバンドタイプであり、測定層厚５０ｃｍ間隔で１０秒間隔１５ｐｉｎｇ（１５回分）の平均流速を収録した。なお、仕様上の標準偏差は２．３ｃｍであった。

ＡＤＣＰを用いた流速観測は三隻の船舶の側面にＡＤＣＰを設置して行った。ＡＤＣＰを用いた流速の測定点は図１４に示す通りである。なお、ＡＤＣＰを用いた流速の測定時刻は、１月３１日は７時から１３時にかけて１２測点３潮時であり、２月１７日は１３時３０分から１６時３０分にかけて１５測点２潮時であった。

ＡＤＣＰを用いた流速観測は、１月３１日は図１４（Ａ）におけるＡ−Ａ’，Ｂ−Ｂ’，Ｃ−Ｃ’ラインを、２月１７日は図１４（Ｂ）におけるＤ−Ｄ’，Ｅ−Ｅ’，Ｆ−Ｆ’ラインを異なる船舶で移動し、Ｓt.Ａ及びＳt.Ｂの二局のレーダのビーム方位が重なる地点（図中の記号○）で行った。そして、一地点の計測では、数１０分間連続した流速結果を収録した。なお、本実施例における本発明による推定結果とＡＤＣＰを用いた測定結果との比較では、ＡＤＣＰを用いた観測によって得られた流速については機器の仕様上最上層である水深約２ｍにおける結果を用いた。

レーダによって得られた実際のスペクトルデータを用いて本発明によって流速の推定を行い、図１５に示す結果が得られた。なお、図１５は１月１８日から３１日までの間の平均流速である。図１５において、図中の矢印の向きは当該矢印の位置における流れの向きを表し、矢印の長さは当該矢印の位置における流れの速さを表す。

そして、本発明によって推定されたビーム方向流速^ＲＶとＡＤＣＰを用いて測定された流速^ＡＶとの一致具体を比較するために両者の散布図を作成して図１６に示す結果が得られた。図１６において横軸は本発明による推定結果の値であって縦軸はＡＤＣＰを用いた測定結果の値であり、両者が等しくなっている場合の目安のために図１６中には縦軸の値＝横軸の値となる直線（ｙ＝ｘとなる破線）が示されている。なお、ＡＤＣＰを用いて測定された流速はレーダの視線方向流速に変換した結果を用いた。また、本発明による推定結果もＡＤＣＰを用いた測定結果も、１５分間の平均流速を用い、１月３１日と２月１７日とのデータを合わせたものである。

図１６から、本発明による推定結果とＡＤＣＰを用いた測定結果とはｙ＝ｘ近傍に分布しており、本発明による推定結果はＡＤＣＰを用いた測定結果と同様の結果を示すことが確認された。また、本発明による推定結果とＡＤＣＰを用いた測定結果との相関係数は０．９８となり、非常に高い相関となった。

本実施例の結果から、本発明による推定結果はＡＤＣＰを用いた実際の測定の結果と非常に近い結果となり、本発明による推定は信頼性が非常に高いことが確認された。したがって、本発明の方法によれば適確な流動情報を検出することができるので、例えばピーク周波数の検出結果のチェック等の煩雑な後処理が必要とされず推定に係る手間を省くことが可能になると共に推定を短時間で行うことが可能になることが確認された。

また、本実施例における本発明による推定について、一時点の推定にかかる処理時間は３０秒程度であった（処理に用いたＰＣは１ＣＰＵ：Ｐｅｎｔｉｕｍ４ ２．４ＧＨｚ；Ｐｅｎｔｉｕｍは登録商標）。この結果から、本発明は、流速の推定を非常に迅速に行うことが可能であり、例えば遠隔流況モニタリングに適していること、且つ、ほとんどリアルタイムでの流況に関する情報提供を行うことが可能であることが確認された。

本発明の表層流速推定方法の実施形態の一例を説明するフローチャートである。 実施形態の表層流速推定方法をプログラムを用いて実施する場合の表層流速推定装置の機能ブロック図である。 実施形態のドップラースペクトルを示す図である。 実施形態のスペクトル偏差を示す図である。（Ａ）は負の領域スペクトル偏差を示す図である。（Ｂ）は正の領域スペクトル偏差を示す図である。 実施形態のスペクトル偏差の部分領域平均値を示す図である。（Ａ）は負の領域スペクトル偏差の部分領域平均値を示す図である。（Ｂ）は正の領域スペクトル偏差の部分領域平均値を示す図である。 実施形態の正の領域スペクトル偏差の部分領域平均値と負の領域スペクトル偏差の部分領域平均値とを足し合わせた合算値を示す図である。 実施形態の部分領域に対応する二値化行列を示す図である。 実施形態の隣接行列を示す図である。 実施形態のフィルタ基礎行列を示す図である。 実施形態のフィルタ行列を示す図である。 実施形態のスペクトル偏差合算値を示す図である。 実施形態のフィルタ適用スペクトルを示す図である。 実施例のレーダの設置点及びビーム方位を示す図である。 実施例のＡＤＣＰを用いた流速の測定点を示す図である。（Ａ）は１月３１日の流速観測における測定点を示す図である。（Ｂ）は２月１７日の流速観測における測定点を示す図である。 実施例の本発明による流速の推定結果を示す図である。 実施例の本発明によって推定された流速とＡＤＣＰを用いて測定された流速との比較のための散布図である。

符号の説明

１０ 表層流速推定装置
１１ 制御部
１２ 記憶部
１３ 入力部
１４ 表示部
１５ メモリ
１６ データサーバ
１７ 表層流速推定プログラム
１８ スペクトルデータベース

Claims (4)

1. ビーム方向距離別・周波数別ドップラースペクトルのピーク周波数を検出して表層流速を推定する方法において、前記ドップラースペクトルの正の領域のスペクトルから該正の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出すると共に前記ドップラースペクトルの負の領域のスペクトルから該負の領域のスペクトルの周波数方向の平均値を引いて算出される偏差について部分領域の平均値を算出し、前記部分領域毎の前記正の領域のスペクトルの偏差の部分領域の平均値と前記負の領域のスペクトルの偏差の部分領域の平均値との合算値が最大である部分領域と連接し且つ前記部分領域毎の合算値が所定の閾値よりも大きい部分領域からなる連接領域の中から前記ピーク周波数を検出することを特徴とする表層流速推定方法。
2. 処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行うステップと、前記ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離するステップと、前記正の領域のスペクトル及び前記負の領域のスペクトルについて周波数方向の平均値を算出するステップと、前記正の領域のスペクトルから前記正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に前記負の領域のスペクトルから前記負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出するステップと、前記正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に前記負の領域のスペクトルの偏差について部分領域の平均値を算出するステップと、前記正の領域のスペクトルの偏差の部分領域の平均値と前記負の領域のスペクトルの偏差の部分領域の平均値とを前記部分領域毎に合算して部分領域毎の合算値を算出するステップと、前記部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成するステップと、前記二値化行列について隣接行列を作成するステップと、前記隣接行列を定常状態になるまで掛け合わせて定常解行列を作成するステップと、前記定常解行列から前記部分領域毎の合算値が最大である部分領域と連接し且つ前記部分領域毎の合算値が前記二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して前記二値化行列からフィルタ基礎行列を作成するステップと、前記フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成するステップと、前記正の領域のスペクトルの偏差と前記負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出するステップと、前記フィルタ行列を前記スペクトル偏差合算値に適用してフィルタ適用スペクトルを作成するステップと、前記フィルタ適用スペクトルからピーク周波数を検出するステップと、前記ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出するステップとを有することを特徴とする表層流速推定方法。
3. 処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行う手段と、前記ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離する手段と、前記正の領域のスペクトル及び前記負の領域のスペクトルについて周波数方向の平均値を算出する手段と、前記正の領域のスペクトルから前記正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に前記負の領域のスペクトルから前記負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出する手段と、前記正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に前記負の領域のスペクトルの偏差について部分領域の平均値を算出する手段と、前記正の領域のスペクトルの偏差の部分領域の平均値と前記負の領域のスペクトルの偏差の部分領域の平均値とを前記部分領域毎に合算して部分領域毎の合算値を算出する手段と、前記部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成する手段と、前記二値化行列について隣接行列を作成する手段と、前記隣接行列を定常状態になるまで掛け合わせて定常解行列を作成する手段と、前記定常解行列から前記部分領域毎の合算値が最大である部分領域と連接し且つ前記部分領域毎の合算値が前記二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して前記二値化行列からフィルタ基礎行列を作成する手段と、前記フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成する手段と、前記正の領域のスペクトルの偏差と前記負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出する手段と、前記フィルタ行列を前記スペクトル偏差合算値に適用してフィルタ適用スペクトルを作成する手段と、前記フィルタ適用スペクトルからピーク周波数を検出する手段と、前記ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出する手段とを有することを特徴とする表層流速推定装置。
4. ビーム方向距離別・周波数別ドップラースペクトルのピーク周波数を検出して表層流速の推定を行う際に、少なくとも、処理対象領域のビーム方向距離別・周波数別ドップラースペクトルデータの読み込みを行う手段、前記ドップラースペクトルデータを周波数が正の領域のスペクトルと負の領域のスペクトルとに分離する手段、前記正の領域のスペクトル及び前記負の領域のスペクトルについて周波数方向の平均値を算出する手段、前記正の領域のスペクトルから前記正の領域のスペクトルの周波数方向の平均値を引いて正の領域のスペクトルの偏差を算出すると共に前記負の領域のスペクトルから前記負の領域のスペクトルの周波数方向の平均値を引いて負の領域のスペクトルの偏差を算出する手段、前記正の領域のスペクトルの偏差について部分領域の平均値を算出すると共に前記負の領域のスペクトルの偏差について部分領域の平均値を算出する手段、前記正の領域のスペクトルの偏差の部分領域の平均値と前記負の領域のスペクトルの偏差の部分領域の平均値とを前記部分領域毎に合算して部分領域毎の合算値を算出する手段、前記部分領域毎の合算値が二値化閾値よりも大きいか否かによって当該部分領域に対応する成分の値を二値化して二値化行列を作成する手段、前記二値化行列について隣接行列を作成する手段、前記隣接行列を定常状態になるまで掛け合わせて定常解行列を作成する手段、前記定常解行列から前記部分領域毎の合算値が最大である部分領域と連接し且つ前記部分領域毎の合算値が前記二値化閾値よりも大きい部分領域からなる連接領域に対応する成分を抽出して前記二値化行列からフィルタ基礎行列を作成する手段、前記フィルタ基礎行列についてフィルタ透過成分の行毎の周波数方向の広がりを制限してフィルタ行列を作成する手段、前記正の領域のスペクトルの偏差と前記負の領域のスペクトルの偏差とをビーム方向距離毎・周波数毎に合算してスペクトル偏差合算値を算出する手段、前記フィルタ行列を前記スペクトル偏差合算値に適用してフィルタ適用スペクトルを作成する手段、前記フィルタ適用スペクトルからピーク周波数を検出する手段、前記ピーク周波数と変換係数とを掛け合わせてビーム方向流速を算出する手段としてコンピュータを機能させるための表層流速推定プログラム。