(57)【特許請求の範囲】
【請求項１】
雨量観測手段によって降雨観測領域を観測して得られ、該降雨観測領域を格子状に分割したときの時系列で連続する該格子毎の降雨強度に関するデータを取得する時系列降雨分布入力工程と、
前記時系列降雨分布入力工程で取得する格子毎の降雨強度に関するデータを、該格子及び該格子に隣接する格子においてゼロでない降雨強度を有するものの降雨強度の平均値に基づいて補正する降雨分布補正工程と、
前記降雨分布補正工程によって補正された格子毎の位置並びに降雨強度に関するデータから移流方程式に基づいて格子毎の移流ベクトルを推定する移流ベクトル推定工程と、
前記時系列降雨分布入力工程で取得する格子毎の降雨強度に関するデータを前記移流ベクトル推定工程で推定する格子毎の移流ベクトルに基づいて移動して得られる降雨分布移動予測データの生成を行なう降雨移動予測工程と、
を含むことを特徴とするコンピュータによる降雨予測方法。
【請求項２】
前記降雨分布補正工程は、前記時系列降雨分布入力工程で取得する格子毎の降雨強度に関するデータを降雨強度によって分類し、該格子及び該格子に隣接する格子においてゼロでない降雨強度を有するものの降雨強度に対して該分類の結果に応じた重み付けを行い、該重み付けのされた降雨強度の平均値に基づいて該データを補正することを特徴とする請求項１に記載のコンピュータによる降雨予測方法。
【請求項３】
前記移流ベクトル推定工程は、最新時刻の降雨強度を含む２以上の観測時刻における同一格子位置の降雨強度に基づいて格子毎の移流ベクトルを推定することを特徴とする請求項１に記載のコンピュータによる降雨予測方法。
【請求項４】
前記移流ベクトル推定工程は、最新時刻の降雨強度を含む２以上の観測時刻における同一格子位置の降雨強度の空間差分値の平均値を格子毎の降雨強度の空間差分値とし、該空間差分値を偏微分項を差分化した移流方程式に代入し、該移流方程式に基づいて格子毎の移流ベクトルを推定することを特徴とする請求項３に記載のコンピュータによる降雨予測方法。
【請求項５】
雨量観測手段によって降雨観測領域を観測して得られ、該降雨観測領域を格子状に分割したときの時系列で連続する該格子毎の降雨強度に関するデータを取得させる時系列降雨分布入力手順と、
前記時系列降雨分布入力手順で取得させる格子毎の降雨強度に関するデータを、該格子及び該格子に隣接する格子においてゼロでない降雨強度を有するものの降雨強度の平均値に基づいて補正する降雨分布補正手順と、
前記降雨分布補正手順によって補正された格子毎の位置並びに降雨強度に関するデータから移流方程式に基づいて格子毎の移流ベクトルを推定させる移流ベクトル推定手順と、
前記時系列降雨分布入力手順で取得させる格子毎の降雨強度に関するデータを前記移流ベクトル推定手順で推定させる格子毎の移流ベクトルに基づいて移動して得られる降雨分布移動予測データの生成を行なわせる降雨移動予測手順と、
をコンピュータに実行させるプログラムを記録したことを特徴とするコンピュータ読み取り可能な記録媒体。
【請求項６】
前記降雨分布補正手順は、前記時系列降雨分布入力工程で取得させる格子毎の降雨強度に関するデータを降雨強度によって分類し、該格子及び該格子に隣接する格子においてゼロでない降雨強度を有するものの降雨強度に対して該分類の結果に応じた重み付けを行い、該重み付けのされた降雨強度の平均値に基づいて該データを補正することを特徴とする請求項５に記載のコンピュータ読み取り可能な記録媒体。
【請求項７】
前記移流ベクトル推定手順は、最新時刻の降雨強度を含む２以上の観測時刻における同一格子位置の降雨強度に基づいて格子毎の移流ベクトルを推定させることを特徴とする請求項５に記載のコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、コンピュータによる降雨予測の精度を向上させる方法に関する。
【０００２】
【従来の技術】
従来、デジタルレーダ雨量計等で得られる降雨分布の時系列データを用いた短時間降雨予測の代表的な手法として、雨域追跡法と移流モデル法とがある。
【０００３】
雨域追跡法に基づく降雨予測の方法は、例えば「レーダ雨量計による短時間降雨予測に関する調査報告書」（建設省土木研究所河川部水文研究室、土木研究所資料第２４０６号）に述べられている。
【０００４】
雨域追跡法による予測方法の一例を図９に示す。同図に基づいて雨域追跡法を簡単に説明する。
時系列降雨分布入力工程１では、降雨観測領域を格子に分割し、その格子毎に観測される連続する時系列降雨データを取得する。ここで入力されるデータは、各格子毎での単位時間あたりの降雨量（降雨強度）に関するデータであり、デジタルレーダ雨量計等の雨量観測手段によって観測されるものである。
【０００５】
移動ベクトル演算工程２では、まず、時系列降雨分布入力工程１で取得した時系列降雨データを、例えば予め設定される雨量の閾値に対する多少によって２値化し、時系列の雨域パターンを得る。次に、ある時刻の雨域パターンを各方向・各距離に移動させた場合の、後の時刻の雨域パターンとの相互相関係数を計算する。そして、この相互相関係数の計算値が最大の値を示すときの雨域パターンの移動の方向・距離を、計算に使用した２つの雨域パターンの観測時刻の間における雨域の移動ベクトルとみなし、これを降雨移動予測工程３へ渡す。
【０００６】
降雨移動予測工程３では、前工程で得られた移動ベクトルに従って雨域が将来も移動するものとし、時系列降雨分布入力工程１で取得した降雨データをこの移動ベクトルに基づいて移動させ、将来の降雨の移動予測データを生成する。
【０００７】
出力工程４では、降雨移動予測工程３で作成した雨域移動の予測データを表示装置等の出力手段に出力する。
雨域追跡法による降雨予測の範疇に含まれる事例としては、次のような報告あるいは実用化例がある。
・連続する時系列降雨パターンを平行移動と仮定し、相互相関関数法により相関係数を計算し、その係数の最大となるものを雨域の移動であるとして移動ベクトルを求め、外挿によって将来の雨域移動を予測する（Austin, G. L. and A. Bellon: Very-short-range forecasting of precipitation by the objective extrapolation of radar and satellite data, Nowcasting(ed. By A. Browning), Academic Press, New York(1982) ）。
・連続する２つの時刻の降水量分布の類似度を計算し、その値が最大となるときの移動方向を移動ベクトルとする。そして、異なる地域の複数の移動ベクトルを用いて降雨予測をする（大林正典、高橋邦夫：気象庁の降水短時間予報、水文・水資源学会誌  Vol.5,No.2(1992)）。
・連続する２つの時系列降雨パターンの重心を求め、その移動方向を移動ベクトルとし、外挿によって将来の雨域移動を予測する。
・連続する２つの時系列降雨パターンの曲面形状を計算し、その特徴点同士の照合により確からしい組み合わせを求めて移動ベクトルとし、外挿によって将来の雨域移動を予測する（森田豊久、田所英之：レーダ雨量計を用いた降雨量予測方式の提案、第３１回下水道研究発表会講演集(1994)）。
【０００８】
次に、移流モデル法について説明する。図１０は従来の移流モデル法による予測方法を示す図である。同図において、時系列降雨分布入力工程１、出力工程４は図９に示す雨域追跡法による予測方法と同一の工程であり、説明は省略する。
【０００９】
移流モデルパラメータ演算工程５は、時系列降雨分布入力工程１で取得した時系列降雨分布データを用いて移流方程式のパラメータを線形最小二乗推定により算出する。この工程について更に説明を加える。
【００１０】
位置（ｘ，ｙ）の観測格子点における時刻ｔの降雨強度をｚ（ｘ，ｙ，ｔ）とし、ｚの空間時間変化が次の移流方程式で表せるとする。
∂ｚ／∂ｔ＋ｕ（∂ｚ／∂ｘ）＋ｖ（∂ｚ／∂ｙ）＝ｗ            （１）
但し、（ｕ，ｖ）は移流ベクトル、ｗは発達衰弱項を表す。ここで
ｕ＝α1 ｘ＋α2 ｙ＋α3 
ｖ＝α4 ｘ＋α5 ｙ＋α6                                       （２）
ｕ＝α7 ｘ＋α8 ｙ＋α9 
と仮定し、式中の移流モデルパラメータα1 、……、α9 を、時系列降雨分布入力工程１で取得する時系列降雨分布データを用いて、以下の線形最小二乗推定で求める。
【００１１】
まず、ｘ方向にΔｘ間隔、ｙ方向にΔｙ間隔でとったＭ×Ｎ個のそれぞれの格子点（ｘi ，ｙj ）、ｉ＝１、…、Ｍ、ｊ＝１、…、Ｎについて、時刻ｔk における移流方程式を差分化する。ここで、時間差分幅をΔｔとし、また、簡単のためにｚ（ｘi ，ｙj ，ｔk ）をｚijk と表すと、
（∂ｚ／∂ｔ）ijk ＝（ｚij(k+1) −ｚijk ）／Δｔ
（∂ｚ／∂ｘ）ijk ＝（ｚ(i+1)jk −ｚ(i-1)jk ）／２Δｘ        （３）
（∂ｚ／∂ｙ）ijk ＝（ｚi(j+1)k −ｚi(j-1)k ）／２Δｙ
と近似できる。なお、時刻ｔk+1 までのｚのデータは前行程で入力されている既知の値とする。（３）式を（１）式に代入することにより、次式を得る。
【００１２】
（α1 ｘi ＋α2 ｙj ＋α3 ）（ｚ(i+1)jk −ｚ(i-1)jk ）／２Δｘ
＋（α4 ｘi ＋α5 ｙj ＋α6 ）（ｚi(j+1)k −ｚi(j-1)k ）／２Δｙ
−（α7 ｘi ＋α8 ｙj ＋α9 ）
＝−（ｚij(k+1) −ｚijk ）／Δｔ−ｄij                          （４）
ここで、ｄijは残差を表すために導入したものであって、この値はこれまでの仮定が正しければ本来ゼロになるべきものであり、よって、この残差の二乗和を最小にするようにα1 、……、α9 を定めるようにする。このためには
【００１３】
【数１】
0003727762.tif 000002
【００１４】
をα1 、……、α9 で偏微分して得られるα1 、……、α9 を未知数とする連立一次方程式を解けばよい。
移流モデルパラメータ演算工程５では、以上の手順に従って移流モデルパラメータα1 、……、α9 を算出する。
【００１５】
図１０の説明に戻る。移流ベクトル演算工程６は、移流モデルパラメータ演算工程５で算出される移流方程式のパラメータα1 、……、α9 を（２）式に代入することで、各観測格子毎の移流ベクトル（ｕ，ｖ）と発達衰弱項ｗとを算出する。以降、移流モデルパラメータ演算工程５と移流ベクトル演算工程６とを総称して移流ベクトル推定工程と呼ぶ。
【００１６】
降雨移動予測工程３では、図９の雨域追跡法での工程と同様に、将来の降雨の移動予測データを作成するので同一の名称を付している。しかしながら、移流モデル法では、前工程で得られた各観測格子毎の移流ベクトルと発達衰弱項に基づいて降雨分布が変動するものとしており、時系列降雨分布入力工程１で取得した最新の時刻の降雨データとこの移流ベクトル及び発達衰弱項とから、各観測格子の将来の降雨強度を算出し、この算出データを総合して将来の雨域移動の予測データを作成している点で、雨域追跡法における工程とは若干の相違がある。
【００１７】
移流モデル法については、「デジタルレーダーデータを利用した実時間降雨・流量予測システムの開発」（高棹琢馬、昭和６２年度科学研究費補助金（試験研究（１））研究成果報告書（昭和６３年３月））で詳細に説明されている。
【００１８】
【発明が解決しようとする課題】
雨域追跡法（及びそれに類する方法）は、予測時刻までの時間内には雨域の形状は変化しないという前提の下で、１もしくは少数で代表させた移動ベクトルに従って雨域を移動させることによって降雨予測を行なうために、雨域の形状変化（回転、歪み、膨張、収縮もしくはこれらの複合した変化等）や降雨量の増加減少を表すことができず、降雨予測の精度低下の原因となっている。このため、雨域内での降雨分布の時間変化を予測する手法を付加する等して予測精度の向上を図ることが行われているが、その結果、予測の方法は複雑化し、予測に要する演算量の増大を招いていた。
【００１９】
一方、移流モデル法においては、空間差分を用いてパラメータ推定を行なうため、観測格子のｘ、ｙ方向の間隔に比較して、連続する時系列降雨分布のパターンの変化は少ないことが望ましい。従って、雨域の移動が速い場合、あるいは空間的予測精度の向上のために格子間隔を小さくした場合において、観測時刻の連続する降雨分布パターン間での重なり部分が少なくなると、降雨分布の移動が平行移動のみであってもパラメータの精度は低下し、場合によってはパラメータの導出が不可能となることもある。また、降雨パターンが回転する等して大きく変化した場合には、その変化に追従できず、降雨パターンの距離的移動の大小とは無関係に予測精度が低下する問題も生じていた。更に、降雨データを２値化データもしくは降雨強度の段階的なレベル情報として与える場合には、本来ゼロでない値を取るはずの空間差分がゼロとなってしまい、パラメータの導出が困難となり、移流モデルを適用できなくなることがある問題もあった。
【００２０】
以上の問題点を鑑み、本発明は、これらの雨域の移動状況や降雨データの観測体制等によって引き起こされる予測精度の低下を改善する降雨予測方法を提供することを目的としている。
【００２１】
【課題を解決するための手段】
前述の目的を達成するために、雨量観測手段によって降雨観測領域を観測して得られ、該降雨観測領域を格子状に分割したときの時系列で連続する該格子毎の降雨強度に関するデータを取得する時系列降雨分布入力工程と、格子毎の位置並びに降雨強度に関するデータから移流方程式に基づいて格子毎の移流ベクトルを推定する移流ベクトル推定工程と、前記時系列降雨分布入力工程で取得する格子毎の降雨強度に関するデータを前記移流ベクトル推定工程で推定する格子毎の移流ベクトルに基づいて移動して得られる降雨分布移動予測データの生成を行なう降雨移動予測工程とを含むコンピュータによる降雨予測方法において、本発明では、前記時系列降雨分布入力工程で取得する格子毎の降雨強度に関するデータを降雨強度によって分類し、予め設定する強降雨地域に該格子が含まれるときには主に降雨強度の観測値の測定誤差を平滑化する補正を行ない、予め設定する弱降雨地域に該格子が含まれるときには主に降雨域の裾を拡大（後述）する補正を行なう降雨分布補正工程を更に含み、前記移流ベクトル推定工程が、前記降雨分布補正工程によって補正された格子毎の降雨強度に関するデータから格子毎の移流ベクトルを推定することを特徴とし、これらの工程により、強降雨地域ではデータの測定誤差に起因している降雨予測の精度低下を改善し、弱降雨地域では時系列降雨分布間での降雨分布パターンの重なりを拡大させて、移流ベクトル等の算出に必要なパラメータ推定の精度を向上させる。
【００２２】
また、前記コンピュータによる降雨予測方法において、本発明では、前記移流ベクトル推定工程が、最新時刻の降雨強度を含む２以上の観測時刻における同一格子位置の降雨強度に基づいて格子毎の移流ベクトルを推定することを特徴とし、これらの工程により、降雨分布パターンの回転等による急激な変化に対する、移流ベクトルの推定精度の向上を図ることができる。
【００２３】
また、前記コンピュータによる降雨予測方法において、本発明では、前記移流ベクトル推定工程が、異なる２時刻の格子毎の降雨強度に関するデータを入力し、該時刻間における雨域の移動を代表して示す基本移流ベクトルを得る基本移流ベクトル演算工程と、前記基本移流ベクトル演算工程で使用するデータの中で移動前の格子毎の降雨強度に関するデータを該基本移流ベクトル演算工程で得られた基本移流ベクトルに基づいて移動して得られる格子毎の仮想降雨分布データの生成を行なう仮想降雨分布演算工程と、前記基本移流ベクトル演算工程で使用する移動後の格子毎の降雨強度に関するデータと前記仮想降雨分布演算工程で生成する格子毎の仮想降雨分布データとの両者の同一格子位置における降雨強度に関するデータを代入する前記移流方程式に基づいて格子毎の補正移流ベクトルを推定する補正移流ベクトル演算工程と、基本移流ベクトル演算工程で得られる基本移流ベクトルと前記補正移流ベクトル演算工程で得られる格子毎の補正移流ベクトルとを格子毎に合成して得られる合成移流ベクトルを生成する合成移流ベクトル演算工程とを含み、前記移流ベクトル推定工程が、前記合成移流ベクトル演算工程で生成する合成移流ベクトルを該移流ベクトル推定工程の推定するベクトルとすることを特徴とし、これらの工程によっても、降雨分布パターンの回転等による急激な変化に対する、移流ベクトルの推定精度の向上を図ることができる。
【００２４】
また、前記コンピュータによる降雨予測方法において、本発明では、前記移流ベクトル推定工程が、該移流ベクトル推定工程の推定する格子毎のベクトルを代表する実績単位時間移流ベクトルを逐次記憶する実績単位時間移流ベクトル記憶工程と、前記実績単位時間移流ベクトル記憶工程で逐次記憶する実績単位時間移流ベクトルを代表する基本移流ベクトルを得る基本移流ベクトル演算工程と、格子毎の降雨強度に関するデータを該基本移流ベクトル演算工程で得られる基本移流ベクトルに基づいて移動して得られる格子毎の仮想降雨分布データの生成を行なう仮想降雨分布演算工程と、前記仮想降雨分布演算工程で使用した格子毎の降雨強度に関するデータと該仮想降雨分布演算工程で生成する格子毎の仮想降雨分布データとの両者の同一格子位置における降雨強度に関するデータを代入する前記移流方程式に基づいて格子毎の補正移流ベクトルを推定する補正移流ベクトル演算工程と、前記基本移流ベクトル演算工程で得られる基本移流ベクトルと前記補正移流ベクトル演算工程で得られる格子毎の補正移流ベクトルとを格子毎に合成して得られる合成移流ベクトルを生成する合成移流ベクトル演算工程とを含み、前記移流ベクトル推定工程が、前記合成移流ベクトル演算工程で生成する合成移流ベクトルを該移流ベクトル推定工程の推定するベクトルとすることを特徴とし、これらの工程により、基本移流ベクトルの算出を簡略化して、算出に必要な演算量を減少させ、予測処理の短時間化が可能となる。
【００２５】
【発明の実施の形態】
以下、本発明の実施の形態を図面に基づいて説明する。
図１は本発明の第一の実施の形態を示す図である。同図において、図１０と同様の工程には同一の符号を付している。本実施例は、従来技術である移流モデル法の移流ベクトル推定工程に、補正演算パラメータ入力工程７と降雨分布補正工程８とが加わった構成である。
【００２６】
補正演算パラメータ入力工程７では、キーボード等の入力手段に入力される降雨分布補正工程８での補正演算のために必要なＲmin 、β、γ1 、γ2 、指定回数の５つのパラメータを取得する。これらのパラメータについては後述する。
【００２７】
降雨分布補正工程８では、補正演算パラメータ入力工程７で取得したパラメータを用いて、時系列降雨分布入力工程１で取得した時系列降雨分布データの補正を行ない、その結果を移流モデルパラメータ演算工程５に渡す。この補正の方法について更に述べる。
【００２８】
図２は降雨分布補正工程８で行なう処理のフローチャートである。同図に基づいて説明を進める。
まず、ステップＳ１では、時系列降雨分布入力工程１で取得した時系列降雨分布データ（但し、ステップＳ９からの繰り返しの途中のときはその補正中の値）から、時刻ｔk における降雨観測範囲の全格子の中での最大降雨強度値Ｒmax を探索し、ステップＳ２へ進む。
【００２９】
ステップＳ２では、β・Ｒmax と降雨強度域判定最小値Ｒmin との大小比較を行なう。ここで、βは後述のステップＳ８で処理する補正演算において、降雨地域を弱い降雨地域（以降弱降雨地域と呼ぶ）とそれよりも強い降雨のある地域（以降強降雨地域と呼ぶ）とに区別するパラメータであり、０．０〜１．０の範囲の値をとり得る。この比較の結果、β・Ｒmax の値がＲmin 以上であるときはステップＳ３へ進み、それ以外のときはステップＳ４へ進む。
【００３０】
ステップＳ３では、Ｒｃ（強降雨地域と弱降雨地域とを区別する降雨強度の閾値）にβ・Ｒmax を代入してステップＳ５へ進む。
ステップＳ４では、ＲｃをＲmin としてステップＳ５へ進む。このステップＳ４は降雨観測範囲において、βで決定される弱降雨地域は存在しない場合のＲｃを設定している。
【００３１】
ステップＳ５では、Ｒ0ij とＲｃとの比較を行なう。ここで、Ｒ0ij とは観測格子（ｉ，ｊ）における降雨強度の補正値（この値をＲijとする）の前回値である。比較の結果、Ｒ0ij がＲｃ以上の場合には、この観測格子はβで決定される強降雨地域に含まれるとしてステップＳ６へ進み、そうでない場合には、この観測格子は弱降雨地域に含まれるとしてステップＳ７へ進む。
【００３２】
ステップＳ６では、γの値をγ1 としてステップＳ８へ進む。ここで、γはステップＳ８で処理する補正演算で使用するパラメータであり、０．０〜１．０の範囲の値をとり得る。また、γ1 はステップＳ８で処理する補正演算を観測値の測定誤差の平滑化を重視する場合に適するγの値である。すなわち、βで決定される強降雨地域に含まれる場合には、測定誤差の平滑化を重視した補正演算を行なって、データの測定誤差に起因する降雨予測の低下を改善するのである。
【００３３】
ステップＳ７では、γの値をγ2 としてステップＳ８へ進む。ここで、γ2 はステップＳ８で処理する補正演算を、降雨域の境界をそれに隣接する降雨のない格子点にまで拡大させて（このことを「降雨域の裾の拡大」と呼ぶ）時系列降雨パターン間の重なりを増加させることを重視する場合に適するγの値である。すなわち、βで決定される弱降雨地域に含まれる場合には、降雨域の裾の拡大を重視した補正演算を行なって、移流ベクトル等の算出に必要なパラメータ推定の精度向上を目指すのである。
【００３４】
ステップＳ８では、今までのステップで決定したγを使用して、以下の式に従って計算を行ない、Ｒijを算出する。
【００３５】
【数２】
0003727762.tif 000003
【００３６】
図３にＲijの位置関係を示す。つまり、（６）式においてｎ＝５とすれば、Ｒijは格子（ｉ，ｊ）における降雨強度補正値の前回値と、それ自身及びその格子のｘ方向若しくはｙ方向に隣接する格子の降雨強度補正値の前回値の平均値との、γに基づく重み付け加算によって算出する。ここで、γの値が１．０に近づけばＲijは隣接格子との平均値に近づき、逆にγの値が０に近づけばＲijは自身の前回値に近づく。
【００３７】
本実施の形態では、Ｒijの補正について、移流ベクトル等の算出に必要なパラメータ推定の精度を更に向上させるために、前述のｎ＝５に固定しない補正を行なう。図４は本実施の形態におけるＲijの算出法を示すフローチャートである。同図を説明する。
【００３８】
ステップＳ０８１では、図３に示す５つの降雨強度補正値の前回値の中で、ゼロでない値を有しているものの個数をｎに代入して、ステップＳ０８２へ進む。ステップＳ０８２では、ｎとゼロとを比較し、ゼロでなければステップＳ０８３へ進み、ゼロならばステップＳ０８４へ進む。
【００３９】
ステップＳ０８３では、ステップＳ０８１で与えられたｎの値と図２のフローチャートで決定されたγとを前述の（６）式に代入してＲijを算出して処理を終了する。つまり、隣接格子の降雨強度補正値の前回値の平均の算出に当たり、特に降雨域の裾の部分での降雨強度補正値が小さな値となって、移流ベクトル等の算出に必要なパラメータ推定の精度が低下することを改善するために、降雨強度補正値の前回値がゼロであるものを平均値の算出から除外するようにするのである。
【００４０】
ステップＳ０８４では、Ｒijを前回の補正値で維持するようにして、Ｒij算出処理を終了する。これは、（６）式において、ｎ＝０のときに発生するゼロの除算を防止している。
【００４１】
図２の説明に戻る。ステップＳ８でＲijを算出した後は、ステップＳ９へ進む。
ステップＳ９では、ここまでに説明した一連の補正計算処理を前行程で取得した指定回数だけ繰り返したかどうかを判断し、繰り返しが完了していれば、この一連の処理を終了し、そうでなければステップＳ１へ戻って以降の処理を繰り返す。
【００４２】
以上説明した処理を降雨分布補正工程８は行ない、得られた降雨分布補正データを移流モデルパラメータ演算工程５に渡す。以降の工程は、図１０に基づいて説明した従来法と同様である。
【００４３】
補正演算のＲmin 、β、γ1 、γ2 、指定回数の５つのパラメータを変化させたときの本実施の形態による予測結果と、実際の降雨分布の移動との比較により、最適なパラメータの組を決定する作業を行なっており、これらの値が、Ｒmin ：２．０程度、β：０．２程度、γ1 ：０．８程度、γ2 ：０．３程度、指定回数：５回程度で良好な結果が得られている。但し、本発明はこれらの値に限定されるものではない。
【００４４】
また、（６）式の補正演算は、本実施の形態に限定されず、例えば指数関数等を用いた関数近似により行なうことも可能である。
図５は本発明の第二の実施の形態を示す図である。同図において、図１０と同一の工程には同一の符号を付している。本実施例は、従来技術である移流モデル法の移流ベクトル推定工程中の移流モデルパラメータ演算工程５が修正移流モデルパラメータ演算工程９に置き換わった構成である。図５について説明する。
【００４５】
修正移流モデルパラメータ演算工程９は、図１０における移流モデルパラメータ演算工程５と同様に、時系列降雨分布入力工程１で取得した時系列降雨分布データを用いて移流方程式のパラメータを線形最小二乗推定により算出するのであるが、移流モデルパラメータ演算工程５では、（３）式に示すように、ｘ、ｙ方向に対しての空間差分は降雨分布の移動前のパターンのみを対象としている。それに対して、図５の修正移流モデルパラメータ演算工程９では、次式に示すように、移動後の降雨パターンに対する空間差分を加味しているのが大きな特徴である。
【００４６】
【数３】
0003727762.tif 000004
【００４７】
上式では、添字ｋが移動前の降雨強度を、そして、添字ｋ＋１が移動後、すなわち最新の降雨強度を表している。つまり、ｘ、ｙ方向に対しての空間差分として、移動前の降雨パターンに対する空間差分と移動後の降雨パターンに対する空間差分との平均値を採用する。こうすることで、降雨パターンの急激な変化に対する予測精度を向上することができる。
【００４８】
他の工程は、図１０に基づいて説明した従来法と同様である。
ここで、（７）式にさらに過去の時系列降雨分布データを加味して移流モデルパラメータを求めても良い。また、図５において、時系列降雨分布入力工程１と修正移流モデルパラメータ演算工程９との間に、本発明の第一の実施の形態で説明した降雨分布の補正を行なう工程（図１の補正演算パラメータ入力工程７及び降雨分布補正工程８）を挿入する構成をとることも可能である。
【００４９】
図６は本発明の第三の実施の形態を示す図である。同図において、図１０と同一の工程には同一の符号を付している。本実施例では、基本移流ベクトル演算工程１０、仮想降雨分布演算工程１１、補正移流ベクトル演算工程１２、そして合成移流ベクトル演算工程１３で前述の移流ベクトル推定工程を構成する。図６について説明する。
【００５０】
基本移流ベクトル演算工程１０では、まず、異なる２つの時刻の時系列降雨分布（以降、時刻の早い方を移動前降雨分布、遅い方を移動後降雨分布と呼び、両者の時間間隔をΔｔa とする）において、降雨強度で重み付けをしたときの降雨分布の重心をそれぞれの降雨分布を代表する点とし、この重心の位置の違いをもってΔｔa の間の降雨分布の基本移流ベクトルとする。すなわち、
Ｘ＝移動後降雨分布重心のｘ座標値−移動前降雨分布重心のｘ座標値
Ｙ＝移動後降雨分布重心のｙ座標値−移動前降雨分布重心のｙ座標値
から基本移流ベクトルＶ0 を次式によって求める。
Ｖ0 ＝（Ｘ／Δｔa ，Ｙ／Δｔa ）                                （８）
ここで、次工程での演算の簡単化のために、Ｘ、Ｙの値を観測格子間隔ΔＸ、ΔＹを基準単位とした最も近い整数値に丸め込むと便利である。この場合の基本移流ベクトルＶ0 は次式で表される。
Ｖ0 ＝（Ｘm ・ΔＸ／Δｔa ，Ｙm ・ΔＹ／Δｔa ）                （９）
但し、
Ｘm ＝（Ｘ＋０．５ΔＸ）／ΔＸ      小数点以下切り捨て
Ｙm ＝（Ｙ＋０．５ΔＹ）／ΔＹ      小数点以下切り捨て
とする。
【００５１】
なお、基本移流ベクトルの算出は、本実施の形態に限定するものではなく、例えば、従来技術で説明した雨域追跡法及びそれに類する方法で採用されている移動ベクトルを求める種々の方法が使用できる。
【００５２】
仮想降雨分布演算工程１１では、前行程で得られた基本移流ベクトルに従って移動前降雨分布のデータの移動演算処理を行なう。前工程でＸm 、Ｙm を算出しているのであれば、移動演算は、降雨分布データをｘ方向にＸm 格子分、ｙ方向にＹm 格子分移動させるだけの極めて容易な演算とすることができる。移動演算により得られたこの降雨分布を仮想降雨分布と呼ぶこととする。
【００５３】
補正移流ベクトル演算工程１２では、前行程で得られた仮想降雨分布と移動後降雨分布とに対し、（７）式を適用して移動後の降雨パターンに対する空間差分を加味した移流方程式のパラメータの算出を行ない、これを基にして各格子毎の移流ベクトルを求める。このベクトルを補正移流ベクトルと呼び、ΔＶijと表す。
【００５４】
ΔＶij＝（Δｕi ，Δｖj ）                            （１０）
合成移流ベクトル演算工程１３では、図７に示すように基本移流ベクトルＶ0 と補正移流ベクトルΔＶijとをベクトル加算して合成移流ベクトルＶijを算出する。Ｖ0 を（９）式で表す場合には、（９）式及び（１０）式より、Ｖijを次式で求める。
Ｖij＝（ｕi ，ｖj ）
＝（Ｘm ・ΔＸ／Δｔa ＋Δｕi ，Ｙm ・ΔＹ／Δｔa ＋Δｖj ）（１１）降雨移動予測工程３及び出力工程４は図１０を用いて説明した従来の技術と同様である。
【００５５】
なお、本実施の形態において、補正移流ベクトルの算出時に、過去の時系列降雨分布データを複数用いて移流モデルパラメータを求めても良い。また、図６において、時系列降雨分布入力工程１と基本移流ベクトル演算工程１０との間に、本発明の第一の実施の形態で説明した降雨分布の補正を行なう工程（図１の補正演算パラメータ入力工程７及び降雨分布補正工程８）を挿入する構成をとることも可能である。
【００５６】
図８は本発明の第四の実施の形態を示す図である。同図において、図６と同一の工程には同一の符号を付している。本実施の形態は、第三の実施の形態における基本移流ベクトルを求める別の方法として、過去の雨域内の各格子の移流ベクトルの平均を基本移流ベクトルとするものであり、基本移流ベクトル算出の簡略化によって演算量を減少させ、予測処理の短時間化を可能とするものである。また、本実施例では、基本移流ベクトル演算工程１０、仮想降雨分布演算工程１１、補正移流ベクトル演算工程１２、合成移流ベクトル演算工程１３に加え、パラメータ入力工程１４並びに実績単位時間移流ベクトル記憶工程１５で前述の移流ベクトル推定工程を構成する。以下、図８に従って、本実施の形態を説明する。
【００５７】
パラメータ入力工程１４では、後述する実績単位時間移流ベクトル記憶工程１５及び基本移流ベクトル演算工程１０での処理に必要なパラメータを、キーボードなどの入力手段から取得する。パラメータの内容は後に説明する。
【００５８】
実績単位時間移流ベクトル記憶工程１５では、合成移流ベクトル演算工程１３で算出した各格子の移流ベクトルのｘ、ｙ成分ごとの平均値を算出し、その値を成分とする実績単位時間移流ベクトルＶ1 ＝（ｕ1 ，ｖ1 ）（ｌ＝ｋ＋１、ｋ、…）として逐次記憶する。このベクトル算出の際の平均値の計算量の減少を優先させるために、例えば、パラメータ入力工程１４で取得する格子選択のためのパラメータに従うなどして、計算に使用する格子を限定することも可能である。
【００５９】
基本移流ベクトル演算工程１０では、基本移流ベクトルを求める点においては図７に示すものと同様であるが、ここでは、実績単位時間移流ベクトル記憶工程１５で記憶した、過去の実績単位時間移流ベクトルＶ1 を取り出して、それらのベクトルのｘ、ｙ成分ごとの平均値を算出し、その値を成分とするベクトルを基本移流ベクトルＶ0 とする。実績単位時間移流ベクトルＶ1 をいつの過去にまで遡って平均値算出の対象とするかは、なるべく多くのデータを使用するほうが一般的には望ましいが、計算量や記憶容量との兼ね合いもあるので、ここでは、パラメータ入力工程１４でｍなるパラメータを取得し、最新の値からｍ回前の値までを平均値算出の対象とすることとする。
【００６０】
以上の説明を式で表すと次式のようになる。
【００６１】
【数４】
0003727762.tif 000005
【００６２】
また、ここで第三の実施の形態での説明と同様に、次工程での演算の簡単化のために、格子間隔ΔＸ、ΔＹを基準単位とした整数値への丸め込みを行なっておくと便利である。すなわち、降雨分布の観測時間間隔をΔｔとして、基本移流ベクトルＶ0 を次式により求める。
Ｖ0 ＝（Ｘm ・ΔＸ／Δｔ，Ｙm ・ΔＹ／Δｔ）                  （１３）
但し、
【００６３】
【数５】
0003727762.tif 000006
【００６４】
ここで、Ｘm 、Ｙm 共に小数部は切捨てる。
以上説明したようにして得られた基本移流ベクトルは、後の工程である仮想降雨分布演算工程１１及び合成移流ベクトル演算工程１３で使用される。これらの工程は第三の実施の形態で説明したものと同一である。
【００６５】
なお、本実施の形態においては、演算初期時に実績移流ベクトルの初期値が必要となる。このベクトルの初期値は例えばパラメータ入力工程１４で取得する。この初期値は任意でよいが、そのため、実際の降雨予測に本実施の形態を使用する際には、基本移流ベクトルが安定するまでのしばらくの間、試運転を行なうことが望ましい。
【００６６】
また、図８において、時系列降雨分布入力工程１と仮想降雨分布演算工程１１との間に、本発明の第一の実施の形態で説明した降雨分布の補正を行なう工程（図１の補正演算パラメータ入力工程７及び降雨分布補正工程８）を挿入する構成をとることも可能である。
【００６７】
以上説明した各実施の形態においては移流モデルとして移流方程式を適用したが、大気の乱流拡散を考慮して移流方程式を拡張した移流拡散方程式
【００６８】
【数６】
0003727762.tif 000007
【００６９】
但し、Ｄx 、Ｄy はｘ方向、ｙ方向の拡散係数
を本発明の各実施の形態に適用することも可能である。
また、本発明の各実施の形態は、その各工程をコンピュータに行なわせるためのプログラムを記録したコンピュータ読み取り可能な記録媒体として構成することも可能である。このような記録媒体としては、図１１に示すように、例えば、ＣＤ−ＲＯＭディスクや、フロッピーディスク、更には図示しないがＭＯディスク、ＤＶＤ−ＲＯＭディスクあるいはリムーバブルディスク等の可搬記録媒体２２や、コンピュータ２１からネットワーク回線２５を介して接続されるプログラムサーバ２３の記憶装置２４や、あるいはコンピュータ２１自身が有する記憶装置２６（ＲＡＭやハードディスクなど）が含まれる。更に、プログラムを記録したＲＯＭをコンピュータ２１に取り付けて直接読ませる構成をとることも可能であり、これも記録媒体に含まれる。
【００７０】
【発明の効果】
本発明は、以上詳細に説明したように構成されているので、従来の予測方法では予測精度が低下してしまうような雨域移動の発生している状況下においても、降雨観測領域の各格子毎の移流ベクトルを精度良く算出することが可能となり、その結果、降雨予測の精度向上の効果が得られる。
【図面の簡単な説明】
【図１】本発明の第一の実施の形態を示す図である。
【図２】降雨分布補正工程８で行なう処理のフローチャートである。
【図３】Ｒijの位置関係を示す図である。
【図４】第一の実施の形態におけるＲijの算出方法を示すフローチャートである。
【図５】本発明の第二の実施の形態を示す図である。
【図６】本発明の第三の実施の形態を示す図である。
【図７】移流ベクトルの合成を説明する図である。
【図８】本発明の第四の実施の形態を示す図である。
【図９】従来の雨域追跡法による予測方法を説明する図である。
【図１０】従来の移流モデル法による予測方法を説明する図である。
【図１１】コンピュータ読み取り可能な記録媒体の例を示す図である。
【符号の説明】
１  時系列降雨分布入力工程
２  移動ベクトル演算工程
３  降雨移動予測工程
４  出力工程
５  移流モデルパラメータ演算工程
６  移流ベクトル演算工程
７  補正演算パラメータ入力工程
８  降雨分布補正工程
９  修正移流モデルパラメータ演算工程
１０  基本移流ベクトル演算工程
１１  仮想降雨分布演算工程
１２  補正移流ベクトル演算工程
１３  合成移流ベクトル演算工程
１４  パラメータ入力工程
１５  実績単位時間移流ベクトル記憶工程
２１  コンピュータ
２２  可搬記憶媒体
２３  プログラムサーバ
２４  記憶装置
２５  ネットワーク回線
２６  記憶装置
Ｓ１〜Ｓ９、Ｓ０８１〜Ｓ０８４  処理ステップ
【図１】
0003727762.tif 000008
【図２】
0003727762.tif 000009
【図３】
0003727762.tif 000010
【図４】
0003727762.tif 000011
【図５】
0003727762.tif 000012
【図６】
0003727762.tif 000013
【図７】
0003727762.tif 000014
【図８】
0003727762.tif 000015
【図９】
0003727762.tif 000016
【図１０】
0003727762.tif 000017
【図１１】
0003727762.tif 000018