JP7069433B1 - 風速予測装置、風速予測方法及びレーダ装置 - Google Patents

Abstract

視線方向と水平方向とのなす角である仰角が互いに異なる複数のビームのそれぞれが空間に放射され、空間で散乱された後のそれぞれのビームである散乱信号を取得する散乱信号取得部（１１）と、空間に放射されたそれぞれのビームの仰角に従って、散乱信号取得部（１１）により取得されたそれぞれの散乱信号の距離分解能としてレンジビン幅を設定し、それぞれの散乱信号から、観測領域を含んでいる２次元平面に対応するレンジビンのドップラー周波数を算出するドップラー周波数算出部（１２）と、ドップラー周波数算出部（１２）により算出された複数のドップラー周波数から、ＶＶＰ法を用いて、２次元平面の風速分布を推定する第１の風速分布推定部（１３）と、第１の風速分布推定部（１３）により推定された２次元平面の風速分布から、２次元ナビエストークス方程式を用いて、観測領域の風速を予測する風速予測部（１５）とを備えるように、風速予測装置（３）を構成した。

Description

本開示は、風速予測装置、風速予測方法及びレーダ装置に関するものである。

観測領域の風速を予測する風速予測装置の中には、気象モデルを用いて風速のシミュレーションを実施することで、観測領域の風速を予測する風速予測装置（以下「従来の風速予測装置」という）がある（特許文献１を参照）。風速のシミュレーションは、３次元の演算処理である。

特開２０１４－２０２１９０号公報

特許文献１に開示されている従来の風速予測装置は、３次元の演算処理を行うため、演算負荷が高くなってしまうことがあるという課題があった。

本開示は、上記のような課題を解決するためになされたもので、３次元の演算処理を実施することなく、風速を予測することができる風速予測装置及び風速予測方法を得ることを目的とする。

本開示に係る風速予測装置は、視線方向と水平方向とのなす角である仰角が互いに異なる複数のビームのそれぞれが空間に放射され、空間で散乱された後のそれぞれのビームである散乱信号を取得する散乱信号取得部と、空間に放射されたそれぞれのビームの仰角に従って、散乱信号取得部により取得されたそれぞれの散乱信号の距離分解能としてレンジビン幅を可変に設定することで、仰角の異なるビームであっても同一高度の観測点を観測できるようにし、それぞれの散乱信号から、観測領域を含んでいる２次元平面に対応するレンジビンのドップラー周波数を算出するドップラー周波数算出部と、ドップラー周波数算出部により算出された複数のドップラー周波数から、ＶＶＰ（Ｖｏｌｕｍｅ Ｖｅｌｏｃｉｔｙ Ｐｒｏｃｅｓｓｉｎｇ）法を用いて、２次元平面の風速分布を推定する第１の風速分布推定部と、第１の風速分布推定部により推定された２次元平面の風速分布から、２次元ナビエストークス方程式を用いて、観測領域の風速を予測する風速予測部とを備えたものである。

本開示によれば、３次元の演算処理を実施することなく、風速を予測することができる。

実施の形態１に係る風速予測装置３を含むレーダ装置１を示す構成図である。 図２Ａは、飛行領域の位置が、ビーム送受信部２の鉛直方向の真上である場合の複数のビームを示す説明図、図２Ｂは、飛行領域の位置が、ビーム送受信部２の鉛直方向の真上からずれている場合の複数のビームを示す説明図である。 実施の形態１に係る風速予測装置３を示す構成図である。 実施の形態１に係る風速予測装置３のハードウェアを示すハードウェア構成図である。 風速予測装置３が、ソフトウェア又はファームウェア等によって実現される場合のコンピュータのハードウェア構成図である。 風速予測装置３の処理手順である風速予測方法を示すフローチャートである。 図７Ａは、仰角θ_ｎが異なれば、同一のレンジビンＲｂ_ｍの複数の観測点の高度が互いに異なることを示す説明図、図７Ｂは、仰角θ_ｎが変化しても、同一のレンジビンＲｂ_ｍの複数の観測点が同じ高度であることを示す説明図である。 風速の観測点の位置（ｘ，ｙ，ｚ）を示す説明図である。 実施の形態２に係る風速予測装置３を示す構成図である。 実施の形態２に係る風速予測装置３のハードウェアを示すハードウェア構成図である。 実施の形態３に係る風速予測装置３を示す構成図である。 実施の形態３に係る風速予測装置３のハードウェアを示すハードウェア構成図である。

以下、本開示をより詳細に説明するために、本開示を実施するための形態について、添付の図面に従って説明する。

実施の形態１．
図１は、実施の形態１に係る風速予測装置３を含むレーダ装置１を示す構成図である。
図１に示すレーダ装置１は、ビーム送受信部２及び風速予測装置３を備えている。
図２は、ビーム送受信部２から空間に放射される複数のビームを示す説明図である。図２の例では、観測領域が空間に存在している。観測領域は、例えば、飛行体が飛行する予定の飛行領域である。飛行体は、例えば、ドローン、又は、ヘリコプターである。
図２Ａは、飛行領域の位置が、ビーム送受信部２の鉛直方向の真上である場合の複数のビームを示す説明図である。
図２Ｂは、飛行領域の位置が、ビーム送受信部２の鉛直方向の真上からずれている場合の複数のビームを示す説明図である。

ビーム送受信部２は、例えば、ビーム生成部、ビーム送信器、放射方向切替部、アンテナ及びビーム受信器を備えている。図１では、ビーム生成部、ビーム送信器、放射方向切替部、アンテナ及びビーム受信器の表記が省略されている。
ビーム送受信部２は、図２に示すように、視線方向と水平方向とのなす角である仰角θ_ｎ（ｎ＝１，・・・，Ｎ）が互いに異なる複数のビームのそれぞれを空間に放射する。Ｎは、２以上の整数である。
即ち、ビーム送受信部２は、図２Ａに示すように、飛行領域が鉛直方向の真上に存在している場合、それぞれの仰角θ_ｎでビームスキャンを実施することで、それぞれの仰角θ_ｎにおいて、互いに異なる方向にビームを放射する。
図２Ａでは、Ｎ＝３であり、ビーム送受信部２から、仰角θ_１のビーム、仰角θ_２のビーム及び仰角θ_３のビームのそれぞれが空間に放射されている。図２Ａでは、それぞれの仰角θ_ｎにおいて、４つの方向にビームが放射されている。図２Ａにおいて、黒丸は、風速の観測点を示している。
また、ビーム送受信部２は、図２Ｂに示すように、飛行領域の位置が鉛直方向の真上からずれている場合、仰角θ_ｎを切り替えながら、それぞれの仰角θ_ｎでビームを放射する。
図２Ｂでは、Ｎ＝１２であり、ビーム送受信部２から、仰角θ_１～θ_１２のビームが空間に放射されている。図面の簡単化のため、図２Ｂには、仰角θ_１～θ_６のみが表記されており、仰角θ_７～θ_１２の表記が省略されている。図２Ｂでは、それぞれの仰角θ_ｎにおいて、１つの方向にビームが放射されている。図２Ｂにおいて、黒丸は、風速の観測点を示している。

図２Ａ及び図２Ｂのそれぞれは、３つの２次元平面の風速分布を求める例を示している。しかし、これは一例に過ぎず、１つ又は２つの２次元平面の風速分布を求めるものであってもよいし、４つ以上の２次元平面の風速分布を求めるものであってもよい。
ビーム送受信部２から放射されるビームとしては、例えば、ＣＷ（Ｃｏｎｔｉｎｕｏｕｓ Ｗａｖｅ）パルス光のほか、周波数変調されたパルス光を用いることができる。
ビーム送受信部２は、空間で散乱された後のそれぞれのビームを散乱信号として受信する。
ビーム送受信部２は、それぞれの散乱信号をアナログ信号からデジタル信号に変換し、デジタル信号を風速予測装置３に出力する。
また、ビーム送受信部２は、空間に放射したそれぞれのビームの仰角θ_ｎを示す角度情報を風速予測装置３に出力する。
図１に示すレーダ装置１では、ビーム送受信部２が、ビーム生成部、ビーム送信器、放射方向切替部、アンテナ及びビーム受信器を備えているものを想定している。しかし、ビーム送受信部２は、仰角θ_ｎが互いに異なる複数のビームのそれぞれを空間に放射し、空間で散乱された後のそれぞれのビームを散乱信号として受信することができればよく、どのような構成であってもよい。

図３は、実施の形態１に係る風速予測装置３を示す構成図である。
図４は、実施の形態１に係る風速予測装置３のハードウェアを示すハードウェア構成図である。
図３に示す風速予測装置３は、散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４及び風速予測部１５を備えている。

散乱信号取得部１１は、例えば、図４に示す散乱信号取得回路２１によって実現される。
散乱信号取得部１１は、空間で散乱された後のそれぞれの散乱信号として、ビーム送受信部２から、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）を取得する。ｊ＝１，・・・，Ｊである。ｊは、仰角θ_ｎにおけるビームスキャンのスキャン番号である。図２Ａの例では、それぞれの仰角θ_ｎにおいて、４つの方向にビームが放射されているため、Ｊ＝４である。図２Ｂの例では、それぞれの仰角θ_ｎにおいて、１つの方向にビームが放射されているため、Ｊ＝１である。ｔは、時刻を示す変数である。
また、散乱信号取得部１１は、ビーム送受信部２により空間に放射されたそれぞれのビームの仰角θ_ｎを示す角度情報を取得する。
散乱信号取得部１１は、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）及び角度情報のそれぞれをドップラー周波数算出部１２に出力する。

ドップラー周波数算出部１２は、例えば、図４に示すドップラー周波数算出回路２２によって実現される。
ドップラー周波数算出部１２は、散乱信号取得部１１から、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）及び角度情報のそれぞれを取得する。
ドップラー周波数算出部１２は、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）をヒット方向にＦＦＴ（Ｆａｓｔ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍａｔｉｏｎ）する前に、角度情報が示す仰角θ_ｎに従って、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）の距離分解能としてレンジビン幅を設定する。即ち、ドップラー周波数算出部１２は、レンジビン幅に相当するＦＦＴ長Ｌ_θｎを設定する。
ＦＦＴ長Ｌ_θｎは、空間に放射されたビームの仰角θ_ｎが広いほど、短く設定され、空間に放射されたビームの仰角θ_ｎが狭いほど、長く設定される。
図２Ａの例では、仰角θ_１のビームに対応するデジタル信号ｓ_θ１，ｊ（ｔ）のＦＦＴ長Ｌ_θ１が、仰角θ_２のビームに対応するデジタル信号ｓ_θ２，ｊ（ｔ）のＦＦＴ長Ｌ_θ２よりも長く、仰角θ_２のビームに対応するデジタル信号ｓ_θ２，ｊ（ｔ）のＦＦＴ長Ｌ_θ２が、仰角θ_３のビームに対応するデジタル信号ｓ_θ３，ｊ（ｔ）のＦＦＴ長Ｌ_θ３よりも長い。
ドップラー周波数算出部１２は、散乱信号取得部１１により取得されたそれぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）をヒット方向にＦＦＴすることで、デジタル信号ｓ_θｎ，ｊ（ｔ）を周波数領域の信号ｒ_θｎ，ｊ（ｆ）に変換する。ｆは、周波数を示す変数である。

ドップラー周波数算出部１２は、周波数領域の信号ｒ_θｎ，ｊ（ｆ）から、複数のレンジビンＲｂ_１～Ｒｂ_Ｍにおけるそれぞれのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を算出する。ｍ＝１，・・・，Ｍであり、Ｍは、１以上の整数である。
図２では、３つの２次元平面の風速分布を求める例を示している。図２のように、３つの２次元平面の風速分布を求める場合、ドップラー周波数算出部１２は、周波数領域の信号ｒ_θｎ，ｊ（ｆ）から、レンジビンＲｂ_１，Ｒｂ_２，Ｒｂ_３におけるそれぞれのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を算出する。
図３に示す風速予測装置３では、Ｍが２以上の整数であるものを想定している。ただし、Ｍが２以上の整数であるものに限るものではなく、Ｍ＝１であってもよい。
Ｍ＝１の場合、ドップラー周波数算出部１２は、周波数領域の信号ｒ_θｎ，ｊ（ｆ）から、観測領域を含んでいる１つのレンジビンＲｂ_１のドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂ１}のみを算出する。
ドップラー周波数算出部１２は、それぞれのレンジビンＲｂ_ｍのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を第１の風速分布推定部１３に出力する。

第１の風速分布推定部１３は、例えば、図４に示す第１の風速分布推定回路２３によって実現される。
第１の風速分布推定部１３は、ドップラー周波数算出部１２により算出されたそれぞれのレンジビンＲｂ_ｍにおける複数のドップラー周波数ｄｐｆ_{θ１，ｊ，Ｒｂｍ}～ｄｐｆ_{θＮ，ｊ，Ｒｂｍ}から、ＶＶＰ（Ｖｏｌｕｍｅ Ｖｅｌｏｃｉｔｙ Ｐｒｏｃｅｓｓｉｎｇ）法を用いて、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を推定する。
第１の風速分布推定部１３は、Ｍが１であれば、観測領域を含んでいる２次元平面として、レンジビンＲｂ_１に対応する２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）を推定し、２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）を風速予測部１５に出力する。
第１の風速分布推定部１３は、Ｍが２以上であれば、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を第２の風速分布推定部１４に出力する。

第２の風速分布推定部１４は、例えば、図４に示す第２の風速分布推定回路２４によって実現される。
第２の風速分布推定部１４は、第１の風速分布推定部１３により推定された複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）から、複数の２次元平面の間の２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する。例えば、２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）は、レンジビンＲｂ_１に対応する２次元平面とレンジビンＲｂ_２に対応する２次元平面との間の２次元平面の風速分布である。２次元平面の風速分布ｕ_Ｒｂ２’（ｘ，ｙ）は、レンジビンＲｂ_２に対応する２次元平面とレンジビンＲｂ_３に対応する２次元平面との間の２次元平面の風速分布である。
第２の風速分布推定部１４は、複数の２次元平面の間の２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）を風速予測部１５に出力する。

風速予測部１５は、例えば、図４に示す風速予測回路２５によって実現される。
風速予測部１５は、Ｍが１であれば、第１の風速分布推定部１３により推定された２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）から、２次元ナビエストークス方程式（以下「２次元Ｎ－Ｓ方程式」という）を用いて、観測領域の風速ｕ（ｔ）を予測する。観測領域は、例えば、飛行体が飛行する予定の飛行領域を含んでいる領域である。
風速予測部１５は、Ｍが２以上であれば、第１の風速分布推定部１３により推定された複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）及び第２の風速分布推定部１４により推定された２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）の中から、飛行領域を観測領域として含んでいる２次元平面の風速分布を選択する。
風速予測部１５は、選択した２次元平面の風速分布から、２次元Ｎ－Ｓ方程式を用いて、飛行領域の風速ｕ（ｔ）を予測する。

図３では、風速予測装置３の構成要素である散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４及び風速予測部１５のそれぞれが、図４に示すような専用のハードウェアによって実現されるものを想定している。即ち、風速予測装置３が、散乱信号取得回路２１、ドップラー周波数算出回路２２、第１の風速分布推定回路２３、第２の風速分布推定回路２４及び風速予測回路２５によって実現されるものを想定している。

散乱信号取得回路２１、ドップラー周波数算出回路２２、第１の風速分布推定回路２３、第２の風速分布推定回路２４及び風速予測回路２５のそれぞれは、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ Ｓｐｅｃｉｆｉｃ Ｉｎｔｅｇｒａｔｅｄ Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ Ｇａｔｅ Ａｒｒａｙ）、又は、これらを組み合わせたものが該当する。

風速予測装置３の構成要素は、専用のハードウェアによって実現されるものに限るものではなく、風速予測装置３が、ソフトウェア、ファームウェア、又は、ソフトウェアとファームウェアとの組み合わせによって実現されるものであってもよい。
ソフトウェア又はファームウェアは、プログラムとして、コンピュータのメモリに格納される。コンピュータは、プログラムを実行するハードウェアを意味し、例えば、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、中央処理装置、処理装置、演算装置、マイクロプロセッサ、マイクロコンピュータ、プロセッサ、あるいは、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）が該当する。

図５は、風速予測装置３が、ソフトウェア又はファームウェア等によって実現される場合のコンピュータのハードウェア構成図である。
風速予測装置３が、ソフトウェア又はファームウェア等によって実現される場合、散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４及び風速予測部１５におけるそれぞれの処理手順をコンピュータに実行させるためのプログラムがメモリ３１に格納される。そして、コンピュータのプロセッサ３２がメモリ３１に格納されているプログラムを実行する。

また、図４では、風速予測装置３の構成要素のそれぞれが専用のハードウェアによって実現される例を示し、図５では、風速予測装置３がソフトウェア又はファームウェア等によって実現される例を示している。しかし、これは一例に過ぎず、風速予測装置３における一部の構成要素が専用のハードウェアによって実現され、残りの構成要素がソフトウェア又はファームウェア等によって実現されるものであってもよい。

次に、図１に示すレーダ装置１の動作について説明する。
図６は、風速予測装置３の処理手順である風速予測方法を示すフローチャートである。
ビーム送受信部２は、図２に示すように、仰角θ_ｎ（ｎ＝１，・・・，Ｎ）が互いに異なる複数のビームのそれぞれを空間に放射する。
即ち、ビーム送受信部２は、図２Ａに示すように、飛行領域が鉛直方向の真上に存在している場合、それぞれの仰角θ_ｎでビームスキャンを実施することで、それぞれの仰角θ_ｎにおいて、互いに異なる方向にビームを放射する。
ビーム送受信部２は、図２Ｂに示すように、飛行領域の位置が鉛直方向の真上からずれている場合、仰角θ_ｎを切り替えながら、それぞれの仰角θ_ｎでビームを放射する。ビーム送受信部２から放射されるビームは、例えば、パルス光である。
ビーム送受信部２から放射されたそれぞれのビームは、空間内に浮遊しているエアロゾル等の微粒子によって後方散乱される。例えば、１つのビームの視線方向に、複数の微粒子が存在していれば、１つのビームは、それぞれの微粒子によって後方散乱される。後方散乱された後のビームは、散乱信号として、ビーム送受信部２に戻ってくる。ビーム送受信部２から複数の微粒子までの距離は互いに異なるため、それぞれの散乱信号が、ビーム送受信部２に戻ってくる時間は、互いに異なる。
ビーム送受信部２は、それぞれの散乱信号をアナログ信号からデジタル信号ｓ_θｎ，ｊ（ｔ）（ｎ＝１，・・・，Ｎ：ｊ＝１，・・・，Ｊ）に変換し、デジタル信号ｓ_θｎ，ｊを風速予測装置３に出力する。
また、ビーム送受信部２は、空間に放射したそれぞれのビームの仰角θ_ｎを示す角度情報を風速予測装置３に出力する。

散乱信号取得部１１は、ビーム送受信部２から、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）を取得する（図６のステップＳＴ１）。
また、散乱信号取得部１１は、ビーム送受信部２から、それぞれのビームの仰角θ_ｎを示す角度情報を取得する。
散乱信号取得部１１は、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）及び角度情報のそれぞれをドップラー周波数算出部１２に出力する。

ドップラー周波数算出部１２は、散乱信号取得部１１から、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）を取得し、散乱信号取得部１１から、それぞれの角度情報を取得する。
ドップラー周波数算出部１２は、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）をヒット方向にＦＦＴする前に、角度情報が示す仰角θ_ｎに従って、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）のＦＦＴ長Ｌ_θｎを設定する（図６のステップＳＴ２）。ＦＦＴ長Ｌ_θｎは、散乱信号のレンジビン幅に相当する。
即ち、ドップラー周波数算出部１２は、ＦＦＴ長Ｌ_θｎに相当するレンジビン幅が１／ｓｉｎ（θ_ｎ）に比例するように、デジタル信号ｓ_θｎ，ｊ（ｔ）のＦＦＴ長Ｌ_θｎを設定する。
それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）のＦＦＴ長Ｌ_θｎが同じである場合、仰角θ_ｎが異なれば、図７Ａに示すように、同一のレンジビンＲｂ_ｍの複数の観測点の高度が互いに異なる。
一方、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）のＦＦＴ長Ｌ_θｎが仰角θ_ｎに従って設定されている場合、仰角θ_ｎが異なっていても、図７Ｂに示すように、同一のレンジビンＲｂ_ｍの複数の観測点の高度が同じ高度になる。
図１に示すレーダ装置１では、同一のレンジビンＲｂ_ｍの複数の観測点が同じ高度となるように、仰角θ_ｎが広いほど、ＦＦＴ長Ｌ_θｎが短く設定され、散乱信号の仰角θ_ｎが狭いほど、ＦＦＴ長Ｌ_θｎが長く設定されている。

図７は、ドップラー周波数の観測点を示す説明図である。
図７Ａは、仰角θ_ｎが異なれば、同一のレンジビンＲｂ_ｍの複数の観測点の高度が互いに異なることを示す説明図である。
図７Ｂは、仰角θ_ｎが変化しても、同一のレンジビンＲｂ_ｍの複数の観測点が同じ高度であることを示す説明図である。
ＦＦＴ長Ｌ_θｎが短く設定されれば、レンジビンＲｂ_ｍの幅が狭くなり、距離分解能が高くなる。一方、ＦＦＴ長Ｌ_θｎが長く設定されれば、レンジビンＲｂ_ｍの幅が広くなり、距離分解能が低くなる。
ドップラー周波数算出部１２は、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）を、ＦＦＴ長がＬ_θｎのＦＦＴをヒット方向に行うことで、それぞれのデジタル信号ｓ_θｎ，ｊ（ｔ）を周波数領域の信号ｒ_θｎ，ｊ（ｆ）に変換する（図６のステップＳＴ３）。

ドップラー周波数算出部１２は、周波数領域の信号ｒ_θｎ，ｊ（ｆ）から、複数のレンジビンＲｂ_１～Ｒｂ_Ｍにおけるそれぞれのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を算出する（図６のステップＳＴ４）。
ドップラー周波数算出部１２は、それぞれのレンジビンＲｂ_ｍのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を第１の風速分布推定部１３に出力する。
レンジビンＲｂ_ｍのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を算出する処理自体は、公知の技術であるため詳細な説明を省略する。

なお、デジタル信号ｓ_θｎ，ｊ（ｔ）のＳＮＲ（Ｓｉｇｎａｌ ｔｏ Ｎｏｉｓｅ Ｒａｔｉｏ）が低いために、ドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}の算出精度が低くなることがある。このような場合、ドップラー周波数算出部１２が、例えば、以下に示す処理を実施することで、ドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}の算出精度を高めるようにしてもよい。
ドップラー周波数算出部１２は、例えば、仰角θ_ｎが同じで、ヒット番号が互いに異なる複数のパルス光に係るデジタル信号ｓ_θｎ，ｊ（ｔ）をそれぞれ取得して、複数のビームに係るデジタル信号ｓ_θｎ，ｊ（ｔ）を周波数領域の信号ｒ_θｎ，ｊ（ｆ）にそれぞれ変換する。そして、ドップラー周波数算出部１２は、複数の周波数領域の信号ｒ_θｎ，ｊ（ｆ）をコヒーレント積分することで、ドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}の算出精度を高める。ドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}の算出精度を高める処理自体は、公知の技術であるため詳細な説明を省略する。

第１の風速分布推定部１３は、ドップラー周波数算出部１２から、それぞれのレンジビンＲｂ_ｍのドップラー周波数ｄｐｆ_{θｎ，ｊ，Ｒｂｍ}を取得する。
第１の風速分布推定部１３は、それぞれのレンジビンＲｂ_ｍにおける複数のドップラー周波数ｄｐｆ_{θ１，ｊ，Ｒｂｍ}～ｄｐｆ_{θＮ，ｊ，Ｒｂｍ}から、ＶＶＰ法を用いて、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を推定する（図６のステップＳＴ５）。

以下、第１の風速分布推定部１３による２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）の推定処理を具体的に説明する。
図８は、風速の観測点の位置（ｘ，ｙ，ｚ）を示す説明図である。
風速の観測点の位置（ｘ，ｙ，ｚ）は、以下の式（１）に示すように、レーダ装置１が備えるビーム送受信部２の位置Ｏを中心とする極座標で表されるものとする。

式（１）において、θは、ビームの仰角である。φは、観測点が存在しているｘ－ｙ平面の中心点Ｏ’と観測点の位置（ｘ，ｙ，ｚ）とを結ぶ線分と、ｙ軸とのなす角である。ｒは、ビーム送受信部２の位置Ｏと観測点の位置（ｘ，ｙ，ｚ）との間の距離である。

例えば、観測点の鉛直下方における地表面の地形が特殊な地形ではなく、一般的な地形であれば、風速予測装置３は、観測点が存在している２次元平面に対して、鉛直方向の風速を観測できなくても、２次元平面と平行な方向の風速を観測できれば、飛行体に影響する風を観測する上で、実用上問題がない。
一般的な地形とは、平坦な地形、傾斜している地形、又は、小さな凹凸がある地形等が該当する。特殊な地形とは、谷底がある地形、又は、標高差が大きな崖がある地形等が該当する。
このため、図３に示す風速予測装置３では、第１の風速分布推定部１３が、鉛直方向の風速を無視し、観測点における真の風速ベクトルｕが、以下の式（２）のように表されるものとする。

式（２）において、ｕ_ｘは、風速ベクトルｕのｘ軸と平行な方向の成分を示し、ｕ_ｙは、風速ベクトルｕのｙ軸と平行な方向の成分を示している。

第１の風速分布推定部１３は、ｘ－ｙ平面の中心点Ｏ’における風速ｕ_０が、以下の式（３）のように表されるものとして、中心点Ｏ’を中心とする風速ベクトルｕを、ｘ，ｙに関してテーラー展開する。例えば、３次の項以降を無視すれば、以下の近似式（４）が成立する。

レーダ装置１によって観測されるのは、観測点における風速の視線方向速度であるため、観測風速ｕ_ｒ（θ，φ）は、以下の式（５）のように表される。

第１の風速分布推定部１３は、Ｋ本のビームによる観測風速ｕ_ｒ（θ_１，φ_１）～ｕ_ｒ（θ_Ｋ，φ_Ｋ）を取得し、取得した観測風速ｕ_ｒ（θ_１，φ_１）～ｕ_ｒ（θ_Ｋ，φ_Ｋ）を以下の式（８）に代入することによって、式（７）に示すｐを推定する。Ｋは、２以上の整数である。Ｋ本のビームは、仰角θ及びなす角φの組が互いに異なるＫ個のパルス光である。

式（７）に示すｐの要素のうち、上から４番目、７番目及び８番目の要素は、２つのパラメータの和として表現されているため、不確定要素である。
第１の風速分布推定部１３は、式（８）によって推定したｐの要素のパラメータを均等に割り振ることで、不確定要素を確定する。
ｐの要素のうち、上から４番目がｐ_４であるとすれば、第１の風速分布推定部１３は、例えば、以下の式（９）のように、ｐ_４を決定する。
ｐの要素のうち、上から７番目がｐ_７であるとすれば、第１の風速分布推定部１３は、例えば、以下の式（１０）のように、ｐ_７を決定する。
ｐの要素のうち、上から８番目がｐ_８であるとすれば、第１の風速分布推定部１３は、例えば、以下の式（１１）のように、ｐ_８を決定する。

以上により、近似式（４）に示す全ての係数ｕ_０、∂ｕ／∂ｘ、∂ｕ／∂ｙ、∂^２ｕ／∂ｘ^２、∂^２ｕ／∂ｘ∂ｙ、∂^２ｕ／∂ｙ^２が決定されるため、第１の風速分布推定部１３は、近似式（４）より、それぞれの観測点の風速ベクトルｕ（ｘ，ｙ）を得ることができる。
第１の風速分布推定部１３は、それぞれの観測点の風速ベクトルｕ（ｘ，ｙ）を得ることができれば、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を得ることができる。
ここでは、第１の風速分布推定部１３が、風速ベクトルｕ（ｘ，ｙ）のテーラー展開を０次～２次の項で近似している。しかし、これは一例に過ぎず、第１の風速分布推定部１３が、風速ベクトルｕ（ｘ，ｙ）のテーラー展開を０次～３次以上の項で近似するようにしてもよい。

第１の風速分布推定部１３は、Ｍが１であれば、レンジビンＲｂ_１に対応する２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）を風速予測部１５に出力する。
第１の風速分布推定部１３は、Ｍが２以上であれば、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を第２の風速分布推定部１４に出力する。

第２の風速分布推定部１４は、第１の風速分布推定部１３から、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を取得する。
第２の風速分布推定部１４は、複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）から、複数の２次元平面の間の２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する（図６のステップＳＴ６）。
第２の風速分布推定部１４は、２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）を風速予測部１５に出力する。

以下、第２の風速分布推定部１４による２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）の推定処理を具体的に説明する。
例えば、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面と、風速分布がｕ_{Ｒｂｍ＋１}（ｘ，ｙ）である２次元平面との間の２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する場合、第２の風速分布推定部１４は、風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する２次元平面から、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面までの鉛直方向の距離Ｌ_ｍ－ｍ’を算出する。
また、第２の風速分布推定部１４は、風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する２次元平面から、風速分布がｕ_{Ｒｂｍ＋１}（ｘ，ｙ）である２次元平面までの鉛直方向の距離Ｌ_{（ｍ＋１）－ｍ’}を算出する。
第２の風速分布推定部１４は、以下の式（１２）から式（１４）に示すように、距離Ｌ_ｍ－ｍ’及び距離Ｌ_{（ｍ＋１）－ｍ’}に基づいて、２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）に対する重み係数ｗ_ｍと、２次元平面の風速分布ｕ_{Ｒｂｍ＋１}（ｘ，ｙ）に対する重み係数ｗ_ｍ＋１とを算出する。

第２の風速分布推定部１４は、以下の式（１５）に示すように、重み係数ｗ_ｍ及び重み係数ｗ_ｍ＋１に基づいて、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面と、風速分布が風速分布ｕ_{Ｒｂｍ＋１}（ｘ，ｙ）である２次元平面との間の２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する。

ここでは、第２の風速分布推定部１４が、風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定する２次元平面からの距離に基づく重み付き平均として、２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定している。しかし、これは一例に過ぎず、例えば、第２の風速分布推定部１４が、距離Ｌ_ｍ－ｍ’と距離Ｌ_{（ｍ＋１）－ｍ’}とを比較する。そして、第２の風速分布推定部１４は、距離Ｌ_ｍ－ｍ’が距離Ｌ_{（ｍ＋１）－ｍ’}よりも短ければ、２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）として、２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を用い、距離Ｌ_ｍ－ｍ’が距離Ｌ_{（ｍ＋１）－ｍ’}以上であれば、２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）として、２次元平面の風速分布ｕ_{Ｒｂｍ＋１}（ｘ，ｙ）を用いるようにしてもよい。
また、ここでは、第２の風速分布推定部１４が、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面と、風速分布がｕ_{Ｒｂｍ＋１}（ｘ，ｙ）である２次元平面との間の２次元平面の風速分布として、１つの２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定している。しかし、これは一例に過ぎず、第２の風速分布推定部１４は、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面と、風速分布がｕ_{Ｒｂｍ＋１}（ｘ，ｙ）である２次元平面との間の２次元平面の風速分布として、鉛直方向の位置が互いに異なる、２つ以上の２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を推定するようにしてもよい。

風速予測部１５は、Ｍが１であれば、第１の風速分布推定部１３から、２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）を取得する。
風速予測部１５は、Ｍが２以上であれば、第１の風速分布推定部１３から、複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）を取得し、第２の風速分布推定部１４から、２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）を取得する。

風速予測部１５は、Ｍが１であれば、２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）から、２次元Ｎ－Ｓ方程式を用いて、当該２次元平面内の観測領域の風速ｕ（ｔ）を予測する（図６のステップＳＴ７）。
風速予測部１５は、Ｍが２以上であれば、複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）及び２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）の中から、いずれかの２次元平面の風速分布を選択する。
風速ｕ（ｔ）の予測対象の観測点が、レンジビンＲｂ_ｍに対応する２次元平面内の観測領域に存在していれば、風速予測部１５は、風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を選択する。風速ｕ（ｔ）の予測対象の観測点が、レンジビンＲｂ_ｍ’に対応する２次元平面内の観測領域に存在していれば、風速予測部１５は、風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）を選択する。
風速予測部１５は、選択した２次元平面の風速分布から、２次元Ｎ－Ｓ方程式を用いて、選択した２次元平面内の観測領域の風速ｕ（ｔ）を予測する（図６のステップＳＴ７）。

以下、風速予測部１５による風速ｕ（ｔ）の予測処理を具体的に説明する。
風は、一般的に、非圧縮流体とみなされる。非圧縮流体に関する２次元Ｎ－Ｓ方程式は、以下に示す式（１６）及び式（１７）によって表現される。

式（１６）において、ｕは風速ベクトル、νは粘性係数、ρは密度、ｐは圧力、ｆは外力ベクトル、ｔは時間変数である。

粘性係数νは、非常に小さいので、風の粘性を無視することができる。また、外力ベクトルＦの項も無視することができる。よって、風の粘性及び外力のそれぞれが存在しないという条件の下では、式（１６）は、以下の式（１８）のように表され、式（１７）は、以下の式（１９）のように表される。

式（１８）の右辺第１項は、移流項と呼ばれ、右辺第２項は、圧力項と呼ばれる。
風速ベクトルｕが２次元ベクトルｕ＝［ｕ_ｘ，ｕ_ｙ］^Ｔであるとすれば、式（１８）は、以下の式（２０）及び式（２１）のように表される。

式（２２）を用いると、式（２０）は、以下の式（２３）のように表され、式（２１）は、以下の式（２４）のように表される。

差分法によって、風速ベクトルｕを、時間に関して離散化すると、式（２３）は、以下の式（２５）のように表され、式（２４）は、以下の式（２６）のように表される。

式（２５）及び式（２６）において、ｇ，ｇ＋１のそれぞれは、時間ステップを示す記号である。

中間変数ｕ_ｘ ^＊を用いれば、式（２５）は、以下の式（２７）及び式（２８）のように、２つの式に分けることができる。
また、中間変数ｕ_ｙ ^＊を用いれば、式（２６）は、以下の式（２９）及び式（３０）のように、２つの式に分けることができる。

式（２７）は、式（３１）のように変形され、式（２９）は、式（３２）のように変形される。

式（２８）は、ｘ方向に微分すれば、以下の式（３３）のように表され、式（３０）は、ｙ方向に微分すれば、以下の式（３４）のように表される。

式（３３）及び式（３４）の左辺同士が加算され、式（３３）及び式（３４）の右辺同士が加算されると、以下の式（３５）のように表される。

式（２２）より、推定される次時刻ステップｇ＋１の風速ベクトルｕに関しては、∇・ｕ^ｇ＋１＝０である必要がある。このため、式（３５）より、圧力に関するポアソン方程式として、以下の式（３６）が得られる。

風速予測部１５は、現在の時刻ステップｇの風速ベクトルｕ^ｇ＝［ｕ_ｘ ^ｇ，ｕ_ｙ ^ｇ］^Ｔを空間方向に離散化する。
風速予測部１５は、離散化したｕ_ｘ ^ｇを式（３１）に代入することで、中間変数ｕ_ｘ ^＊を算出し、離散化したｕ_ｙ ^ｇを式（３２）に代入することで、中間変数ｕ_ｙ ^＊を算出する。
風速予測部１５は、中間変数ｕ_ｘ ^＊及び中間変数ｕ_ｙ ^＊を用いて、式（３６）から圧力ｐを算出する。
風速予測部１５は、中間変数ｕ_ｘ ^＊及び圧力ｐを用いて、式（２８）から次時刻ステップｇ＋１の風速ベクトルｕ_ｘ ^ｇ＋１を算出する。
風速予測部１５は、中間変数ｕ_ｙ ^＊及び圧力ｐを用いて、式（３０）から次時刻ステップｇ＋１の風速ベクトルｕ_ｙ ^ｇ＋１を算出する。
風速予測部１５により予測される風速ｕ（ｔ）は、［ｕ_ｘ ^ｇ＋１，ｕ_ｙ ^ｇ＋１］^Ｔであり、風速ｕ（ｔ）は、例えば、図示せぬ表示器に表示される。

従来の風速予測装置は、観測領域の風速を予測するために、気象モデルを用いるシミュレーションを実施しており、気象モデルを用いるシミュレーションは、３次元の演算処理である。例えば、３次元の観測領域内の観測点の数が、Ｃｘ×Ｃｙ×Ｃｚであり、Ｃ＝Ｃｘ＝Ｃｙ＝Ｃｚであるとすれば、演算量オーダは、概ねＣ^３に比例するオーダとなる。
これに対して、図３に示す風速予測装置３の演算処理は、２次元の演算処理である。例えば、２次元平面の数がＭであり、２次元平面に存在している観測領域内の観測点の数が、Ｃｘ×Ｃｙであり、Ｃ＝Ｃｘ＝Ｃｙであるとすれば、演算量オーダは、概ねＣ^２×Ｍに比例するオーダとなる。Ｍ＝Ｃであれば、図３に示す風速予測装置３の演算量オーダは、従来の風速予測装置の演算量オーダと、概ね同じになる。
しかしながら、３次元の演算処理は、鉛直風を無視しない３次元Ｎ－Ｓ方程式を用いることが想定される。当該３次元Ｎ－Ｓ方程式では、鉛直方向の観測点の間隔を大幅に広げることが困難であり、３次元Ｎ－Ｓ方程式を解くことが可能な、観測点の間隔の最大値は、鉛直方向のクーラン条件として決まっている。
一方、２次元の演算処理である、鉛直風を無視する２次元Ｎ－Ｓ方程式では、鉛直方向のクーラン条件が存在していないので、Ｃｚ≫Ｍとすることが可能である。例えば、（Ｍ／Ｃｚ）＝１／Ｑであるとすれば、図３に示す風速予測装置３の演算量オーダは、従来の風速予測装置の演算量オーダの概ね１／Ｑとなる。Ｑは、１０～１０００程度の値である。
上述したように、例えば、観測点の鉛直下方における地表面の地形が特殊な地形ではなく、一般的な地形であれば、風速予測装置３は、観測点が存在している２次元平面に対して、鉛直方向の風速を観測できなくても、２次元平面と平行な方向の風速を観測できれば、飛行体に影響する風を観測する上で、実用上問題がない。
仮に、観測点の鉛直下方における地表面の地形が特殊な地形であっても、２次元平面と平行な方向の風速は、飛行体に影響する風を観測する上で、有用な情報として利用可能である。

以上の実施の形態１では、視線方向と水平方向とのなす角である仰角が互いに異なる複数のビームのそれぞれが空間に放射され、空間で散乱された後のそれぞれのビームである散乱信号を取得する散乱信号取得部１１と、空間に放射されたそれぞれのビームの仰角に従って、散乱信号取得部１１により取得されたそれぞれの散乱信号の距離分解能としてレンジビン幅を設定し、それぞれの散乱信号から、観測領域を含んでいる２次元平面に対応するレンジビンのドップラー周波数を算出するドップラー周波数算出部１２と、ドップラー周波数算出部１２により算出された複数のドップラー周波数から、ＶＶＰ法を用いて、２次元平面の風速分布を推定する第１の風速分布推定部１３と、第１の風速分布推定部１３により推定された２次元平面の風速分布から、２次元ナビエストークス方程式を用いて、観測領域の風速を予測する風速予測部１５とを備えるように、風速予測装置３を構成した。したがって、風速予測装置３は、３次元の演算処理を実施することなく、風速を予測することができる。

図３に示す風速予測装置３では、第１の風速分布推定部１３が、それぞれのレンジビンＲｂ_ｍに対応する２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を推定している。しかしながら、第１の風速分布推定部１３により推定された風速分布ｕ_Ｒｂｍ（ｘ，ｙ）は、観測ノイズの影響によって、実際の風速分布から多少のずれを生じることがある。第１の風速分布推定部１３は、観測ノイズの影響によるずれを抑えるため、風速分布ｕ_Ｒｂｍ（ｘ，ｙ）を時間方向に平滑化するようにしてもよい。
風速分布ｕ_Ｒｂｍ（ｘ，ｙ）の時間方向の平滑化は、２次元Ｎ－Ｓ方程式のような流体移動モデルと、推定した風速分布ｕ_Ｒｂｍ（ｘ，ｙ）とを融合することで、実現することが可能である。流体移動モデルと風速分布ｕ_Ｒｂｍ（ｘ，ｙ）との融合は、データ同化と呼ばれる。
例えば、前時間ステップｇ－１の真の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）がｕ_{ｔ＝ｔｇ－１}で表され、現在の時間ステップｇの真の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）がｕ_ｔ＝ｔｇで表されるものとする。２次元Ｎ－Ｓ方程式は、ｕ_{ｔ＝ｔｇ－１}とｕ_ｔ＝ｔｇとの関係を非線形に記述することが可能な流体モデルであり、非線形関数がｆ（・）で表されるとすれば、以下の式（３７）が成立する。

データ同化処理は、以下の式（３８）に示すように、データ同化関数Ａｓ（・）を用いて実施する処理である。即ち、データ同化処理は、現在の時間ステップｇの風速分布ｕ_Ｒｂｍ（ｘ，ｙ）であるｕ（ｘ，ｙ；_ｔ＝ｔｇ）と、前時間ステップｇ－１の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）が現在の時間ステップｇまで時間発展されたｆ（ｕ（ｘ，ｙ；_{ｔ＝ｔｇ－１}）バー）とを足し合わせることで、現在の時間ステップｇの平滑化されたｕ（ｘ，ｙ；_ｔ＝ｔｇ）バーを算出するものである。明細書の文章中では、電子出願の関係上、文字の上に“－”の記号を付することができないので、例えば、ｕ（ｘ，ｙ；_ｔ＝ｔｇ）バーのように表記している。

ここでは、第１の風速分布推定部１３が、データ同化関数Ａｓ（・）を用いて、データ同化処理を実施している。しかし、これは一例に過ぎず、第１の風速分布推定部１３が、例えば、２次元Ｎ－Ｓ方程式を用いたパーティクルフィルタ、又は、２次元Ｎ－Ｓ方程式を用いたアンサンブルカルマンフィルタによって、データ同化処理を実施するようにしてもよい。また、第１の風速分布推定部１３が、２次元Ｎ－Ｓ方程式を線形関数に近似した拡張カルマンフィルタ、又は、２次元Ｎ－Ｓ方程式を用いた変分法によって、データ同化処理を実施するようにしてもよい。

実施の形態２．
実施の形態１に係る風速予測装置３では、風速予測部１５が、複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）及び２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）の中から、いずれかの２次元平面の風速分布を選択している。
実施の形態２では、複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）と２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）とから、飛行体が飛行する予定の飛行領域を観測領域として含んでいる２次元平面の風速分布を推定する風速予測部１６を備える風速予測装置３について説明する。

図９は、実施の形態２に係る風速予測装置３を示す構成図である。図９において、図３と同一符号は同一又は相当部分を示すので説明を省略する。
図１０は、実施の形態２に係る風速予測装置３のハードウェアを示すハードウェア構成図である。図１０において、図４と同一符号は同一又は相当部分を示すので説明を省略する。
図９に示す風速予測装置３は、散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４及び風速予測部１６を備えている。

風速予測部１６は、例えば、図１０に示す風速予測回路２６によって実現される。
風速予測部１６は、第１の風速分布推定部１３により推定された複数の２次元平面の風速分布ｕ_Ｒｂ１（ｘ，ｙ）～ｕ_ＲｂＭ（ｘ，ｙ）と、第２の風速分布推定部１４により推定された２次元平面の風速分布ｕ_Ｒｂ１’（ｘ，ｙ）～ｕ_{ＲｂＭ－１’}（ｘ，ｙ）とから、飛行領域を観測領域として含んでいる２次元平面の風速分布を推定する。
風速予測部１６は、推定した２次元平面の風速分布から、２次元Ｎ－Ｓ方程式を用いて、飛行領域の風速ｕ（ｔ）を予測する。

図９では、風速予測装置３の構成要素である散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４及び風速予測部１６のそれぞれが、図１０に示すような専用のハードウェアによって実現されるものを想定している。即ち、風速予測装置３が、散乱信号取得回路２１、ドップラー周波数算出回路２２、第１の風速分布推定回路２３、第２の風速分布推定回路２４及び風速予測回路２６によって実現されるものを想定している。

散乱信号取得回路２１、ドップラー周波数算出回路２２、第１の風速分布推定回路２３、第２の風速分布推定回路２４及び風速予測回路２６のそれぞれは、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ、ＦＰＧＡ、又は、これらを組み合わせたものが該当する。
風速予測装置３の構成要素は、専用のハードウェアによって実現されるものに限るものではなく、風速予測装置３が、ソフトウェア、ファームウェア、又は、ソフトウェアとファームウェアとの組み合わせによって実現されるものであってもよい。

風速予測装置３が、ソフトウェア又はファームウェア等によって実現される場合、散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４及び風速予測部１６におけるそれぞれの処理手順をコンピュータに実行させるためのプログラムが図５に示すメモリ３１に格納される。そして、図５に示すプロセッサ３２がメモリ３１に格納されているプログラムを実行する。
また、図１０では、風速予測装置３の構成要素のそれぞれが専用のハードウェアによって実現される例を示し、図５では、風速予測装置３がソフトウェア又はファームウェア等によって実現される例を示している。しかし、これは一例に過ぎず、風速予測装置３における一部の構成要素が専用のハードウェアによって実現され、残りの構成要素がソフトウェア又はファームウェア等によって実現されるものであってもよい。

次に、図９に示す風速予測装置３の動作について説明する。風速予測部１６以外は、図３に示す風速予測装置３と同様であるため、ここでは、風速予測部１６の動作のみを説明する。
ここでは、説明の便宜上、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面と、風速分布がｕ_Ｒｂｍ’（ｘ，ｙ）である２次元平面との間の２次元平面に、飛行領域が存在しているものとする。

風速予測部１６は、飛行領域が存在している２次元平面から、風速分布がｕ_Ｒｂｍ（ｘ，ｙ）である２次元平面までの鉛直方向の距離Ｌ_ｐ－ｍを算出する。
また、風速予測部１６は、飛行領域が存在している２次元平面から、風速分布がｕ_Ｒｂｍ’（ｘ，ｙ）である２次元平面までの鉛直方向の距離Ｌ_ｐ－ｍ’を算出する。
風速予測部１６は、以下の式（３９）から式（４１）に示すように、距離Ｌ_ｐ－ｍ及び距離Ｌ_ｐ－ｍ’に基づいて、２次元平面の風速分布ｕ_Ｒｂｍ（ｘ，ｙ）に対する重み係数ｗ_ｐ－ｍと、２次元平面の風速分布ｕ_Ｒｂｍ’（ｘ，ｙ）に対する重み係数ｗ_ｐ－ｍ’とを算出する。

風速予測部１６は、以下の式（４２）に示すように、重み係数ｗ_ｐ－ｍ及び重み係数ｗ_ｐ－ｍ’に基づいて、飛行領域が存在している２次元平面の風速分布ｕ_Ｒｂｐ（ｘ，ｙ）を推定する。

風速予測部１６は、推定した２次元平面の風速分布ｕ_Ｒｂｐ（ｘ，ｙ）から、２次元Ｎ－Ｓ方程式を用いて、飛行領域の風速ｕ（ｔ）を予測する。風速予測部１６による風速ｕ（ｔ）の予測処理自体は、図３に示す風速予測部１５による風速ｕ（ｔ）の予測処理と同様である。

実施の形態３．
実施の形態３では、風速予測部１５により予測された風速ｕ（ｔ）が閾値Ｔｈ以上となる時刻ｔを算出し、時刻ｔを通知する時刻通知部１７を備える風速予測装置３について説明する。

図１１は、実施の形態３に係る風速予測装置３を示す構成図である。図１１において、図３及び図９と同一符号は同一又は相当部分を示すので説明を省略する。
図１２は、実施の形態３に係る風速予測装置３のハードウェアを示すハードウェア構成図である。図１２において、図４及び図１０と同一符号は同一又は相当部分を示すので説明を省略する。
図１１に示す風速予測装置３は、散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４、風速予測部１５及び時刻通知部１７を備えている。

時刻通知部１７は、例えば、図１２に示す時刻通知回路２７によって実現される。
時刻通知部１７は、風速予測部１５により予測された風速ｕ（ｔ）が閾値Ｔｈ以上となる時刻ｔを特定する。
時刻通知部１７は、特定した時刻ｔを、例えば、飛行体、又は、飛行体の飛行状態を監視している管制塔に通知する。

図１１は、時刻通知部１７が図３に示す風速予測装置３に適用されている例を示している。しかし、これは一例に過ぎず、時刻通知部１７が図９に示す風速予測装置３に適用されているものであってもよい。

図１１では、風速予測装置３の構成要素である散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４、風速予測部１５及び時刻通知部１７のそれぞれが、図１２に示すような専用のハードウェアによって実現されるものを想定している。即ち、風速予測装置３が、散乱信号取得回路２１、ドップラー周波数算出回路２２、第１の風速分布推定回路２３、第２の風速分布推定回路２４、風速予測回路２５及び時刻通知回路２７によって実現されるものを想定している。

散乱信号取得回路２１、ドップラー周波数算出回路２２、第１の風速分布推定回路２３、第２の風速分布推定回路２４、風速予測回路２５及び時刻通知回路２７のそれぞれは、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ、ＦＰＧＡ、又は、これらを組み合わせたものが該当する。
風速予測装置３の構成要素は、専用のハードウェアによって実現されるものに限るものではなく、風速予測装置３が、ソフトウェア、ファームウェア、又は、ソフトウェアとファームウェアとの組み合わせによって実現されるものであってもよい。

風速予測装置３が、ソフトウェア又はファームウェア等によって実現される場合、散乱信号取得部１１、ドップラー周波数算出部１２、第１の風速分布推定部１３、第２の風速分布推定部１４、風速予測部１５及び時刻通知部１７におけるそれぞれの処理手順をコンピュータに実行させるためのプログラムが図５に示すメモリ３１に格納される。そして、図５に示すプロセッサ３２がメモリ３１に格納されているプログラムを実行する。
また、図１２では、風速予測装置３の構成要素のそれぞれが専用のハードウェアによって実現される例を示し、図５では、風速予測装置３がソフトウェア又はファームウェア等によって実現される例を示している。しかし、これは一例に過ぎず、風速予測装置３における一部の構成要素が専用のハードウェアによって実現され、残りの構成要素がソフトウェア又はファームウェア等によって実現されるものであってもよい。

次に、図１１に示す風速予測装置３の動作について説明する。時刻通知部１７以外は、図３に示す風速予測装置３と同様であるため、ここでは、時刻通知部１７の動作のみを説明する。
時刻通知部１７は、風速予測部１５が風速ｕ（ｔ）を予測する毎に、風速予測部１５から、予測された風速ｕ（ｔ）を取得する。
時刻通知部１７は、予測された風速ｕ（ｔ）を取得する毎に、風速ｕ（ｔ）と閾値Ｔｈとを比較する。閾値Ｔｈは、例えば、飛行体の飛行に影響を及ぼす可能性がある風速の下限値である。閾値Ｔｈは、時刻通知部１７の内部メモリに格納されていてもよいし、外部から与えられるものであってもよい。

時刻通知部１７は、風速ｕ（ｔ）と閾値Ｔｈとの比較結果に基づいて、風速ｕ（ｔ）が閾値Ｔｈ以上となる時刻ｔを特定する。
時刻通知部１７は、特定した時刻ｔを、例えば、飛行体、又は、飛行体の飛行状態を監視している管制塔に通知する。

以上の実施の形態３では、風速予測部１５により予測された風速が閾値以上となる時刻を特定し、特定した時刻ｔを通知する時刻通知部１７を備えるように、図１１に示す風速予測装置３を構成した。したがって、図１１に示す風速予測装置３は、図３に示す風速予測装置３と同様に、３次元の演算処理を実施することなく、風速を予測することができる。また、飛行体の飛行に影響を及ぼす可能性がある風が吹く時刻を知らせることができる。

なお、本開示は、各実施の形態の自由な組み合わせ、あるいは各実施の形態の任意の構成要素の変形、もしくは各実施の形態において任意の構成要素の省略が可能である。

本開示は、風速予測装置、風速予測方法及びレーダ装置に適している。

１ レーダ装置、２ ビーム送受信部、３ 風速予測装置、１１ 散乱信号取得部、１２ ドップラー周波数算出部、１３ 第１の風速分布推定部、１４ 第２の風速分布推定部、１５，１６ 風速予測部、１７ 時刻通知部、２１ 散乱信号取得回路、２２ ドップラー周波数算出回路、２３ 第１の風速分布推定回路、２４ 第２の風速分布推定回路、２５，２６ 風速予測回路、２７ 時刻通知回路、３１ メモリ、３２ プロセッサ。

Claims (7)

1. 視線方向と水平方向とのなす角である仰角が互いに異なる複数のビームのそれぞれが空間に放射され、前記空間で散乱された後のそれぞれのビームである散乱信号を取得する散乱信号取得部と、
   前記空間に放射されたそれぞれのビームの仰角に従って、前記散乱信号取得部により取得されたそれぞれの散乱信号の距離分解能としてレンジビン幅を可変に設定することで、仰角の異なるビームであっても同一高度の観測点を観測できるようにし、それぞれの散乱信号から、観測領域を含んでいる２次元平面に対応するレンジビンのドップラー周波数を算出するドップラー周波数算出部と、
   前記ドップラー周波数算出部により算出された複数のドップラー周波数から、ＶＶＰ（Ｖｏｌｕｍｅ Ｖｅｌｏｃｉｔｙ Ｐｒｏｃｅｓｓｉｎｇ）法を用いて、前記２次元平面の風速分布を推定する第１の風速分布推定部と、
   前記第１の風速分布推定部により推定された２次元平面の風速分布から、２次元ナビエストークス方程式を用いて、前記観測領域の風速を予測する風速予測部と
   を備えた風速予測装置。
2. 前記ドップラー周波数算出部は、それぞれの散乱信号から、複数のレンジビンにおけるそれぞれのドップラー周波数を算出することで、観測する３次元領域を複数の２次元平面の積み上げとして表現し、
   前記第１の風速分布推定部は、前記ドップラー周波数算出部により算出された高度の異なる各２次元平面に対応するレンジビンにおける複数のドップラー周波数から、ＶＶＰ法を用いて、それぞれのレンジビンに対応する高度の異なる２次元平面の風速分布を推定することを特徴とする請求項１記載の風速予測装置。
3. 前記風速予測部は、前記第１の風速分布推定部により推定された複数の高度の異なる２次元平面の風速分布の中から、飛行体が飛行する予定の飛行高度に最も近い高度の２次元平面の風速分布を選択し、選択した風速分布から、２次元ナビエストークス方程式を用いて、前記飛行体が飛行する予定の飛行領域の風速を予測することを特徴とする請求項２記載の風速予測装置。
4. 前記第１の風速分布推定部により推定された高度の異なる複数の２次元平面の風速分布から、前記複数の２次元平面の間の２次元平面の風速分布を推定して、２次元平面の風速分布を高度方向に内挿補間することにより、指定された高度の２次元風速分布を推定する第２の風速分布推定部を備え、
   前記風速予測部は、飛行体が飛行する予定の飛行高度における２次元平面の風速分布を前記第２の風速分布推定部を用いて推定し、推定した風速分布から、２次元ナビエストークス方程式を用いて、前記飛行体が飛行する予定の飛行領域の風速を予測することを特徴とする請求項２記載の風速予測装置。
5. 前記風速予測部により予測された風の速度が閾値以上となる時刻を特定し、前記時刻を通知する時刻通知部を備えたことを特徴とする請求項１記載の風速予測装置。
6. 視線方向と水平方向とのなす角である仰角が互いに異なる複数のビームのそれぞれが空間に放射されたとき、
   散乱信号取得部が、前記空間で散乱された後のそれぞれのビームである散乱信号を取得し、
   ドップラー周波数算出部が、前記空間に放射されたそれぞれのビームの仰角に従って、前記散乱信号取得部により取得されたそれぞれの散乱信号の距離分解能としてレンジビン幅を可変に設定することで、仰角の異なるビームであっても同一高度の観測点を観測できるようにし、それぞれの散乱信号から、観測領域を含んでいる２次元平面に対応するレンジビンのドップラー周波数を算出し、
   第１の風速分布推定部が、前記ドップラー周波数算出部により算出された複数のドップラー周波数から、ＶＶＰ（Ｖｏｌｕｍｅ Ｖｅｌｏｃｉｔｙ Ｐｒｏｃｅｓｓｉｎｇ）法を用いて、前記２次元平面の風速分布を推定し、
   風速予測部が、前記第１の風速分布推定部により推定された２次元平面の風速分布から、２次元ナビエストークス方程式を用いて、前記観測領域の風速を予測する
   風速予測方法。
7. 視線方向と水平方向とのなす角である仰角が互いに異なる複数のビームのそれぞれを空間に放射し、前記空間で散乱された後のそれぞれのビームである散乱信号を受信するビーム送受信部と、
   前記ビーム送受信部により受信されたそれぞれの散乱信号を取得する散乱信号取得部と、
   前記ビーム送受信部から空間に放射されたそれぞれのビームの仰角に従って、前記散乱信号取得部により取得されたそれぞれの散乱信号の距離分解能としてレンジビン幅を可変に設定することで、仰角の異なるビームであっても同一高度の観測点を観測できるようにし、それぞれの散乱信号から、観測領域を含んでいる２次元平面に対応するレンジビンのドップラー周波数を算出するドップラー周波数算出部と、
   前記ドップラー周波数算出部により算出された複数のドップラー周波数から、ＶＶＰ（Ｖｏｌｕｍｅ Ｖｅｌｏｃｉｔｙ Ｐｒｏｃｅｓｓｉｎｇ）法を用いて、前記２次元平面の風速分布を推定する第１の風速分布推定部と、
   前記第１の風速分布推定部により推定された２次元平面の風速分布から、２次元ナビエストークス方程式を用いて、前記観測領域の風速を予測する風速予測部と
   を備えたレーダ装置。

Images