JP2010096593A - 風予測装置及びプログラム - Google Patents

Abstract

【課題】
短時間後の強風域を局所的に予測する。
【解決手段】
風向・風速算出装置（１４）は、複数のドップラーレーダ（１０−１〜１０−ｎ）により風観測されたドップラー速度データとレーダ反射因子データから風向・風速分布データを計算する。データ同化装置（３０）は、風向・風速分布データを雲解像数値気象モデルにデータ同化し、所定格子間隔の風速場データを生成する。移動ベクトル推定装置（３６）は、現在及び過去の一定時間内の複数の風速場データから風の時間変化を示す移動ベクトルを推定する。上空風速場予測装置（３８）は、装置（３６）による移動ベクトルに従い、現在の風速場データを時間方向に外挿して上空風速場を予測する。地上強風域予測装置（４０）は、装置（３８）の予測結果に従い地上強風風域を予測する。
【選択図】 図１

Description

本発明は、風予測装置及びプログラムに関し、より具体的には、地上の風向・風速、特に強風域を予測する装置及びコンピュータプログラムに関する。

風（風向と風速）を予測することは、例えば、局所的に発生する強風への対策として有用であり、また、降雨予測の精度向上にも寄与する。風を予測するには、現在の風向・風速を精密に観測する技術又は体制と、観測値から精度よく予測可能なシステムの両方が必要である。

気象の時間変化を包括的に予測する手法として、３次元空間の所定距離間隔で想定されるメッシュ又は格子点上の気象データに一定の物理方程式を適用し、その時間変換をコンピュータシミュレーションする方法が知られている。特許文献１には、気象データサーバＤＳから定期的に提供される広域気象予測データを初期値として気象予測モデルを作成し、このモデルに気象状況の観測結果をデータ同化し、予測モデルから風分布を含む極細密気象予測演算を行い、その演算結果を気象予測情報として提供するシステムが記載されている。

また、風向・風速を計測する方法として、ゾンデやウインドプロファイラを用いた高層気象観測、地上の風向・風速計による観測、及び、ドップラーレーダによる観測がある。
特開２００７−０１７３１６号公報

ゾンデやウインドプロファイラによる観測は、空間的にも時間的にも粗いために、強風の予測、特に、数分とか数十分後等の短時間後の予測（ナウキヤスト）には利用できない。

その点、ドップラーレーダによる観測は空間的時間的に十分な範囲のデータを得ることができる。しかし、１台のドップラーレーダから求まる情報は、動径線上の、近づく方向から遠ざかる方向での速度分布であり、そのままでは風向・風速の分布を求めることはできない。ＶＡＤやＶＶＰなどの従来手法は、風の場が一様あるいは線形であるという仮定のもとで１台のドップラーレーダの観測結果から風向・風速を推定するものであり、求まる風向・風速は、ある領域の平均的な値になる。すなわち、風が急変するような現象は捉えられない。

このような状況から、現在まで、５００ｍというような密な空間分解能で風のナウキヤストを行っている例はない。

本発明は、局所的な強風の発生を的確に予測可能な風予測装置及びプログラムを提示することを目的とする。

本発明に係る風予測装置は、複数のドップラーレーダにより風観測されたドップラー速度データを記憶する記憶手段と、当該ドップラー速度データを折り返し補正する折り返し補正手段と、当該折り返し補正手段により補正されたドップラー速度データから、地形データを参照して風向・風速分布データを計算する風向・風速分布計算手段と、当該風向・風速分布データを数値気象モデルにデータ同化し、所定格子間隔の風速場データを生成するデータ同化手段と、所定時間内の複数の風速場データから、風の移動ベクトルを推定する移動ベクトル推定手段と、当該移動ベクトルに従い、現在の風速場データを時間方向に外挿して上空風速場を予測する上空風速場予測手段と、当該上空風速場予測手段の予測結果に従い所定の地上風域を予測する地上風域予測手段とを具備することを特徴とする。

本発明に係る風予測プログラムは、複数のドップラーレーダにより風観測されたドップラー速度データをコンピュータの記憶手段に格納する機能と、当該コンピュータに、当該ドップラー速度データを折り返し補正させる折り返し補正機能と、当該コンピュータに、当該折り返し補正機能により補正されたドップラー速度データから、地形データを参照して風向・風速分布データを計算させる風向・風速分布計算機能と、当該コンピュータに、当該風向・風速分布データを数値気象モデルにデータ同化させ、所定格子間隔の風速場データを生成させるデータ同化機能と、当該コンピュータに、所定時間内の複数の風速場データから、風の移動ベクトルを推定させる移動ベクトル推定機能と、当該コンピュータに、当該移動ベクトルに従い、現在の風速場データを時間方向に外挿して上空風速場を予測させる上空風速場予測機能と、当該コンピュータに、当該上空風速場予測機能の予測結果に従い所定の地上風域を予測させる地上風域予測機能とを具備することを特徴とする。

本発明によれば、短時間後の地上風、例えば、強風域を局所的に精度よく予測することが可能になる。

以下、図面を参照して、本発明の実施例を詳細に説明する。

図１は、本発明の一実施例の概略構成ブロック図を示す。本実施例は、理解を容易にするために、１台のコンピュータ上で動作するとしているが、勿論、複数のコンピュータを協調動作させることで、大量のデータを扱うことが可能になり、高速に結果を得ることができることは明らかである。

ドップラーレーダは、電波を空中に放射し、その反射波のドップラー効果による周波数変化から、空中浮遊物質、例えば、降水粒子の動径方向での移動速度（ドップラー速度）を計測するものである。降水粒子は、水平方向にはその場の水平風速度で移動し、鉛直方向にはその場の鉛直風速度と大気に相対的な落下速度の和で移動すると考えられる。ドップラー速度は、測定位置の風速度と落下速度の和の動径方向成分である。一般に、降水粒子の落下速度とレーダ反射因子（受信電波の強度）との関係式が知られており、レーダ反射因子から推定された降水落下速度の動径方向成分を求め、これをドップラー速度から差し引くことにより、風速度の動径方向成分を知ることができる。

通常、ドップラーレーダの観測では、ある仰角に固定した状態で方位角方向に一回転させるＰＰＩ（Plan Position Indicator）観測を基本としつつ、一定時間内に複数の仰角のＰＰＩ観測を行う（ボリュームスキャン）ことで、三次元的データを取得する。ドップラーレーダで観測されたドップラー速度およびレーダ反射因子データは、測定位置の仰角、方位角及び動径距離の情報とともに保存される。

１台のドップラーレーダでは、電波放射軸線上の速度成分（動径方向成分）しか計測できないが、複数のドップラーレーダの放射ビームを所定範囲の角度で交差させ、それぞれの電波で計測された速度成分をベクトル合成することで、その交差点上の風速度（風向・風速）を計測できる。このベクトル合成では、交差角度が９０度からプラスマイナス６０度の範囲で有意な結果を得ることができる。風速度を求めたい領域がこの交差角の範囲に収まるようにレーダを設置する必要がある。例えば観測範囲を８０ｋｍ四方とする場合、１辺の距離が４０ｋｍ程度の正三角形の頂点にドップラーレーダを配置すればよい。

本実施例では、ｎ台のドップラーレーダ１０−１〜１０−ｎが計測対象空間の風向・風速を計測できるように、分散して配置されている。各ドップラーレーダ１０−１〜１０−ｎは、放射電波（本実施例では、Ｘバンド）のドップラー効果による反射波を解析して、動径方向のドップラー速度成分を出力する。反射の遅れを加味することで、動径線に沿った速度分布も計測可能である。各ドップラーレーダ１０−１〜１０−ｎのボリュームスキャンによる観測データ（ドップラー速度およびレーダ反射因子データ）は、各観測時の仰角、方位角及び動径距離の情報と共に、データネットワーク１２を介して本実施例の風向・風速算出装置１４に入力する。

風向・風速算出装置１４は、データネットワーク１２を介して各ドップラーレーダ１０−１〜１０−ｎから供給されるドップラー速度データ等から、風向・風速に関する数値気象予測データ（本実施例では、メソスケール気象予測モデル（ＭＳＭ）の出力データを採用した）と、地形データ（標高・地表面粗度データ）を参照して、対象三次元空間の風向・風速分布を算出する。

降水観測で用いられるドップラーレーダは、一般にパルスレーダであり、観測されるドップラー速度を折り返し補正する必要がある。すなわち、パルスレーダでは、測定できるドップラー速度には、パルス繰り返し周波数と波長で決定される最大速度（ナイキスト速度）が存在する。このナイキスト速度を超えるドップラー速度は折り返して表現される。実際のドップラー速度Ｖと、検出されたドップラー速度Ｖｏの間には、
Ｖnyq = ｆλ/４ （１）
Ｖ＝Ｖｏ＋２ｎVnyq （２）
という関係が成立する。但し、ナイキスト速度をＶnyqはナイキスト速度、ｆはパルス繰り返し周波数、λは電波の波長、ｎは整数である。

例えば、パルス繰り返し周波数ｆが２０００Ｈｚ、波長λが３．２ｃｍのドップラーレーダの場合、表現できるドップラー速度は−１６ｍ／ｓから＋１６ｍ／ｓである。計測されたドップラー速度Ｖｏが＋１２ｍ／ｓのとき、実際のドップラー速度Ｖは＋１２＋３２×ｎ（ｍ／ｓ）となり、即ち、．．．-５２、-２０、＋１２、＋４４、＋７６、．．．という不確定性が存在する。これは、ある信号を標本化するときに、サンプリング周波数の１／２を超える周波数は折り返されて表現されることに由来する。

そこで、風向・風速算出装置１４は、数値気象予測モデルの予報データを参照してドップラー速度の折り返しを補正するドップラー速度折り返し補正装置１６と、補正装置１６により補正されたドップラー速度データから、地形（標高・地表面粗度データ）を参照して、風向・風速の分布を計算する風向・風速分布計算装置１８とからなる。図２は、風向・風速算出装置１４の概略構成ブロック図を示す。

ハードディスク２０ａには、データネットワーク１２を介して入力する各ドップラーレーダ１０−１〜１０−ｎからのドップラー速度データが格納され、ハードディスク２０ｂには、各ドップラーレーダ１０−１〜１０−ｎからのレーダ反射因子データが格納される。また、ハードディスク２０ｃには、対象地域に対する数値気象モデルによる予報データが格納される。この予報データは、水平垂直の３次元空間で一定距離に位置する格子点上の値（Grid Point Value）で表現されている。ハードディスク２０ｄには、対象地域の標高・地表面粗度のデータが格納される。

ドップラー速度折り返し補正装置１６は、ハードディスク２０ａに格納されるドップラー速度データを、ハードディスク２０ｂに格納されるレーダ反射因子データ及びハードディスク２０ｃに格納される数値気象予報データを参照して、折り返し補正する。図３は、ドップラー速度折り返し補正装置１６における折り返し補正のフローを示す。

ドップラー速度折り返し補正装置１６は、ハードディスク２０ａからドップラー速度データを取り込むと共に、ハードディスク２０ｂに格納されるレーダ反射因子データを取り込み、一定の品質管理基準の下で、明らかな異常値を除外する（Ｓ１）。具体的には、観測されたドップラー速度には、地形による電波の反射や干渉波の影響によるノイズが含まれるので、これらの異常なデータを検出し削除する。地形は移動しないので、その影響を受けたドップラー速度はほぼゼロとなる。そこで、ゼロに近いドップラー速度データをまず削除する。次に、干渉波の影響を取り除くために、動径方向にドップラー速度データの移動平均偏差を計算し、その値がある設定した閾値を超える場合はそのドップラー速度データを削除する。

異常データ又はノイズデータを除外した後（Ｓ１）、本実施例では、複数の方法で折り返し補正を実行し、精度を高める。

まず、数値気象予報モデルとの比較により、折り返しを補正する（Ｓ２）。気象庁のような現業の気象機関は、数値気象モデル（例えば気象庁メソ数値モデル：ＭＳＭ）により計算された気温、気圧、風向・風速などの予報値を三次元格子データとしてリアルタイム配信している。ハードディスク２０ｃには、このような数値気象モデルの予報データが格納され、この数値気象モデルの予報データには、ある一定の時間間隔及び空間間隔で風速度データが含まれている。

ドップラー速度折り返し補正装置１６は、予報データの風速度データを、ドップラーレーダ１０−１〜１０−ｎによる測定位置（レーダ１０−１〜１０−ｎからの仰角・方位角・動径距離で記述される）に時空間内挿する。次に、内挿された風速度と、レーダ観測されたレーダ反射因子により推定される降水の落下速度とから、レーダ１０−１〜１０−ｎ動径方向の速度を計算する。そして、式（２）で与えられるドップラー速度Ｖがこの計算された動径方向の速度に一番近くなるｎを決定する。この操作を、品質管理された全てのドップラー速度データＶｏに対して実行し、それぞれに対するｎ、即ち、ドップラー速度Ｖを決定する。

本実施例では、折り返し補正に関して、複数の簡便な手法を順次組み合わせて使用している。そこで、各手法によりどの程度、折り返し補正が正しく行われたかを示す指標として、下記式に従いあるＰＰＩデータに関する下記式のＤ値を定義する。すなわち、

ｉは動径方向に関する添え字、ｊは方位角方向に関する添え字である。正しく折り返し補正が行われれば、ドップラー速度は空間的になめらかに変化する。そこで、動径方向・方位角方向に隣り合うドップラー速度データの差はナイキスト速度Ｖnyqよりも小さくなると考えられ、Ｄ値は減少する。

数値気象予報データとの比較により折り返し補正したドップラー速度データＶに対して、上記Ｄ値を算出する。この算出結果をＤ１とする。

次に、ドップラー速度折り返し補正装置１６は、ＶＡＤ（Velocity Azimuth Display）法による水平風との比較に基づき、折り返しを補正する（Ｓ３）。ＶＡＤ法とは、一つのＰＰＩ観測データから、観測領域を代表する水平風の高度分布を推定する手法である。まず、ドップラー速度折り返し補正装置１６は、ＶＡＤ法により、観測時刻における水平風速度の高度分布を計算する。次に、ドップラーレーダ１０−１〜１０−ｎによる測定位置の高度を求め、ＶＡＤ法で求めておいたその高度の水平風速度と、レーダ反射因子により推定される降水の落下速度とから、ドップラーレーダ１０−１〜１０−ｎの動径方向速度を計算する。ただし、仰角が１０度程度までならば、鉛直方向の風速成分はゼロとみなすことができる。

ステップＳ２における数値気象予報モデルとの比較と同様に、式（２）で、ここで計算された動径方向速度に一番近くなるｎを決定する。すなわち、ステップＳ２で折り返し補正されたドップラー速度データをＶｏとして式（２）に代入し、式（２）で得られるドップラー速度Ｖがここで計算された動径方向の速度に一番近くなるｎを決定する。この操作を、ステップＳ２で折り返し補正された全てのドップラー速度データに対して実行し、それぞれに対するｎ、即ち、ドップラー速度Ｖを決定する。

ステップＳ３での折り返し補正の結果に対して、式（３）によりＤ値を計算し、その結果をＤ２とする。Ｄ２がＤ１以下の場合、ステップＳ３による折り返し補正を有効とする。Ｄ２がＤ１よりも大きくなっている場合は、折り返し補正が正しく行われていないと判断し、ドップラー速度データをステップＳ３の適用前のものに戻し、Ｄ２にＤ１を代入する。すなわち、ステップＳ３の折り返し補正をキャンセルする。

次に、時間的な連続性に基づき、折り返しを補正する（Ｓ４）。上空の風速は時間に関して連続的に変化しており、その変動時間スケールよりもドップラーレーダの観測間隔が短い場合（例えば５分間程度）、直前の折り返し補正済みの速度データとの比較は有意である。そこで、先ず、直前の時間の同じ仰角における折り返し補正済みドップラー速度データを、ステップＳ３からの折り返し補正されたドップラー速度データの測定位置と同じ位置に空間的に内挿する。そして、ステップＳ３で折り返し補正されたドップラー速度データをＶｏとして式（２）に代入し、式（２）で得られるドップラー速度Ｖが、空間的に内挿されたドップラー速度に一番近くなるｎを決定する。この操作を、ステップＳ３で折り返し補正された全てのドップラー速度データに対して実行し、それぞれに対するｎ、即ち、ドップラー速度Ｖを決定する。

ステップＳ４での折り返し補正の結果に対して、式（３）によりＤ値を計算し、その結果をＤ３とする。Ｄ３がＤ２以下の場合、ステップＳ４による折り返し補正を有効とする。Ｄ３がＤ２よりも大きくなっている場合は、折り返し補正が正しく行われていないと判断し、ドップラー速度データをステップＳ４の適用前のものに戻し、Ｄ３にＤ２を代入する。すなわち、ステップＳ４の折り返し補正をキャンセルする。

最後に，空間連続性を考慮して、折り返しを補正する（Ｓ５）。すなわち、同時刻の隣接する位置の風速度は近い値を取ることを利用する。上空の風速は空間に関して連続的に変化しており、ドップラー速度に関しても、方位角・動径方向に隣り合うデータは連続に分布している。

まず、ステップＳ４の折り返し補正結果のうち、より補正結果の確からしいデータ点を抽出する。その抽出条件は例えば、ステップＳ２，Ｓ３，Ｓ４において、それぞれの補正済みドップラー速度と、それぞれの比較に用いられた速度データの差が常にある閾値（例えば１ｍ／ｓ）よりも小さいこととし、この条件に合致するデータ点を抽出する。次に、その確からしい点データと方位角・動径方向に隣り合う点のドップラー速度データに対して、確からしいデータ点のドップラー速度に一番近くなるｎを決定する。このように補正されたデータも確からしいデータ点として、この点に隣り合うデータを補正する。このような再帰的な処理により、全体のドップラー速度データの折り返しを補正する。

ステップＳ５での折り返し補正の結果に対して、式（３）によりＤ値を計算し、その結果をＤ４とする。Ｄ４がＤ３以下の場合、ステップＳ５による折り返し補正を有効とする。Ｄ４がＤ３よりも大きくなっている場合は、折り返し補正が正しく行われていないと判断し、ドップラー速度データをステップＳ５の適用前のものに戻す。

ドップラー速度折り返し補正装置１６は、このように複数の手法で折り返しを補正したドップラー速度データをハードディスク２２に格納する。

風向・風速分布計算装置１８は、ハードディスク２０ｂのレーダ反射因子データと、後述する標高・地表面データとを参照し、ハードディスク２２の折り返し補正済みドップラー速度データから上空の風向・風速分布を計算し、計算結果をハードディスク２４に格納する。

一般にレーダ観測では地球の曲率及び山岳・構造物等による電波の遮蔽により、地表面付近の観測が困難である。しかし、防災上は、地表面付近の風情報は非常に重要である。これに対し、本実施例では、対数則と呼ばれる風速の鉛直プロファイルモデルを仮定することにより、対流圏下層の指定した海抜高度Ｚ_１（たとえば、海抜高度１ｋｍ）における風向・風速と、地表面付近（たとえば、地上高度１０ｍ）における風向・風速を計算する。

計算に先立ち、海抜高度Ｚ_１おいて、緯度幅及び経度幅が一定で格子間隔が約５００ｍ程度となるような水平格子を定義する。次に、標高・地表面粗度データを、定義した格子点と同じ緯度・経度上に空間内挿する。ただし、地表面粗度データの空間内挿に関しては、その格子点の値ではなく、格子点から２５００ｍ程度風上の地点までの間で平均化した地表面粗度を使用する。たとえば、その格子点から風上方向に、０ｍ、１００ｍ、２００ｍ、…、２４００ｍ、２５００ｍというように１００ｍ置きに離れた地点の地表面粗度を空間内挿により計算し、それらを平均したものをその格子点における地表面粗度として使用する。これは、対流圏下層風速の鉛直プロファイルが直下の地表面状態のみでなく、風上側の地表面状態の影響を受けていることを考慮するためである。また、風上方向を定義する際に使用する風向・風速は、本実施例により直前に計算されたもの、即ち，ハードディスク２４の風向・風速分布データから読み出したデータと、計算毎に更新する。

定義した格子点上の風向・風速を以下のように計算する。ここで、各格子点上の風速のｘ軸（東西）、ｙ軸（南北）及びｚ軸（鉛直）各成分をｕ，ｖ，ｗとする。まず、各格子点の風速三成分（ｕ，ｖ，ｗ）と折り返し補正済みドップラー速度から計算するコストファンクションＪを定義し、このコストファンクションＪを最小とするような（ｕ、ｖ、ｗ）を求める。コストファンクションＪは、

と記述される。ここで、ｉは格子点を示す添え字である。ｊは格子点ｉを中心とする、ある大きさの楕円球（たとえば水平軸の半径が１．５ｋｍ、鉛直軸の半径が５００ｍ）の内側に分布するレーダ観測点を示す添え字である。ωは格子点ｉからレーダ観測点までの距離に依存する重みであり、距離が長くなると小さくなる関数で与えられる。α，β，γはそれぞれ、レーダビームがｘ軸（東西）、ｙ軸（南北）及びｚ軸（鉛直）方向となす角度を示す。ｆはレーダ反射因子により推定される降水の落下速度である。Ｋ_ｉは、格子点ｉの周りに分布するレーダ観測点におけるドップラー速度Ｖと格子点ｉにおける風速度の動径方向成分の差についての重み付き二乗平均である。

最終的に計算される風速三成分（ｕ，ｖ，ｗ）が、観測されたドップラー速度と矛盾しなければ、Ｋ_ｉの値は０に近づき、全ての格子点についてのＫ_ｉの総和であるコストファンクションＪの第一項目の値も０に近づくはずである。コストファンクションＪの第二項目は、空間に関するローパスフィルタの機能を持ち、減衰させる波長は係数Ｃによって調節される。

地上に設置したドップラーレーダによる観測の場合、得られるドップラー速度には、鉛直風成分の寄与が小さいので、コストファンクションＪを最小化させる（ｕ，ｖ，ｗ）を求める際に、鉛直風速成分ｗが不安定な挙動を示す。すなわち、鉛直風速成分ｗを確定するのが難しい。そこで、本実施例では、質量保存の法則と風速の対数則の仮定を適用し、鉛直風速成分ｗをｘ軸風速成分ｕとｙ軸風速成分ｖの関数として記述する。

地表面付近の海抜高度ｚにおける風速Ｕ（ｚ）、Ｖ（ｚ）は、風速の対数則により、

と表現できると仮定する。ただし、ｕ＊、ｖ＊は摩擦速度、ｋはカルマン定数、ｚ_ｈは地表面の海抜高度、ｚ_０は地表面粗度である。海抜高度ｚ_１における水平面内の風速を（ｕ、ｖ）とすると、式（１）、（２）から摩擦速度ｕ＊，ｖ＊を消去でき、

と変形できる。

式（４），（５）を地表面高度ｚ_ｈからｚ_１まで積分することにより、高度ｚ_１以下での質量フラックス（Ｆｘ，Ｆｙ）は、

となる。ただし、ρは大気密度である。ここでは、地表面付近の現象のみを扱っているので、大気密度ρは、高度によらず一定であると仮定している。Ｌは下記式

のように定義される。

式（９）及び式（１０）の質量フラックス（Ｆｘ，Ｆｙ）から、高度ｚ_１における鉛直風速成分ｗは、質量保存の法則により、

となる。ただし、ａは地球の半径、λは経度、φは緯度である。水平発散ｄは、

と表される。式（１２）により鉛直風速成分ｗは、水平速度成分（ｕ，ｖ）の関数として表すことができる。これを式（４）に代入することにより、コストファンクションＪは、ｕとｖの関数になる。

以上の計算により、最小にすべきコストファンクションＪが決定する。このコストファンクションＪが最小になる各格子点の風速成分（ｕ，ｖ）が決定されれば、式（１２）により鉛直風速成分ｗが求まり、式（５），（６）により高度ｚ_１以下の風速の鉛直プロファイルが決定される。

コストファンクションＪを最小にする各格子点の水平風速成分（ｕ、ｖ）の値を決定するには、非線形計画法で取り扱われる様々な解法を利用できる。たとえば共役勾配法の一種であるＰｏｌａｋ−Ｒｉｂｉｅｒｅ法を用いると、コストファンクションＪとそのｕ，ｖについての偏微分を与えることにより、各格子点における（ｕ，ｖ）の各値が求められる。このようにして得られた各格子点上の水平風速成分（ｕ，ｖ）の値、及び、これらから式（１２）により求めた鉛直方向風速成分（ｗ）は、ハードディスク２４に格納され、データ同化装置３０に供給される。

データ同化装置３０は、風向・風速算出装置１４からの風向・風速分布データ（風観測値データ）と、雲解像モデル（本実施例では、Cloud Resolving Storm Simulator（ＣＲｅＳＳ）を採用した）の予測結果である第一推定値を三次元変分法により同化し、観測値と物理的に整合する広範囲の風速場（これを客観解析データと呼ぶ）を作成する。同化手法として三次元変分法を用いることで、計算コストが甚大な四次元変分法やアンサンブルカルマンフィルタとは異なり、実時間における風速の客観解析データを提供することが可能となる。

図４は、データ同化装置３０の動作フローを示す。データ同化装置３０は、風向・風速算出装置１４により算出された風速場（ハードディスク）２４）のデータを風観測値として取り込み、レーダ反射因子データを降水粒子の観測値として取り込む。取り込んだ観測値を、その観測場所付近の雲解像数値モデル（ＣＲｅＳＳモデル）の格子点に内挿する（Ｓ２１）。

次に、リアルタイム運用されている雲解像数値モデルの予測データから、観測時刻に対応するデータを抽出し、それを第一推定値（観測以外の手段により推定される、各格子点の風向・風速、気圧、気温、水蒸気混合比、降水粒子の混合比）とする。

こうして、観測領域内の雲解像数値モデルの格子点には、観測値と第一推定値の２つの風データを持つことになる。以下の手順を繰り返すことで、観測値と第一推定値の両者にできるだけ近く、かつ、質量保存則をできるだけ満たす風の推定値を求める。すなわち、風の推定値と観測値との差分、風の推定値と第一推定値との差分、及び式（１４）で得られる収束発散項であるＤｖをすべての格子点について合計する（Ｓ２３）。すなわち、

ただし、ρは大気密度である。この合計値を最も小さくする風の推定値が最適な推定値である。総和が一定範囲に収束しない場合（Ｓ２４）、共役勾配法を使って風の推定値を更新し（Ｓ２５）、合計値を再計算する（Ｓ２３）。

合計値が最小値に収束したと判定したら（Ｓ２４）、データ同化を終了し、データ同化による風速場データ（客観解析データ）をハードディスク３２に保存する。ハードディスク３２に保存される風速場データは、ＣＲｅＳＳモデルで想定される地域の海抜１ｋｍ以下の空間における風向・風速を所定距離間隔の格子点上で示すデータである。この風速場の移動方向と速度を予測することで、強風の発生とその地域を予測できる。

予測装置３４は、データ同化装置３０により算出された現在及び過去の一定時間内の風速場の時間変化から風速場の移動ベクトル（時間変化）を推定し、１時間程度先までの範囲で５分毎の、地上での強風の発生とその地域を推定する。なお、水平解像度を５００ｍメッシュとし、更新間隔を５分とした。図５は、予測装置３４の移動ベクトル推定装置３６と上空風速場予測装置３８の動作フローを示す。

移動ベクトル推定装置３６が、ハードディスク３２の５００ｍメッシュの現在及び過去の一定時間内の風速場データから、強風域の移動方向と移動速度（以降、移動ベクトルと呼ぶ）を推定する。なお、本実施例の目的が強風域の予測なので、ハードディスク３２から現在及び過去３時刻の風速場データを読み込み（Ｓ３１），その内から、一定以上の強風域、ここでは１５ｍ／ｓ以上を抽出する（Ｓ３２）。抽出した風速場データに対し、移流モデルを用いて移動ベクトルを算出する（Ｓ３３）。強風域と予測されなかった領域の値は、不定（風速は弱いかもしれないし強いかもしれない）として扱う。

本実施例では、移動ベクトルの推定に移流モデルを採用する（Ｓ３３）。移流モデルの概念を簡単に説明する。風速をＺ（ｘ，ｙ，ｔ）のように位置（ｘ，ｙ）と時刻ｔの関数で表したとき、Ｚ（ｘ，ｙ，ｔ）は、時刻を固定すると１つの曲面を表す。風速予測は、この曲面の変化を予測することに相当する。

そこで、風速Ｚの変化を、

と表す。式（１５）はいわば質量保存の式であり、ある微小面積内でのＺの時間変化と、その微小面積に周囲から流れ込み又は周囲に流れ出す量との和が、その場所でのＺの増減と等しいことを意味している。ここで、ｕはＺの東西方向の移動速度を示し、ｖはＺの南北方向の移動速度を示す。ｕ，ｖは、Ｚの移動方向と移動速度を表すので、これを移動ベクトルと呼ぶ。また、ｗはＺの発達・衰弱を表す。

（１５）式における（ｕ，ｖ，ｗ）を、パラメータＣ_１〜Ｃ_９を使って、位置座標の一次式で表すことができるものと仮定すると、下記のように表される。

式（１６）により、風速Ｚの平行移動だけでなく、変形も表現できる。

式（１６）のパラメータＣ_１〜Ｃ_９を同定できれば、風速Ｚの移動ベクトルと発達・衰弱量を推定できる。そこで、連続する複数時刻分の風速データ（Ｚ）を使用して、点（ｘ，ｙ）、時刻ｔ_ｋにおける（１５）式の残差（左辺と右辺の差）をＶ_ｉｊｋとおくと、

が得られる。ここで、ｉ，ｊ，ｋは三次元の格子点を示す添え字である。

パラメータＣ_１〜Ｃ_９が正しく同定されれば、残差Ｖ_ｉｊｋはゼロとなるはずである。これは、下記式

で与えられる評価値Ｊ_ｃを最小にすることに相当する。従って、式（１８）を最小にするように、Ｃ_１〜Ｃ_９を決めれば良く、具体的には、式（１８）のＪ_ｃをＣ_１〜Ｃ_９で偏微分して０とおいて得られる連立一次方程式を解くことにより求めることができる。

このモデルでは、風速域の回転まで表現することができる。ただし、精度上問題がある例も見られるので、予測範囲をいくつかの小領域に分け、小領域ごとに単一の移動ベクトルを算出した後、それらを内挿してメッシュ単位の移動ベクトルに置き換えて、場所による風速域の移動方向の変化を表現する。小領域毎の移動ベクトルは、パラメータＣ_１〜Ｃ_９のうち、Ｃ_３およびＣ_６以外はすべて０とする。これは、各小領域では風速域の移動が平行移動であると仮定していることに相当する。

上空風速場予測装置３８は、移動ベクトル推定装置により推定された移動ベクトルが１時間先まで継続すると仮定して、ハードディスク３２の現在の風速場データの強風域を外挿する（Ｓ３４）。これにより、１時間先までの強風域の位置又は移動を予測する。

上空風速場予測装置３８により予測される強風域は、上空（例えば、高度１０００ｍ）での強風域である。そこで、地上強風域予測装置４０が、上空風速場予測装置３８により予測される各時刻の上空強風域に対し地上強風域を予測する。

対象地域の風向・風速の観測値が存在する場合には、現時刻の地上風速観測値と上空風速値の差を、上空風速場予測装置３８により得られる上空風速予測値に加算し、その加算結果を地上風速予測値とする。

対象地域の観測値が存在しない場合、対数則と呼ばれる風速の鉛直プロファイルに基づき、上空風速場予測装置３８により得られる上空風速予測値から地上風速値を予測する。

地上強風域予測装置４０は、このようにして得られた各時刻の地上風向・風速データをハードディスク４２に格納する。強風域の表示を目的としているので、ハードディスク４２に格納される地上風向・風速データは、上空で強風になっている地域のデータである。

なお、移動ベクトル推定装置３６の推定単位、並びに、上空風速場予測装置３８及び地上強風域予測装置４０の予測単位を、予測対象地域を分割した小領域とすることで、計算量を削減でき、計算時間を短縮できる。

表示装置４４は、ハードディスク４２に格納される風向・風速分布情報を、グラフィカル・ユーザ・インターフェースの画面上に、強度を色で区分して表示するグラフィックとして表示する。

本実施例により、風、特に強風を局所的に予測できる。得られる風の予測情報を基にして、降雨の発達・衰弱を考慮できるようになるので、降雨の短時間予測精度も向上する。

本実施例は、主として、コンピュータプログラムにより実現されるが、その機能の一部を専用ハードウエアに置換しても同様の作用効果を奏することができることは明らかである。また、単一のコンピュータ上で動作するコンピュータプログラムのみならず、多数のコンピュータ上でそれぞれ動作するコンピュータプログラムを協調動作させることでも、同様の作用効果を奏することができる。これらのいずれの構成も、本発明の技術的範囲に属するものである。

特定の説明用の実施例を参照して本発明を説明したが、特許請求の範囲に規定される本発明の技術的範囲を逸脱しないで、上述の実施例に種々の変更・修整を施しうることは、本発明の属する分野の技術者にとって自明であり、このような変更・修整も本発明の技術的範囲に含まれる。

本発明の一実施例の概略構成ブロック図である。 風向・風速算出装置１４の概略機能ブロック図である。 ドップラー速度折り返し補正装置１６における折り返し補正のフローを示す。 データ同化装置３０の動作フローである。 移動ベクトル推定装置３６及び上空風速場予測装置３８の動作フローである。

符号の説明

１０−１〜１０−ｎ：ドップラーレーダ
１２：データネットワーク
１４：風向・風速算出装置
１６：ドップラー速度折り返し補正装置
１８：風向・風速分布計算装置
２０ａ：ハードディスク（ドップラー速度データ）
２０ｂ：ハードディスク（レーダ反射因子データ）
２０ｃ：ハードディスク（数値気象モデル予報データ）
２０ｄ：ハードディスク（標高・地表面粗度データ）
２２：ハードディスク（折り返し補正済みドップラー速度データ）
２４：ハードディスク（風向・風速分布データ）
３０：データ同化装置
３２：ハードディスク（風速場データ）
３４：予測装置
３６：移動ベクトル推定装置
３８：上空風速場予測装置
４０：地上強風域予測装置
４２：ハードディスク（風向・風速分布予測データ）
４４：表示装置

Claims (8)

1. 複数のドップラーレーダにより風観測されたドップラー速度データを記憶する記憶手段（２０ａ）と、
   当該ドップラー速度データを折り返し補正する折り返し補正手段（１６）と、
   当該折り返し補正手段により補正されたドップラー速度データから、地形データを参照して風向・風速分布データを計算する風向・風速分布計算手段（１８）と、
   当該風向・風速分布データを数値気象モデルにデータ同化し、所定格子間隔の風速場データを生成するデータ同化手段（３０）と、
   所定時間内の複数の風速場データから、風の移動ベクトルを推定する移動ベクトル推定手段（３６）と、
   当該移動ベクトルに従い、現在の風速場データを時間方向に外挿して上空風速場を予測する上空風速場予測手段（３８）と、
   当該上空風速場予測手段の予測結果に従い所定の地上風域を予測する地上風域予測手段（４０）
   とを具備することを特徴とする風予測装置。
2. 当該折り返し補正手段が、数値気象予報データとの比較による第１の折り返し補正、水平風との比較による第２の折り返し補正、時間方向での比較による第３の折り返し補正、及び空間連続性を考慮した第４の折り返し補正の内の２以上を実行することを特徴とする請求項１に記載の風予測装置。
3. 当該記憶手段が、当該複数のドップラーレーダからのレーダ反射因子データを記憶し、
   当該当該折り返し補正手段の、当該第１及び第２の折り返し補正において、当該レーダ反射因子を参照する
   ことを特徴とする請求項２に記載の風予測装置。
4. 当該地上風域予測手段は、当該上空風速場予測手段の予測結果に従い地上の強風域を予測する地上強風域予測手段（４０）であることを特徴とする請求項１乃至３の何れか１項に記載の風予測装置。
5. 複数のドップラーレーダにより風観測されたドップラー速度データをコンピュータの記憶手段に格納する機能と、
   当該コンピュータに、当該ドップラー速度データを折り返し補正させる折り返し補正機能（１６）と、
   当該コンピュータに、当該折り返し補正機能により補正されたドップラー速度データから、地形データを参照して風向・風速分布データを計算させる風向・風速分布計算機能（１８）と、
   当該コンピュータに、当該風向・風速分布データを数値気象モデルにデータ同化させ、所定格子間隔の風速場データを生成させるデータ同化機能（３０）と、
   当該コンピュータに、所定時間内の複数の風速場データから、風の移動ベクトルを推定させる移動ベクトル推定機能（３６）と、
   当該コンピュータに、当該移動ベクトルに従い、現在の風速場データを時間方向に外挿して上空風速場を予測させる上空風速場予測機能（３８）と、
   当該コンピュータに、当該上空風速場予測機能の予測結果に従い所定の地上風域を予測させる地上風域予測機能（４０）
   とを具備することを特徴とする風予測プログラム。
6. 当該折り返し補正機能が、数値気象予報データとの比較による第１の折り返し補正機能、水平風との比較による第２の折り返し補正機能、時間方向での比較による第３の折り返し補正機能、及び空間連続性を考慮した第４の折り返し補正機能の内の２以上を有することを特徴とする請求項５に記載の風予測プログラム。
7. 当該記憶手段が、当該複数のドップラーレーダからのレーダ反射因子データを記憶し、
   当該第１及び第２の折り返し補正機能が当該レーダ反射因子を参照する
   ことを特徴とする請求項６に記載の風予測プログラム。
8. 当該地上風域予測機能は、当該コンピュータに当該上空風速場予測機能の予測結果に従い地上の強風域を予測させる地上強風域予測機能（４０）であることを特徴とする請求項５乃至７の何れか１項に記載の風予測プログラム。

Images