JP2008269097A - 水面波挙動検出装置、水面波挙動検出方法、及び、水面波挙動検出プログラム - Google Patents

Abstract

【課題】水面波を捉えた画像から水面波の挙動を検出することを課題とする。
【解決手段】水面波を含む画像を入力する画像入力部１００と、入力された画像を時系列に蓄積する画像蓄積部１１０と、波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、当該輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出する目的関数導出部１２０と、連続する２枚の画像を画像蓄積部１１０から読み出して、この２枚の画像の強度値を目的関数に代入し、速度ベクトルを推定するオプティカルフロー推定部１３０と、推定された速度ベクトルを表示する表示部１４０とを備える。
【選択図】図１

Description

本発明は、水面波の挙動を検出する技術に関する。

川や海などでは、水害事故を防止する為、水位センサや水量センサなどを用いて、特異な現象を事前に察知する環境センシングが行われている。特に、海では、ミリ波レーダを用いて、海面に浮かばせたブイからの反射波に基づき、波の状態を把握する観測方法も採用されている。

最近では、監視カメラの導入が急増しており、監視員や管理人が監視センターで川や海の様子を目視する方法も利用されている。監視カメラの設置台数は増加傾向にあるものの、人件費に係る費用を抑制するため、必ずしも監視員等の増員には至っていない。そのため、川や海の水面の状態を映像や画像から検出し、自動で認識可能な観測方法が切望されている。

その一つが画像処理を用いた観測方法である。例えば、ＭＰＥＧ（Moving Picture Experts Group）に代表される符号化の分野で用いられている相互相関法（「ＣＣ法」、「パターンマッチング法」とも称される）が知られている。ＣＣ法は、連続する２枚の画像間の類似性を相互相関関数等により評価して、類似する領域を画像間の移動距離（ベクトル）とみなす方法である。しかしながら、ＣＣ法は、画像間で輝度変動がないことを前提条件とする測定方法であるため、輝度変動が激しい水面観測には適用性に乏しい。そこで、以下のような様々な観測方法や推定方法が利用されている。

非特許文献１では、カメラで撮影した映像や画像の輝度変動を利用して、オプティカルフロー法に基づく動き推定方法が試みられている。

また、非特許文献２では、ビデオ映像から動きに関するいくつかの物理パラメータを推定するための枠組みが提案されている。具体的には、例えば、映像中の揺れる木や浮かんでいるヨットの動きを解析することによって、風力による外力から調和振動方程式を構成する各パラメータを求めている。

更に、非特許文献３では、河川工学や海洋工学の分野において、川や海に関する様々な波の動きモデルが提案されている。
J.L. Barron、外２名、「Performance of Optical Flow Techniques」、IJCV、1994、vol.12、no.1、p.43-77 M. Sun、外２名、「Video input driven animation (VIDA)」、IEEE International Conference on Compute Vision (ICCV’03)、2003、vol.1、p.96-103 堀川 清司、「湾岸工学」、東京大学出版会、1991、p.42-81

しかしながら、非特許文献１に記載されたオプティカルフロー法は、映像や画像の輝度変動のみ用いて動きを検出する方法であり、波動現象に関する動きモデルが考慮されていないため、水面波の挙動を検出できないという問題があった。

また、非特許文献２に記載された推定方法によれば、推定結果から風と波数を間接的に推定できるものの、規則的な波動を前提とする波モデルであるため、振幅や位相などの波に関する固有の動き物理パラメータ、即ち、波の不規則性が検出されないという問題があった。

更に、非特許文献３で提案された波生成方程式は、不規則性を表現する為に多くの線形和で表現されているため、各パラメータを演算する時間が膨大になるという問題があった。

本発明は、上記を鑑みてなされたものであり、水面波を捉えた画像から水面波の挙動を検出することを課題とする。

第１の本発明に係る水面波挙動検出装置は、水面波を含む画像を入力する画像入力手段と、入力された前記画像を時系列に蓄積する画像蓄積手段と、波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、当該輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出する目的関数導出手段と、連続する２枚の画像を前記画像蓄積手段から読み出して、当該２枚の画像の強度値を前記目的関数に代入し、前記速度ベクトルを推定するオプティカルフロー推定手段と、推定された当該速度ベクトルを表示する表示手段と、を有することを要旨とする。

本発明にあっては、波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、この輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出し、連続する２枚の画像の強度値をこの目的関数に代入して速度ベクトルを推定するため、輝度変動に基づくオプティカルフロー法を用いて、水面波を捉えた画像から不規則な水面波の挙動を検出することができる。

第２の本発明に係る水面波挙動検出装置は、前記波生成式が、波の振幅，波数，周波数，及び，位相を変数とする波生成式の数を項数とした場合に、当該項数の総和で示された波生成式であって、前記オプティカルフロー推定手段は、前記２枚の画像を複数のサブブロックにそれぞれ分割する分割手段と、当該２枚の画像において同じ場所に位置するサブブロックの強度値を前記目的関数に代入して、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する第１の推定手段と、所定の前記項数の範囲で前記第１の推定手段の推定を繰り返し、推定された前記周波数に基づいて波のスペクトル分布を求め、当該波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する第２の推定手段と、を有することを要旨とする。

本発明にあっては、２枚の画像を複数のサブブロックにそれぞれ分割し、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する第１の推定手段と、所定の項数の範囲で第１の推定手段の推定を繰り返し、推定された周波数に基づいて波のスペクトル分布を求め、この波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する第２の推定手段とを有するため、水面波の挙動をより正確に検出することができる。即ち、第１の推定手段で推定された局所的な速度ベクトルは、隣接するサブブロック間の連続性を考慮していないが、第２の推定手段において、推定された周波数から波のスペクトル分布を求め、この波のスペクトル分布と所定のスペクトル分布との類似度を判定し、所定の項数の範囲で最も類似度の高いスペクトル分布の速度ベクトル推定するため、画像全体の連続性を考慮した水面波の挙動を検出することが可能となる。

第３の本発明に係る水面波挙動検出装置は、前記第２の推定手段が、前記波のスペクトル分布の面積を１になるように正規化し、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布の速度ベクトルを推定することを要旨とする。

本発明にあっては、波のスペクトル分布の面積を１になるように正規化し、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布の速度ベクトルを推定するため、より正確に水面波の挙動を検出することができる。

第４の本発明に係る水面波挙動検出装置は、前記第２の推定手段が、前記速度拘束値及び／又は前記サブブロックの大きさを更に変化させて前記第１の推定手段の推定を繰り返し、最も類似度の高いスペクトル分布の速度ベクトルを推定することを要旨とする。

本発明にあっては、項数に加えて、速度拘束値及び／又はサブブロックの大きさを更に変化させて速度ベクトルなどを推定し、類似度の高いスペクトル分布の速度ベクトルを推定するため、より正確に水面波の挙動を検出することができる。

第５の本発明に係る水面波挙動検出装置は、前記所定のスペクトル分布が、ブレットシュナイダーのスペクトル分布であることを要旨とする。

本発明にあっては、波の統計分布として信頼性の高いブレットシュナイダーのスペクトル分布を用いて類似度を判定するため、信頼性の高い水面波の挙動を検出することができる。

第６の本発明に係る水面波挙動検出装置は、前記波生成式が、振幅を前記項数で除算し、波数，周波数，及び，位相を前記項数で乗算した波生成式であることを要旨とする。

本発明にあっては、振幅を項数で除算し、波数，周波数，及び，位相を項数で乗算した波生成式であるため、高い周波数の波を少数項でも確実に得ることができる。

第７の本発明に係る水面波挙動検出装置は、前記オプティカルフロー推定手段が、共役勾配法を用いて前記速度ベクトルを推定すること要旨とする。

本発明にあっては、共役勾配法を用いて速度ベクトルを推定するため、目的関数の最小化をより効率的に計算することができる。

第８の本発明に係る水面波挙動検出装置は、前記目的関数導出手段が、前記ロバスト関数に対してローレンツ関数を用いることを要旨とする。

本発明にあっては、ロバスト関数に対してローレンツ関数を用いるため、目的関数の計算速度を向上することができる。

第９の本発明に係る水面波挙動検出方法は、水面波を含む画像を入力する第１のステップと、入力された前記画像を時系列に画像蓄積手段に蓄積する第２のステップと、波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、当該輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出する第３のステップと、連続する２枚の画像を前記画像蓄積手段から読み出して、当該２枚の画像の強度値を前記目的関数に代入し、前記速度ベクトルを推定する第４のステップと、推定された当該速度ベクトルを表示する第５のステップと、を有することを要旨とする。

第１０の本発明に係る水面波挙動検出方法は、前記波生成式が、波の振幅，波数，周波数，及び，位相を変数とする波生成式の数を項数とした場合に、当該項数の総和で示された波生成式であって、前記第４のステップは、前記２枚の画像を複数のサブブロックにそれぞれ分割する第６のステップと、当該２枚の画像において同じ場所に位置するサブブロックの強度値を前記目的関数に代入して、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する第７のステップと、所定の前記項数の範囲で前記第７のステップの推定を繰り返し、推定された前記周波数に基づいて波のスペクトル分布を求め、当該波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する第８のステップと、を有することを要旨とする。

第１１の本発明に係る水面波挙動検出方法は、前記第８のステップが、前記波のスペクトル分布の面積を１になるように正規化し、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布の速度ベクトルを推定することを要旨とする。

第１２の本発明に係る水面波挙動検出方法は、前記第８のステップが、前記速度拘束値及び／又は前記サブブロックの大きさを更に変化させて前記第７のステップの推定を繰り返し、最も類似度の高いスペクトル分布の速度ベクトルを推定することを要旨とする。

第１３の本発明に係る水面波挙動検出方法は、前記所定のスペクトル分布が、ブレットシュナイダーのスペクトル分布であることを要旨とする。

第１４の本発明に係る水面波挙動検出方法は、前記波生成式が、振幅を前記項数で除算し、波数，周波数，及び，位相を前記項数で乗算した波生成式であることを要旨とする。

第１５の本発明に係る水面波挙動検出方法は、前記第４のステップが、共役勾配法を用いて前記速度ベクトルを推定することを要旨とする。

第１６の本発明に係る水面波挙動検出方法は、前記第３のステップが、前記ロバスト関数に対してローレンツ関数を用いることを要旨とする。

第１７の本発明に係る水面波挙動検出プログラムは、請求項９乃至１６のいずれか１項に記載の水面波挙動検出方法における各ステップをコンピュータによって実行させることを要旨とする。

本発明によれば、水面波を捉えた画像から水面波の挙動を検出することができる。

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

〔水面波挙動検出装置の構成及び処理の流れについて〕
図１は、本発明の実施の形態に係る水面波挙動検出装置１０を示した構成図である。水面波挙動検出装置１０は、水面波を含む画像を入力する画像入力部１００と、入力された画像を時系列に蓄積する画像蓄積部１１０と、目的関数を導出する目的関数導出部１２０と、速度ベクトルを推定するオプティカルフロー推定部１３０と、推定された速度ベクトルを表示する表示部１４０とを備えた構成である。

なお、水面波挙動検出装置１０は、コンピュータ本体を構成するＣＰＵ，メモリ，ハードディスクなどにより実現され、以下に説明する各ステップの処理はプログラムで実行される。

続いて、水面波挙動検出装置１０を構成する各部の処理の流れについて説明する。なお、画像の輝度値，強度値または濃淡値の意味するところは全て等価であり、計算する際に使用される画像データは画像の強度値として以下説明する。

最初に、画像入力部１００には、水面波が記録された複数の画像が入力される（ステップＳ１）。この水面波とは、例えば、海上の波や、川の波などであり、入力される画像とは、例えば、図２（ａ）で示すような画像である。

次に、画像蓄積部１１０は、画像入力部１００で入力された複数の画像を時系列に蓄積する（ステップＳ２）。時系列とは、時間とともに変動する現象の観測値を時間とともに記録した並びを意味するものである。

続いて、目的関数導出部１２０は、波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、導出した輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出する（ステップＳ３）。

なお、上記波生成式とは、波の振幅，波数，周波数，及び，位相を変数とする波生成式の数を項数とした場合に、この項数の総和で示された波生成式である。目的関数の具体的な導出方法については後述する。

そして、オプティカルフロー推定部１３０は、連続する２枚の画像を画像蓄積部１１０から読み出して、読み出した２枚の画像の強度値をステップＳ３で導出した目的関数に代入し、速度ベクトルを推定する（ステップＳ４）。

オプティカルフロー推定部１３０は、画像をサブブロックに分割する分割部１３１、速度ベクトルを局所的に推定する第１の推定部１３２、速度ベクトルを大局的に推定する第２の推定部１３３を備えている。

まず、分割部１３１は、画像蓄積部１１０から読み出した連続する２枚の画像を、図２（ｂ）に示すように、複数のサブブロックにそれぞれ分割する。

次に、第１の推定部１３２は、分割された２枚の画像において、同じ場所に位置するサブブロックの強度値を目的関数に代入して、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する。

そして、第２の推定部１３３は、所定の項数の範囲で第１の推定部１３２の推定を繰り返し、推定された周波数に基づいて波のスペクトル分布を求め、求めた波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する。

なお、第２の推定部１３３は、所定の項数の範囲のみではなく、後述する速度拘束値やサブブロックの大きさを更に変化させて、その変化させた範囲内で最も類似度の高い速度ベクトルを推定することも可能である。オプティカルフロー推定部１３０に関するより詳細な処理の流れについては後述する。

最後に、表示部１４０は、ステップ４で推定された速度ベクトルを、水面波の動きベクトルとして矢印記号で画面に表示する（ステップＳ５）。

〔目的関数の導出方法について〕
次に、ステップＳ３で説明した目的関数の具体的な導出方法について、その理由を踏まえながら説明する。

最初に、波生成モデルについて説明する。波生成式は、流体力学や海洋工学において数多くモデル化されている不規則で多方向に変化しながら動く波の基本特性を記述する式であって、多方向性および不規則性の波生成（ＭＩＷＧ：Multidirectionality-Irregularity-Wave-Generation）式を適用する。このＭＩＷＧ方程式は、線形理論に基づいているので、他のモデルと組み合わせを行う上で数理的に都合がよく、実際に、港湾設計や沖合いの波浪のシュミレーションなどに用いられている。

ＭＩＧＷ方程式は、振幅ａ_ｍ，波数ｋ_ｍ，周波数ｆ_ｍ，及び，位相θ_ｍを用いて、式（１）で表現することができる。

式（１）の左辺であるηは、時間ｔにおける画像の２次元位置（ｘ，ｙ）の波の高さである。また、右辺のＭは、いわゆる波らしさを表現するものであり、右辺のΣを除いた波生成式の総和の項数である。

まず、式（１）についての項数を決定する必要がある。通常、波生成のシュミレーションにおける項数は１００前後又はそれ以上である。しかしながら、本実施の形態では、（ｉ）低解像度の画像を用いており、（ｉｉ）波の統計的な特性として利用する後述のＢ−スペクトラム分布は単一のピークを持ち、そのピークは低周波数側に存在し、（ｉｉｉ）波における相対的な振幅は、高周波数，高い波数ほど小さいことが知られているため、３〜１０の項数で十分に波を近似することができる。項数を大幅に減らすことで計算効率の向上することができる。

ここで、特に上記（ｉｉｉ）の理由から、式（１）を式（２）で再定義（近似）する。

式（２）は、式（１）の波生成式において、振幅ａ_ｍを項数ｍで除算し、且つ、波数ｋ_ｍ，周波数ｆ_ｍ，及び，位相θ_ｍを項数ｍで乗算したものである。これにより、例えば、ｃｏｓθがｃｏｓ４θになるので、項数ｍが大きくなると、高周波数の波が現れることになる。即ち、式（１）では、高周波の項が得られるまでランダムに多くの項数を用いる必要があるが、式（２）を用いることで、高い周波数の波が少数項でも確実に得ることができる。

次に、光学反射モデルについて説明する。画像間の水面変化に伴い、透明，白波など、多くの要因が輝度変化を生じさせている。特に、太陽光と水面との間の鏡面反射と陰影との影響を考慮した場合、この光学反射は局所的に急激な変化をもたらすことになる。故に、式（３）で示す水面の輝度の指数関数を更に適用する必要がある。

なお、κは指数係数、Iは２次元画像、δｔはサンプリング時間である。

続いて、画像から速度ベクトルを推定するため、輝度変動の拘束条件として、式（４）で示す輝度拘束条件式を用いて、上記にて説明した波生成モデルと光学反射モデルとを統合する。

なお、∇は空間１次微分、ｕはオプティカルフローの速度ベクトル、Ｉ_ｔは時間１次微分である。また、ｑは物理現象の数を意味している。

そして、式（２）〜式（４）に基づいて、式（２）の時間差分を計算して、式（５）に示す輝度変動方程式を導出する。

なお、（ｕ，ｖ）は、オプティカルフローの２次元成分であるｘ軸方向およびｙ軸方向の速度ベクトルである。

式（５）を用いて速度ベクトル，振幅，波数，周波数，及び，位相を推定する際に、項数Ｍに加えて、波の物理的な性質を加味する必要がある。特に、波は、波数ｋと角周波数ωとの関係を、式（６）で示す分散関係式で表現できることが知られている。

なお、ω＝２πｆ，ｋ，ｈ，ｇは、角振動数，波数，水深，重力加速度である。また、ｆは周波数，Ｔ＝１／ｆである。式（６）は、波長Ｌ＝２π／ｋ毎で、波の速度ｕ（＝Ｌ／Ｔ＝ω／ｋ）が異なることを示している。

また、波の時系列画像には、その他、様々な複雑な動き成分が含まれている。そこで、流体と波との基本的な特性に基づき、速度に対して拘束条件を課すことで、式（５）の計算を安定化させることができる。具体的には、波の粘性の影響により、波の高周波成分が急速に減衰する傾向があるので、式（７）で示す速度と周波数との反比例関係の速度拘束条件式を更に用いて、式（５）の計算を行う。

なお、ｕは速度、γは速度拘束値、ｆは周波数である。

以上説明した式（４）〜式（７）に基づいて、オプティカルフローの速度ベクトル（ｕ，ｖ）と、５つの物理パラメータ（振幅ａ_ｍ，波数ｋ_ｍ，周波数ｆ_ｍ，位相θ_ｍ，指数係数κ）とを推定するため、式（８）に示すように、輝度変動項（ｉｍｇ＿ｖａｒ）と速度拘束条件（ｖｅｌ＿ｃｏｎｓｔｒａｉｎｔ）の２つの項からなる目的関数を導出する。

なお、ρはロバスト統計学で用いられているロバスト関数である。また、σは分散である。ロバスト関数ρには、ローレンツ関数∂ρ／∂ｅを用いることも可能であり、これにより目的関数の計算速度を向上することができる。

最後に、式（８）で示された目的関数を用いて速度ベクトルを推定する際の最小化計算について説明する。この目的関数は非線形方程式なので、最急勾配法（ＳＧ法）や、Ｆｌｅｔｃｈｅｒ−Ｒｅｅｖｅｓの共役勾配法（ＣＧ法）を用いて最小化計算をすることができる。

ＣＧ法を用いて計算する際、速度ベクトルと５つの物理パラメータとに対する偏微分方程式を導出し、時間ｔと時間ｔ＋１との２枚の時系列画像の強度値をそれぞれ入力して、反復回数をｐとする式（９）を用いて反復計算を繰り返す。

なお、ＣＧ法を用いて計算する方法は、ＳＧ法よりも少ない数回で収束し、極小値が小さい値となる利点を持つので、目的関数の最小化をより効率的に計算することができる。

〔オプティカルフロー推定部の処理の流れについて〕
図３は、オプティカルフロー推定部１３０の処理の流れを示すフロー図である。式（９）を用いて最小化計算を行う際に、項数Ｍ，速度拘束値γ，サブブロックの大きさのそれぞれを決定する必要がある。項数Ｍの大きさを変えることは波の高周波数を変えることになり、サブブロックの大きさを変えることは動き（フロー）の平滑化に関係することになる。また、速度拘束値γについては、フローの急な変化（不連続性）に関係することになる。

最初に、オプティカルフロー推定部１３０は、式（９）に項数Ｍ＝１，速度拘束値γ＝０．０の初期値を代入する（ステップＳ４−１）。

次に、分割部１３１において、画像蓄積部１１０から読み出した連続する２枚の画像を、例えば、３×３の大きさのサブブロックにそれぞれ分割する（ステップＳ４−２）。

続いて、第１の推定部１３２において、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する（ステップＳ４−３）。

ステップＳ４−３により、各サブブロックにおける局所的な速度ベクトルが推定されるので、サブブロックごとに細かい波の動きを推定することができる。しかしながら、実際の環境における波の現象は複雑であり、隣接するサブブロック同士の波の連続性を把握できない可能性がある。

そこで、波の統計分布として知られたブレットシュナイダー（Ｂｒｅｔｓｃｈｎｅｉｄｅｒ）のスペクトル分布（以降、「Ｂ−スペクトル分布」と称する）と比較して類似度を検討し、全体としての動きを大局的に捉える必要がある。Ｂ−スペクトル分布は式（１０）で示され、単一のピークを持つガンマ関数に近い分布を有している。

なお、有義波とは、波群中、波高の大きい方から数えて、全体の波の数の１／３の数の波を選び出したものである。

即ち、第２の推定部１３３において、ステップＳ４−３で推定された周波数の値を式（１０）に代入し、図２（ｃ）に示すように全てのサブブロックにおける波のスペクトル分布を求め、求めた波のスペクトル分布をＥｓｔ（ｆ）_{Ｍ，γ，ｓｕｂｂｌｏｃｋ}とした場合に、求めた波のスペクトル分布とＢ−スペクトル分布との類似度を式（１１）を用いて計算する。そして、例えば、項数Ｍを１〜１０の範囲、速度拘束値γを０．０〜３．０の範囲、サブブロックの大きさを３×３〜２０×２０の範囲でそれぞれ変化させて同様に類似度を計算し、式（１１）の値が最小となる項数Ｍ，速度拘束値γ，サブブロックの大きさのスペクトル分布を最も類似度が高いスペクトル分布と判定し、そのスペクトル分布の速度ベクトルを推定する。（ステップＳ４−４）。

これにより、画像全体の連続性を考慮した水面波の挙動を検出することが可能となる。

更に、Ｂ−スペクトル分布は、前述したように、単一のピークを持ち、そのピークは低周波数側に存在するという特性を利用して、次に説明する方法で制限された範囲で速度ベクトルを推定することも可能である。

まず、第２の推定部１３３は、上記ステップＳ４−３で推定された周波数の値を用いて全てのサブブロックのスペクトル分布を求め、求めた波のスペクトル分布の面積を１になるように正規化する。

その後、第２の推定部１３３は、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の第１角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の第２角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布を有する速度ベクトルを推定する。

なお、上記第１角度，第２角度とは、Ｂ−スペクトル分布におけるピーク値に対応するＸ軸の値を基準値とした場合に、その基準値からＸ軸方向に所定の距離だけ離れた場所におけるＢ−スペクトル曲線の接線と、Ｘ軸とのなす角度である。ここでは角度を例に説明したが、上記基準値からの所定の距離を有効範囲とするものであってもよい。

〔他の速度拘束条件式について〕
式（７）を用いて説明した速度拘束条件式の他に、例えば、式（１２）で示す拘束条件を用いることもできる。即ち、周波数ｆ_ｍが所定の閾値としての１．０より大きい場合には、速度ベクトル（ｕ，ｖ）をＷｅｉｇｈｔで除算することで小さくするように重み付けを行い、１．０より小さい場合には、速度ベクトルをＷｅｉｇｈｔで乗算することで大きくするように重み付けを行う。

〔推定結果について〕
図４は、３種類の波の画像を示す図である。図４（ａ）は、弱い波の画像であり、図４（ｂ）は、穏やかな波の画像であり、図４（ｃ）は、強い波の画像である。これらの画像は、３次元世界が２次元画像座標系へ投影された画像であり、水面に対して斜め上方向からの画像は、波の山と谷との検出を可能とする。

図５は、図４で示す３種類の波の画像を第１の従来法で推定した速度ベクトルを矢印で表示した図である。この第１の従来法は、流体力学で用いる連続方程式を輝度変動方程式の代わりに適用し、速度拘束条件式に発散（ｄｉｖｅｒｇｅｎｃｅ）と回転（ｃｕｒｌ）を適用した方法である。

図６は、図４で示す３種類の波の画像を第２の従来法で推定した速度ベクトルを矢印で表示した図である。この第２の従来法は、流体力学で用いる連続方程式のみを適用した方法である。

図７は、図４で示す３種類の波の画像を本実施の形態で推定した速度ベクトルを矢印で表示した図である。また、図８は、図４（ｃ）で示す強い波の場合における本実施の形態で推定した指数係数，波数，周波数，位相，振幅の２次元分布を示した図である。更に、図９は、図４で示す３種類の波の画像に対する本実施の形態で推定した振幅の２次元分布を示した図である。

図５及び図６で示された全ての波は、全体が滑らかな動きで表現されており、波の山や谷といった特徴を把握することができない。一方、図７で示された水面波は、いずれの波についても、上下運動を中心に、はっきりとした波の動きを確認することができる。また、図８では、強い波の特徴が明確に把握することはできる。更に、図９では、弱い波，穏やかな波，強い波の順に、波の特徴が表現されていることがはっきりと確認できる。

図１０（ａ）は、Ｂ−スペクトル分布を示す分布図である。図１０（ａ）では、自然界における波のスペクトル分布が示されており、低周波数側に一つのピークが存在し、高周波数側に広がりを有することを特徴としている。また、図１０（ｂ）で示す分布図は、異なる項数における強い波のスペクトル分布図であって、全てのサブブロックの周波数を面積が１になるように正規化した結果である。実線２１０で示す分布の方が、破線２２０で示す分布よりもＢ−スペクトル分布に類似していることが確認することができる。また、実線２１０および破線２２０で示す分布は、共に、低周波数側にピークがあり、更に、１つのピークを有する分布図であるが、実線２１０の方が、ピークから低周波数側の分布特性が急な角度で下降し、且つ、ピークから高周波数側の分布特性が緩やかな角度で下降しているので、図１０（ｂ）の場合には、実線２１０の項数を有する速度ベクトルを用いた水面波が検出されることになる。

本実施の形態によれば、波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、この輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出し、連続する２枚の画像の強度値をこの目的関数に代入して速度ベクトルを推定するので、輝度変動に基づくオプティカルフロー法を用いて、水面波を捉えた画像から不規則な水面波の挙動を検出することができる。

本実施の形態によれば、２枚の画像を複数のサブブロックにそれぞれ分割し、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する第１の推定部１３２と、所定の項数の範囲で第１の推定部１３２の推定を繰り返し、推定された周波数に基づいて波のスペクトル分布を求め、この波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する第２の推定部１３３とを有するので、水面波の挙動をより正確に検出することができる。即ち、第１の推定部１３２で推定された局所的な速度ベクトルは、隣接するサブブロック間の連続性を考慮していないが、第２の推定部１３３において、推定された周波数から波のスペクトル分布を求め、この波のスペクトル分布と所定のスペクトル分布との類似度を判定し、所定の項数の範囲で最も類似度の高いスペクトル分布の速度ベクトル推定するため、画像全体の連続性を考慮した水面波の挙動を検出することが可能となる。

本実施の形態によれば、波のスペクトル分布の面積を１になるように正規化し、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布の速度ベクトルを推定するので、より正確に水面波の挙動を検出することができる。

本実施の形態によれば、項数に加えて、速度拘束値及び／又はサブブロックの大きさを更に変化させて速度ベクトルなどを推定し、類似度の高いスペクトル分布の速度ベクトルを推定するので、より正確に水面波の挙動を検出することができる。

本実施の形態によれば、波の統計分布として信頼性の高いブレットシュナイダーのスペクトル分布を用いて類似度を判定するので、信頼性の高い水面波の挙動を検出することができる。

本実施の形態によれば、振幅を項数で除算し、波数，周波数，及び，位相を項数で乗算した波生成式であるので、高い周波数の波を少数項でも確実に得ることができる。

本実施の形態によれば、共役勾配法を用いて速度ベクトルを推定するので、目的関数の最小化をより効率的に計算することができる。

本実施の形態によれば、ロバスト関数に対してローレンツ関数を用いるので、目的関数の計算速度を向上することができる。

本発明の実施の形態に係る水面波挙動検出装置を示した構成図である。 本発明の実施の形態に係る水面波挙動検出装置の処理の流れを示すフロー図である。 オプティカルフロー推定部の処理の流れを示すフロー図である。 ３種類の波の画像を示す図である。 図４で示す３種類の波の画像を第１の従来法で推定した速度ベクトルを矢印で表示した図である。 図４で示す３種類の波の画像を第２の従来法で推定した速度ベクトルを矢印で表示した図である。 図４で示す３種類の波の画像を本実施の形態で推定した速度ベクトルを矢印で表示した図である。 図４（ｃ）で示す強い波の場合における本実施の形態で推定した指数係数，波数，周波数，位相，振幅の２次元分布を示した図である。 図４で示す３種類の波の画像に対する本実施の形態で推定した振幅の２次元分布を示した図である。 スペクトル分布を示す分布図である。

符号の説明

１０…動きベクトル検出装置
１００…画像入力部
１１０…画像蓄積部
１２０…目的関数導出部
１３０…オプティカルフロー推定部
１３１…分割部
１３２…第１の推定部
１３３…第２の推定部
１４０…表示部
２１０…実線
２２０…破線
Ｓ１〜Ｓ４，Ｓ４−１〜Ｓ４−４…ステップ

Claims (17)

1. 水面波を含む画像を入力する画像入力手段と、
   入力された前記画像を時系列に蓄積する画像蓄積手段と、
   波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、当該輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出する目的関数導出手段と、
   連続する２枚の画像を前記画像蓄積手段から読み出して、当該２枚の画像の強度値を前記目的関数に代入し、前記速度ベクトルを推定するオプティカルフロー推定手段と、
   推定された当該速度ベクトルを表示する表示手段と、
   を有することを特徴とする水面波挙動検出装置。
2. 前記波生成式は、波の振幅，波数，周波数，及び，位相を変数とする波生成式の数を項数とした場合に、当該項数の総和で示された波生成式であって、
   前記オプティカルフロー推定手段は、
   前記２枚の画像を複数のサブブロックにそれぞれ分割する分割手段と、
   当該２枚の画像において同じ場所に位置するサブブロックの強度値を前記目的関数に代入して、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する第１の推定手段と、
   所定の前記項数の範囲で前記第１の推定手段の推定を繰り返し、推定された前記周波数に基づいて波のスペクトル分布を求め、当該波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する第２の推定手段と、
   を有することを特徴とする請求項１に記載の水面波挙動検出装置。
3. 前記第２の推定手段は、前記波のスペクトル分布の面積を１になるように正規化し、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布の速度ベクトルを推定することを特徴とする請求項２に記載の水面波挙動検出装置。
4. 前記第２の推定手段は、前記速度拘束値及び／又は前記サブブロックの大きさを更に変化させて前記第１の推定手段の推定を繰り返し、最も類似度の高いスペクトル分布の速度ベクトルを推定することを特徴とする請求項２又は３に記載の水面波挙動検出装置。
5. 前記所定のスペクトル分布は、ブレットシュナイダーのスペクトル分布であることを特徴とする請求項２乃至４のいずれか１項に記載の水面波挙動検出装置。
6. 前記波生成式は、振幅を前記項数で除算し、波数，周波数，及び，位相を前記項数で乗算した波生成式であることを特徴とする請求項２乃至５のいずれか１項に記載の水面波挙動検出装置。
7. 前記オプティカルフロー推定手段は、共役勾配法を用いて前記速度ベクトルを推定することを特徴とする請求項１乃至６のいずれか１項に記載の水面波挙動検出装置。
8. 前記目的関数導出手段は、前記ロバスト関数に対してローレンツ関数を用いることを特徴とする請求項１乃至７のいずれか１項に記載の水面波挙動検出装置。
9. 水面波を含む画像を入力する第１のステップと、
   入力された前記画像を時系列に画像蓄積手段に蓄積する第２のステップと、
   波の不規則性を示す波生成式と水面の輝度を示す指数関数とに基づいて速度ベクトルを変数とする輝度変動方程式を導出し、当該輝度変動方程式と波数及び角周波数で表された分散関係式と所定の速度拘束値で波の速度を拘束する速度拘束条件式とに基づいて、ロバスト統計学で用いるロバスト関数を介した目的関数を導出する第３のステップと、
   連続する２枚の画像を前記画像蓄積手段から読み出して、当該２枚の画像の強度値を前記目的関数に代入し、前記速度ベクトルを推定する第４のステップと、
   推定された当該速度ベクトルを表示する第５のステップと、
   を有することを特徴とする水面波挙動検出方法。
10. 前記波生成式は、波の振幅，波数，周波数，及び，位相を変数とする波生成式の数を項数とした場合に、当該項数の総和で示された波生成式であって、
    前記第４のステップは、
    前記２枚の画像を複数のサブブロックにそれぞれ分割する第６のステップと、
    当該２枚の画像において同じ場所に位置するサブブロックの強度値を前記目的関数に代入して、全てのサブブロックに対する速度ベクトル，振幅，波数，周波数，及び，位相を推定する第７のステップと、
    所定の前記項数の範囲で前記第７のステップの推定を繰り返し、推定された前記周波数に基づいて波のスペクトル分布を求め、当該波のスペクトル分布と所定のスペクトル分布との類似度を判定し、最も類似度の高いスペクトル分布の速度ベクトルを推定する第８のステップと、
    を有することを特徴とする請求項９に記載の水面波挙動検出方法。
11. 前記第８のステップは、前記波のスペクトル分布の面積を１になるように正規化し、正規化された全てのスペクトル分布において、低周波数側にピークがあるスペクトル分布のみを残し、且つ、複数のピークを有するスペクトル分布を排除したスペクトル分布のうち、ピークから低周波数側の分布特性が所定の角度よりも急な角度で下降し、且つ、ピークから高周波数側の分布特性が所定の角度よりも緩やかな角度で下降しているスペクトル分布に該当するスペクトル分布のうち、最も類似度の高いスペクトル分布の速度ベクトルを推定することを特徴とする請求項１０に記載の水面波挙動検出方法。
12. 前記第８のステップは、前記速度拘束値及び／又は前記サブブロックの大きさを更に変化させて前記第７のステップの推定を繰り返し、最も類似度の高いスペクトル分布の速度ベクトルを推定することを特徴とする請求項１０又は１１に記載の水面波挙動検出方法。
13. 前記所定のスペクトル分布は、ブレットシュナイダーのスペクトル分布であることを特徴とする請求項１０乃至１２のいずれか１項に記載の水面波挙動検出方法。
14. 前記波生成式は、振幅を前記項数で除算し、波数，周波数，及び，位相を前記項数で乗算した波生成式であることを特徴とする請求項１０乃至１３のいずれか１項に記載の水面波挙動検出方法。
15. 前記第４のステップは、共役勾配法を用いて前記速度ベクトルを推定することを特徴とする請求項９乃至１４のいずれか１項に記載の水面波挙動検出方法。
16. 前記第３のステップは、前記ロバスト関数に対してローレンツ関数を用いることを特徴とする請求項９乃至１５のいずれか１項に記載の水面波挙動検出方法。
17. 請求項９乃至１６のいずれか１項に記載の水面波挙動検出方法における各ステップをコンピュータによって実行させることを特徴とする水面波挙動検出プログラム。