JP2023086449A - 情報処理装置、情報処理方法及びコンピュータプログラム - Google Patents

Abstract

【課題】雲の状態を高精度に推定する情報処理装置、情報処理方法及びコンピュータプログラムを提供する。【解決手段】本実施形態に係る情報処理装置は、雲をセンシングしたデータに基づき、前記雲の高さ情報を推定する第１推定部と、前記高さ情報に基づき、前記雲を複数の部分に分解する分解部と、前記部分ごとに対象時刻における前記雲の状態を推定する第２推定部と、を備える。【選択図】図１

Description

本実施形態は、情報処理装置、情報処理方法及びコンピュータプログラムに関する。

衛星画像から雲の動きを表す移動ベクトルを計算し、移動ベクトルに基づき、雲の動きを予測（外挿）することによって雲の将来位置を推定する方法がある。この方法では複数の雲が立体的に積み重なっている場合に、雲の移動ベクトルを正確に推定できない。

橋本篤, 小林隆久, 田村英寿, 平口博丸. 太陽光発電出力予測のための衛星画像データを用いた日射量推定・予測モデルの開発. Technical Report N13003, 一般財団法人 電力中央研究所 (CRIEPI), 東京都千代田区大手町 1-6-1, aug 2013. 電力中央研究所研究報告.

本実施形態は、雲の状態を高精度に推定する情報処理装置、情報処理方法及びコンピュータプログラムを提供する。

本実施形態に係る情報処理装置は、雲をセンシングしたデータに基づき、前記雲の高さの情報を推定する第１推定部と、前記高さの情報に基づき、前記雲を複数の部分に分解する分解部と、前記部分ごとに対象時刻における前記雲の状態を推定する第２推定部と、を備える。

第１の実施形態に係る情報処理装置として雲移動推定装置の一例を示すブロック図。 雲構造推定部の詳細ブロック図。 雲分解部により雲を高さのレベルに応じて複数の雲塊に分解する例を示す図。 最近接点間の距離が閾値以下の雲塊同士を同じ雲塊とみなす処理を行う例を示す図。 雲時間変化推定部の一例を示すブロック図。 （A）現在の時間フレームにおける複数の雲塊に対して移流ベクトルを推定した例と、（Ｂ）並進移動される領域を広くして、移流ベクトルを推定した例とを示す図。 補間により雲塊が存在しない位置に対して移流ベクトルを推定した例を示す図。 関連技術において、雲塊がない位置に対して移流ベクトルの推定値が0として算出される例を示す図。 複数の雲塊が重なっている場合における反射強度の計算例を説明する図。 第１の実施形態に係る処理の一例のフローチャート。 時間フレームごとに時間Δτ毎に雲の移流推定が行われる様子を示す図。 第２の実施形態に係る雲分解部の一例のブロック図。 全体として同じ高さのレベルにある雲（雲画像）を複数の雲塊に分解する例を示す図。 隠れた雲の推定の考え方を示す図。 第３の実施形態における雲分解部の一例のブロック図。 マルコフ確率場モデルを説明する図。 第４の実施形態の概要を説明する図。 第４の実施形態に係る雲分解部の一例のブロック図。 第５の実施形態に係る雲時間変化推定部の一例のブロック図。 第６の実施形態に係る情報処理装置として雲移動推定装置の一例のブロック図。 分割パラメータの違いにより雲の分解の相違を示す図。 第７の実施形態の概要を説明する図。 第７の実施形態に係る情報処理装置として雲移動推定装置の一例のブロック図。 衛星画像データベースの一例を示す図。 数値気象計算履歴データベースの一例を示す図。 第８の実施形態に係る情報処理装置として雲移動推定システムの一例のブロック図。 第９の実施形態に係る情報処理装置として日射両推定装置の一例のブロック図。 本発明の実施形態に係る情報処理装置のハードウェアブロック図。

以下、図面を参照しながら本発明の実施形態を説明する。

（第１の実施形態）
図１は、第１の実施形態に係る情報処理装置として雲移動推定装置１００（以下、本装置１００）の一例を示すブロック図である。本装置１００は、気象衛星により検出された気象衛星画像（以下、衛生画像）を用いて、将来の雲の状態（例えば位置）を高精度に予測する技術に関する。衛星画像は雲をセンシングしたデータの一例である。以下、衛星画像について本実施形態の説明に必要な範囲で説明を行う。衛星画像の詳細についは本明細書の末尾に記載した参考文献リストの参考文献１に記載されている。

衛星画像は気象衛星が備えるセンシング装置により雲をセンシングすることにより取得されるデータである。衛星画像は、大きく、可視画像と赤外画像とに分けることができる。可視画像は、一例としてＲＧＢカメラ、モノクロカメラ等の撮像装置により取得される。赤外画像は、一例として赤外線カメラにより取得される。可視画像とは可視光、あるいは近赤外光など、太陽光からくる光の反射光をとらえた画像である。従って、可視画像からは雲の(雲が無い場合は地表面の)反射強度を知ることができる。雲の反射強度は、雲に含まれる水滴(雲粒子)の数、水滴の大きさ、水滴の性質(水か氷かなど)によって影響される。本実施形態では雲が移流する範囲において、雲粒子の大きさ及び性質が大きく変わらないと考え、雲の反射強度は雲粒子の密度と置き換え可能な量であると見なす。

赤外画像は、雲それ自体が持っている熱で発生する輻射光をとらえた画像である。つまり、赤外画像は雲の温度をとらえたものと言ってもよく、赤外画像は雲の温度分布情報を含んでいる。雲の温度のことを輝度温度と呼ぶこともある。雲の温度は、雲の周囲の大気の温度と一致している。そして、大気の鉛直方向の温度プロファイルは、大気の運動を流体力学的、および熱力学的に扱った数値気象モデルを用いて計算できるので、輝度温度は高さ（高度）に変換できる。すなわち、赤外画像は雲の高度を表していると言える。

気象衛星は静止軌道上にあり、一定時間間隔で半球(地球のうち衛星から見える側)、あるいは半球の一部を、複数のバンド（周波数）でそれぞれ撮影している。一つの画像をとる時間幅（コマ）を、以下では時間フレームと呼ぶことにする。本実施形態ではある１つのバンドで撮影した可視画像と、別の１つのバンドで撮影した赤外画像のみを使用する。例えば気象庁のひまわり8号の場合、一つの時間フレームとして10分の間に、全16のバンド（周波数）でそれぞれ半球画像を撮影する。この場合、本実施形態で使用するバンドの画像を可視画像及び赤外画像それぞれについてあらかじめ決めておくものとする。

以下の説明で衛星画像の画素をインデックスの組(i, j)で表す。衛星画像は経度、緯度の方向に沿って画素が並んだものを想定し、i, j はそれぞれ経度、緯度に対応するものとみなす。緯度及び経度の組による位置を表すのに単にｘと表記することもある。ｘは(i, j)と同じ意味であり、互いに置き換え可能なものとみなす。

本実施形態では雲を含む領域（第１領域）から検出されるデータ（センシング情報）として可視画像及び赤外画像を用いるが、雲の位置及び高さを検出可能であれば、可視画像及び赤外画像に限られない。例えばレーザにより雲の位置及び高さを検出可能であれば、レーザ又はデプスカメラなどにより検出されるデータを用いてもよい。

入力部１１０は、衛星から無線を介して、衛星画像を受信する。衛星画像は可視画像と赤外画像とを含む。入力部１１０は、衛星画像を時間フレーム毎に受信する。入力部１１０は、衛星画像を雲構造推定部１０（第１推定部）に提供する。入力部１１０は、受信した衛星画像を記憶部１３０に格納してもよい。

雲構造推定部１０は、衛星画像に基づき、雲構造推定部１０は、雲の水平方向の位置と、鉛直方向における位置（高さ）とを推定する。雲構造推定部１０は、雲の水平方向の位置情報（例えば雲が存在する範囲の情報）と、雲の高さ情報（例えば位置毎の高さ）との少なくとも一方を推定する第１推定部に対応する。

図２は、雲構造推定部１０の詳細ブロック図である。雲構造推定部１０は、雲の水平方向の位置を推定する雲水平位置推定部１１と、雲の立体構造、すなわち水平方向の位置毎の雲の鉛直方向の位置（高さ）を、雲の高さの分布として推定する雲立体構造推定部１２とを備える。

雲水平位置推定部１１は、衛星画像において、雲に対応する画素を特定し、特定した画素に対応する位置を雲の位置として推定する。例えば、可視画像における画素の画素値である反射強度（例えば明るさ及び色の少なくとも一方など）に基づき、当該画素が雲に対応するかを判定できる。その他、赤外画像を用いて雲に対応する画素を特定することも可能である。例えば赤外画像から得られる画素の画素値である輝度温度が地表面温度よりも低ければ、当該画素には雲が映っているものと判定できる。なお、ある画素が雲に対応するかどうかを判定する手法は雲マスクと呼ばれる技術として知られ、上述した方法以外にも様々な方法が提案されている。どのような方法を用いてもかまわない。

雲立体構造推定部１２は、雲水平位置推定部１１により雲に対応すると判定された画素に対して、当該画素に対応する雲の鉛直方向における位置（高さ）を推定する。赤外画像から得られる雲の画素の輝度温度は間接的に雲の高さを示しており、よって、画素の輝度温度に基づき、当該画素に対応する雲の鉛直方向の位置（高さ）を特定できる。輝度温度から高さへの変換は、気温と高さとの関係を対応付けた情報に基づき、輝度温度に対応する高さを、雲の高さとして推定できる。気温と高さとの関係を対応付けた情報として、例えば、数値気象計算の結果として得られる気温の鉛直方向のプロファイルを用いることができる。

雲水平位置推定部１１及び雲立体構造推定部１２で行う処理の結果として、雲構造推定部１０は以下の情報を得る。
H(t, i, j): 時間フレームt,画素(i,j)に対応する位置の雲の高さ
A(t, i, j): 時間フレームt,画素(i,j)に対応する位置の雲からの反射強度
A(t, i, j)≧0である。A(t, i,j)が0又は閾値未満の場合はその画素に対応する位置に雲が無いことを表す。A(t, i,j)が0より大きい又は閾値以上の場合はその画素に対応する位置に雲が有ることを表す。

雲分解部２０は、衛星画像における雲を、雲の各画素に対応する位置の高さに応じて複数の高さのレベルに分解する。同じ緯度・経度であっても高度が違うと大気の流れ（風向・風速）は大きく異なる。従って、異なったレベルの高さにある雲は異なった移動モードにあると考えることができる。同じレベルの高さにある雲画素の集合、あるいは、同じレベルの高さにありかつ互いに隣接する雲画素の集合を雲塊画像と呼び、雲塊画像が示す雲を、雲塊と呼ぶ。

具体的な処理例として、分解したいレベルの数（Lとする）に応じて輝度温度の複数の閾値T₁, T₂,…, T_L (T₁ <T₂ <…<T_L)を設定する。すなわち、複数の温度範囲を設定する。赤外画像から取得される各画素の輝度温度Tに応じて、雲に含まれる画素を次のように分類する。
T < T₁ → レベル1
T₁ ≦T <T₂ → レベル2
…
T_L-1≦T < T_L → レベルL

処理の簡単のため、分解したいレベルの数Lを２としてもよい。この場合、T₁を273K, T₂を地表面の平均的温度としてもよい。T₁を氷の凝固点とすることで、粒子が水の雲と氷の雲とに分離されることを期待できる。この二種類の雲は、でき方及び性質が異なるため、高度的に離れた位置に存在することがしばしばある。

図３は、雲分解部２０により雲を高さのレベルに応じて複数の雲塊に分解する例を示す。この例では、対象となる雲を輝度温度がT₁ 以下のレベル１と、T₁ より大きいレベル２との２つの雲塊（部分）に分解する例が示される。すなわち、対象となる雲を、設定した２つの温度範囲ごとに分解することで、雲を分解した２つの雲塊（部分）を得る。レベル1の雲塊の真下ではレベル2の雲塊が隠れているとみなして、レベル２の雲塊を特定している（隠れている雲を特定する動作については後述する他の実施形態の説明で記載する）。

雲の種類によっては分解によって非常に多くの雲塊（雲塊画像）が生成され、処理が複雑になる可能性がある。雲塊の個数を抑制するため、高さのレベルが同じであり、かつ空間的に近い雲塊は同じ雲塊として扱ってもよい。距離の閾値（dとする）を導入し、最近接点間の距離が閾値d以下の雲塊同士を同じ雲塊と見なす処理を行ってもよい。

図４は、最近接点間の距離が閾値d以下の雲塊同士を同じ雲塊とみなす処理を行う例を示す。図の左には、距離が閾値d以下の雲塊が１つの雲塊Ａとしてまとめられ、図の右には距離が閾値d以下の雲塊が１つの雲塊Ｂとしてまとめられた例を示す。なお、雲塊A及び雲塊Bともに同じ１つの衛星画像に含まれる雲塊である。

処理の結果、雲分解部２０は以下の情報を得る。
N:雲塊の数
n ∈{1, 2, ,,, N}を、雲塊を指すインデックスとして、
H_n(t, i, j): 時間フレームtにおける画素(i,j)に対応する位置の雲塊nの高さ
A_n(t, i, j): 時間フレームtにおける画素(i,j)に対応する位置の雲塊nの反射強度
A_n(t, i, j)≧0である。A_n(t,i,j)=0又は閾値未満の場合は、画素(i,j)に対応する位置に雲塊nが無いことを表す。A_n(t,i,j)が0より大きい又は閾値以上場合は、画素(i,j)に対応する位置に雲塊nが有ることを表す。

なお、各雲塊は画像の一部にのみ存在するため、A_n, H_nはすべての(i, j)に対して値を保持する必要はない。A_n, H_nは、雲塊が存在する位置に応じて、適切にトリミングされていてもよい。

推定処理部３０（第２推定部）は、雲分解部２０が取得した個々の雲（雲塊）の状態の時間的な変化を推定する。推定処理部３０は、個々の雲塊の対象時刻における状態を推定する第２推定部に対応する。推定処理部３０は、個々の雲（雲塊）に対応する複数の雲時間変化推定部３０＿１～３０＿Ｎを備える。雲時間変化推定部３０＿１～３０＿Ｎを特に区別する必要がないときは雲時間変化推定部３０を記載する。各雲時間変化推定部３０は、対応する雲（雲塊）の状態（本実施形態では位置）の時間的な変化を互いに独立に推定する。雲の位置の時間変化は、風により雲が流されること（移流と呼ぶ）に対応する。雲の時間的な変化は、移流以外にも雲の発達・減衰などがあるが、これについては第５の実施形態で記載する。

図５は、雲時間変化推定部３０の一例を示すブロック図である。雲時間変化推定部３０は、移流ベクトル推定部３１と移流計算部３２とを備える。

移流ベクトル推定部３１は、連続した複数の時間フレームの画像から、雲（雲塊）の動く速度(速さと方向)を表すベクトルを推定する。このベクトルを移流ベクトル又は移動ベクトルと呼ぶ。推定の時間間隔（リードタイム）は時間フレームと同じ時間間隔でもよいし、時間フレームよりも短い時間間隔でもよいし、時間フレームよりも長い時間間隔でもよい。また、推定は、対象となる（現在の）時間フレームに対して、当該時間フレームの時刻以降で上記リードタイムごとの時刻（対象時刻）に対して行う。次時間フレーム（次時間フレームの衛星画像はまだ受信されていない）以降も、次時間フレームに対する推定とは独立して継続して行う。すなわち、次時間フレームの衛星画像が受信された後も、現在の時間フレームに対する移流ベクトルの推定を継続してリードタイムごとに行う。次時間フレームに対する推定は、現在の時間フレームとは独立して同様にして行う。

移流ベクトルの推定例として、２つの時間フレームにおける雲（雲塊）の位置間のベクトル、すなわち、１つ前の時間フレームにおける雲（雲塊）の位置と、現在の時間フレームにおける雲（雲塊）の位置間のベクトルを移流ベクトルとしてもよい。雲（雲塊）の位置として例えば雲の重心位置を用いる。

図６（Ａ）は、現在の時間フレームにおける複数の雲塊に対して移流ベクトルを推定した例を示す。

移流ベクトルの他の推定方法として、時間的に前後した画像を並進移動させて、重なる面積が最大化する位置を見つけることで移流ベクトルを推定する技術（並進移動最適化）を用いることもできる(非特許文献１参照)。並進移動最適化については、後述する図６（Ｂ）の説明を参照されたい。

また、画像処理により時間的に連続した画像から物体の移動をとらえるオプティカルフローと呼ばれる技術を適用する方法もある(参考文献２参照)。この場合、例えば、２つ以上の時間フレームに対して雲（雲塊）のオプティカルフローを算出し、オプティカルフローに沿った雲（雲塊）の位置の変化に基づき移流ベクトルを算出する。オプティカルフローでは、計算効率上の工夫や、対象となる物体のスケールが様々に変わりうることへの対応等、技術的な洗練度は高い。根本にある考え方は並進移動の最適化と同じである。

雲（雲塊）が存在しない位置の移流ベクトルを、見えている雲（雲塊）の位置又は移流ベクトルに基づき、空間的な補間により算出してもよい。補間の方法として、カーネルスプライン法を用いてもよい。この方法は、空間的に格子状に配列されていない複数の地点に対して値が与えられている場合に、当該複数の地点以外の地点での値を補間で求める方法である(参考文献３)。各雲塊の位置に対して移流ベクトルが与えられている場合に、雲（雲塊）が存在しない位置の移流ベクトルを、緯度方向、経度方向成分をそれぞれ独立にカーネルスプライン法により補完することができる。このように雲（雲塊）が存在しない位置の移流ベクトルも算出することで、さらに、より高精度な雲（雲塊）の移流推定（位置の時間変化の推定）が可能になる。

図７は、補間により雲塊が存在しない位置に対して移流ベクトルを推定した例を示す。雲塊Ｃと雲塊Ｄとの間の位置Ｅに対して、雲塊Ｃの移流ベクトルと雲塊Ｄの移流ベクトルとの平均が採用されている。

関連技術では、雲の移流推定において、雲（雲塊）の高さが考慮されておらず、二次元的な画像処理が行われていた。すなわち、雲が鉛直方向に立体的に配置され、場合によっては重なり合っていることを考慮していなかった。したがって、高さの異なった雲が重なるように存在し、かつ、それぞれが異なった移流ベクトルを持っている場合には、雲の移流（位置の時間変化）を適正に推定できなかった。本実施形態では、異なるレベルの高さの雲（雲塊）ごとに移流ベクトルを算出するため、後述する移流計算部３２で雲の移流を適正に推定できる。

また、関連技術では、雲が存在しない位置の移流ベクトルを考慮できていなかった。並進移動最適化（図６（Ｂ）参照）において並進移動される領域をある程度細かくすれば、雲塊の移流ベクトルを高精度に推定することができるが、当該領域外の位置に対する移流ベクトルの推定値は0になる。

図８は、関連技術において、雲塊がない位置Ｃに対して移流ベクトルの推定値が0として算出される例を示す。

移流（雲の移動）が風によって引き起こされることを考えると、雲塊が無い位置であっても、風があるのであればその位置の移流ベクトルは0ではなく、正しく推定されなければならない。さもないと、当該位置に、別の雲が移動してきたときに、そこから雲が動かなくなってしまうからである。

上述した並進移動最適化において、並進移動される領域を広くすれば、その領域内はすべて同じ移流ベクトルの値になり、移流ベクトルが0になる問題を低減することは可能である。

図６（Ｂ）は、並進移動される領域を広くして、移流ベクトルを推定した例を示す。設定した領域に属する各雲塊は同じ移流ベクトルを有すると見なす。結果的に雲塊が無い位置（空間）の移流ベクトルも推定したことになる。

しかしながら、並進移動される領域を広くとると、領域内部の個々の雲塊等の移流ベクトルの違い（図６（Ａ）参照）を見逃すことになり、雲塊の移動の推定精度を落とすことになる。この問題は、空間的な雲の位置と、雲間の間隙とを区別して取り扱えていないことにより発生している。

本実施形態では、並進移動最適化を用いて移流ベクトルを算出した場合であっても、雲が存在しない位置の移流ベクトルも補間により算出するため（図７参照）、雲の移流を高精度に推定できる。

移流計算部３２は、移流ベクトル推定部３１で推定された移流ベクトルを用いて、雲（雲塊）の状態の時間変化として、雲の位置の時間変化（移流）を計算する。簡単な方法として、移流ベクトルに従って、雲（雲塊）が移動し続ける、すなわち同じ方向に雲（雲塊）が移動し続けるとして、雲（雲塊）の位置の時間変化を計算する。流体中にある粒子の密度分布が流体の流れによって変化する様子を計算する方法（流体力学の分野では伝統的に扱われてきた移流計算の方法）を衛星画像に適用して、雲粒子の密度の時間発展を計算することで雲（雲塊）の位置の時間変化を計算してもよい（参考文献６参照）。この場合、上述した雲が存在しない位置の移流ベクトルも用いることで高精度な計算が可能となる。また移流後に雲（雲塊）の雲粒子の密度が変わることで反射強度も変わるため、変更後の反射強度も推定できる。

処理の結果として、雲時間変化推定部３０（移流ベクトル推定部３１、移流計算部３２）は雲塊nに対して、τを推定のリードタイムとして、以下の情報を得る。

H’_n(t,τ i, j): 時間フレームtにおけるリードタイムτ(時刻t+τ)の推定における画素(i,j)に対応する位置の雲塊nの高さ（移流ベクトルの推定において、高さは変わらないと仮定して、水平方向成分のみを算出してもよく、この場合、高さの値は一定としてよい）
A’_n (t, τ,i, j): 時間フレームtにおけるリードタイムτ(時刻t+τ)の推定における画素(i,j)に対応する位置の雲塊nの反射強度（移流の前後で反射強度は変わらないと仮定してもよい。この場合、反射強度の値は一定である）
A’_n(t, τ, i, j)≧0である。A’_n(t,τ,i,j)=0又は閾値未満の場合は、画素(i,j)に対応する位置に雲塊nの雲が無いことを表す。A’_n(t,τ,i,j)=0より大きい又は閾値より大きい場合は、画素(i,j)に対応する位置に雲塊nの雲が有ることを表す。

雲合成部４０は、雲時間変化推定部３０が推定した各雲（雲塊）の位置及び高さに関する情報の集合を、移動後の雲の雲情報として取得する。例えば、移動後の雲の雲情報は、雲時間変化推定部３０が得た各雲塊nに対する情報（H’_n(t,τ i, j)、A’_n (t, τ,i, j)）を含む。雲情報は現在の時間フレームにおけるリードタイムごとの時刻（対象時刻）に対して得られる。雲情報は、対象時刻における雲の立体的な構造を表しているといえる。

雲合成部４０は、計算された雲の立体的な構造から可視画像を再構成してもよい。すなわち、計算された雲の立体的な構造が存在すると仮定した場合に、当該構造に対して衛星によって取得されるであろう可視画像を推定又は計算により取得してもよい。このためには、各位置の反射強度の情報を含む画像を再構成画像として生成すればよい。具体的には、対象とする位置に雲塊が存在しない場合はその位置の雲反射強度は0とし、雲塊が一つしかない場合はその雲塊の反射強度をその位置の反射強度（A’_n (t, τ,i, j)）とする。対象とする位置に複数の雲塊が重なっている場合には、以下のように反射強度を計算することができる。

図９は、複数の雲塊が重なっている場合における反射強度の計算例を説明する図である。上層、下層に反射強度がそれぞれA₁,A₂の雲（雲塊）があるとする。このとき、下層の雲で反射された上向きの光が再度、上層の雲で反射される場合がある。このため、反射の過程はやや複雑になる。二つの雲の反射は次の式で与えられる反射強度

と等価になる。

三個以上の雲塊が重なった場合は、式(1)を繰り返し適用することで、反射強度を計算できる。

処理の結果として、再構成画像として、以下の情報が得られる。
A’(t, τ, i, j): 時間フレームtにおけるリードタイムτ(時刻t+τ)の推定における画素(i,j)に相当する位置の反射強度

図１０は、第１の実施形態に係る処理の一例のフローチャートである。

最初に初期の時間フレームt_initを時間フレームtとして設定する（Ｓ１０１）。この時点で、時間フレームが一つ前の画像はすでに取得ずみであるもとする。

入力部１１０が時間フレームtの衛星画像（雲を含む領域のセンシング情報）を受信する（Ｓ１０２）。

雲構造推定部１０が、衛星画像に基づき、雲の存在する位置及び位置毎の高さ（雲の3次元的位置）を推定する（Ｓ１０３）。

雲分解部２０が、位置毎の雲の高さ（雲の高さの分布）に応じて、雲を複数の雲塊（部分）に分解する（Ｓ１０４）。

リードタイムτを０に設定する（Ｓ１０５）。

雲時間変化推定部３０が、個々の雲塊の状態の時間変化（本例では位置の時間変化）、すなわちリードタイムτ後の位置を計算する（Ｓ１０６）。

雲合成部４０が、推定された位置（３次元位置）の雲塊の情報の集合を出力情報として得る（Ｓ１０７）。また雲合成部４０は、推定された位置の雲塊を再合成した画像（可視化画像）を生成する（同Ｓ１０７）。

τにΔτを加算し（Ｓ１０８）、τが最大値τ_maxになるまで（Ｓ１０９）、ステップＳ１０６～Ｓ１０８を繰り返す。Δτは移流推定を行う時間間隔を表す。

tにΔtを加算し（Ｓ１１１）、tが終了時刻（tendとする）に達するまで（Ｓ１１０）、ステップＳ１０１～Ｓ１０９を繰り返す。Δtは衛星画像の時間フレーム間隔である。例えば、ひまわり8号の半球画像の場合、10分である。

移流推定の時間間隔Δτは、時間フレーム間隔Δtとは独立に決めてよい。ただし、Δτ=Δt (あるいは Δτ=Δt /m :m整数)とすると、移流推定した時刻と同じ時刻に対する衛星画像が得られるため、移流推定の精度を向上させやすいなどの利点がある。

図１１は、時間フレームごとに衛星画像（可視化画像及び赤外画像）が入力され、時間フレームごとに時間Δτ毎に雲の移流推定が行われる様子を示す。ここではΔτ=Δtとしている。時間フレームt-Δtで行われたリードタイムΔτ後の移流推定結果が、時間フレームtのリードタイム0における雲の状態に対応するものであることが、縦の点線で示されている。

以上、本実施形態によれば、衛星画像に基づき雲を高さのレベルが異なる複数の部分に分解し、部分ごとに雲の位置の時間変化を推定することで、雲の位置の計算精度を向上させることができる。

（第２の実施形態）
同じ性質を共有する一つの雲の場合は、雲の中の個々の部分で高さのレベル（高度）のばらつきがあっても、これらの部分を分解せずに、同じ雲（雲塊）として扱うことが好ましい。あるいは逆に、同じ輝度温度でも、性質の違う複数の部分を含む雲の場合は、当該部分を別々の雲（雲塊）として扱うことが好ましい。第２の実施形態では、異なる性質を有する複数の部分を雲が含む場合は、これらの部分が同じ高さのレベルの雲であっても複数に分解することを実現する。異なった性質を持つ雲は、画像としてのテキスチャー（文様）の違いとなって現れることがあるため、本実施形態では、テキスチャーの違いを用いて、雲を分解する。

図１２は、第２の実施形態に係る雲分解部２０の一例のブロック図である。その他の要素は第１の実施形態と同様である。雲分解部２０は、画像特徴量計算部２１と、画像セグメンテーション部２２とを備える。

画像特徴量計算部２１は、雲を含む領域の撮像画像である衛星画像（赤外画像）の画素ごとに、画素の輝度温度と周辺画素の輝度温度とに基づき１つ以上の特徴量を計算する。特徴量としてはＧＬＣＭ（参考文献４参照）などが広く知られており、本実施形態でもＧＬＣＭを用いることができる。

画像セグメンテーション部２２は、特徴量の空間に各画素の特徴量をマッピングし、特徴量を表す座標（点）間の距離に基づいて画素の類似性を判断し、画素のクラスタリングを行う。生成された各クラスタに対応する画素の集合が雲塊の画像に対応する。クラスタリングの例として、距離の近い画素同士から順番に同じクラスタに統合する方法がある（参考文献５参照）。

図１３は、全体として同じ高さのレベルにある雲（雲画像）を複数の雲塊に分解する例を示す。クラスタリングにより領域（クラスタ）Ｆと領域（クラスタ）Ｇとに分割される。クラスタＦの明暗の変化パターン（テキスチャー）は均質であり、クラスタＧにおける明暗の変化パターン（テキスチャー）も均質であるが、領域Ｆ及び領域Ｇ間でテキスチャーの内容は異なっている。よってクラスタＦが１つの雲塊を表し、クラスタＧが別の１つの雲塊を表す。このようにクラスタＦとクラスタＧは近接あるいは重なっており、同じ高さのレベルを有するものの、テキスチャーの違いを利用して、別個の雲塊として扱うことができる。

以上、本実施形態によれば、同じ高さのレベルにある雲塊を複数に分解することでより精度の高い移流推定が可能となる。

（第３の実施形態）
上層にある雲に隠れているため画像では見えない下層の雲を推定する新たな機能を雲分解部２０に追加する。

図１４は、本実施形態における隠れた雲を推定する考え方を示す。図１４（Ａ）は、高いレベル（高層）の雲Ｊ１と、低いレベル（低層）の雲Ｋ１、Ｋ２とを示す。雲Ｊ１は、雲Ｋ１，Ｋ２と分離している。すなわち、雲Ｊ１の周りにほかの雲が隣接していない。この場合、高層の雲Ｊ１の下には他の雲は無いと推定するのがもっともらしい。一方、図１４（Ｂ）は、高いレベル（高層）の雲Ｊ２と、低いレベル（低層）の雲Ｋ３とを示す。高層の雲Ｊ２に低層の雲Ｋ３が隣接している。この場合、高層の雲Ｊ２の下にも低層の雲Ｋ３が存在すると推定するのがもっともらしい。すなわち、隠れた部分の雲を推定する基本の考え方は、雲がどのような形に広がるかの統計的なモデルをベースに、隠れた部分に雲があると考えるのがもっともらしいか、あるいは無いと考えるのがもっともらしいかを決めることである。

図１５は、本実施形態における雲分解部２０の一例のブロック図である。雲分解部２０は、雲広がり統計モデル部２３と、雲広がり推定部２４とを含む。

雲広がり統計モデル部２３は、雲広がりの統計モデルを設定する。統計モデルは記憶部１３０に格納されており、雲広がり統計モデル部２３は、記憶部１３０から統計モデルを読み出す。ここでは統計モデルとして、確率モデルであるマルコフ確率場モデルを設定する(参考文献７参照)。マルコフ確率場モデルは、画像解析で良く利用されるモデルである。

図１６は、マルコフ確率場モデルを説明する図である。雲が有る画素を画素ＰＲ、雲が無い画素を画素ＰＧ、雲の有無を判定したい画素（対象画素）を画素ＰＴとする。マルコフ確率場モデルでは、ある画素に雲が存在する確率は、隣接する画素の雲の有無によって決まる。例えば、図１６（Ａ）のように、対象画素ＰＴの周りに、雲の有る画素ＰＲが多ければ、対象画素ＰＴにも雲が有る確率は高いと考えられる。逆に、図１６（Ｂ）のように、対象画素ＰＴの周りに雲がある画素ＰＲが少なければ対象画素ＰＴに雲のある確率は小さい。このように、対象画素ＰＴにおける雲存の在確率を、隣接する画素における雲の有無に応じて決まるように設定するモデルがマルコフ確率場モデルである。

高層の雲（雲塊）の下に隠れた部分の領域を、当該領域に下層の雲（他の雲又は他の雲塊）があるかどうか推定する対象領域に設定する。対象領域に含まれる各画素を対象画素ＰＴとして、整数の組(i, j)で表す。対象領域に含まれる画素(i, j)に雲があるかどうかを表す変数を以下で表す。変数は、確率モデルの変数であり、確率的に推定される対象となる。

対象領域の周囲との境界領域では、下層雲があるかが分かっている。この境界領域に含まれる画素は整数kで表す。境界領域に含まれる画素に雲があるかどうかを、以下で表す。

bkは、衛星画像から一意に確定される。例として図１６（Ｃ）には、図１４（Ｂ）の高い雲の下にある領域における低い雲の有無を推定する対象画素ＰＴが示される。
境界領域において、低い雲が無い画素ＰＧ（b_k=0）、低い雲がある画素ＰＲ（b_k=1）が示される。

マルコフ確率場モデルによると、対象領域における雲の有無は、以下の式(2)で示す変数のセットで与えられる。

当該変数のセットに基づき、確率モデル（統計モデル）は、例えば以下の式で与えられる。

γ(≧0)は、隣接する画素に雲が存在する確率の強さを表す統計モデルのパラメータである。過去のデータに当てはめるなどして予め決めておく。

雲広がり推定部２４は、複数の雲塊（部分）のうちの第１雲塊（第１部分）の境界において第１部分に隣接する他の雲塊（他の部分）と隣接する範囲に基づいて、第１雲塊（第１部分）の下層において他の雲が存在するか否か又は存在する他の雲の範囲を推定する。

より詳細には、雲広がり推定部２４は、雲広がり統計モデル部２３により設定される確率モデル（式(3)参照）に基づいて、対象領域における各画素ＰＴに雲があるかどうか推定を行う。確率モデルに含まれる変数の値を推定する手法は任意でよいが。簡便な手法として、MAP(Maximum a posteriori)推定がある。MAP推定の式を以下に示す。

MAP推定を行う場合、式(3)を最大化又は準最大化するように、確率モデルに含まれる変数の値を推定する。
雲広がり推定部２４は、対象領域（上層の雲に隠された領域）に存在すると推定された下層の雲を、図１４（Ｂ）における低い雲（雲Ｋ３）に統合する。

以上、本実施形態によれば、高層の雲に隠れた他の雲を検出することができる。

（第４の実施形態）
第３の実施形態と異なる手法により高層の雲に隠れた領域に低層の雲が存在するかどうかを推定する。当該推定する機能を雲分解部２０に追加する。

図１７は、第４の実施形態の概要を説明する図である。下層の雲が移動して下層の雲の一部又は全部が上層の雲に隠れた場合、あるいは上層の雲が移動して下層の雲の一部又は全部を隠した場合、前の時刻の画像から隠れた雲の有無を推定できる。時間フレーム(t-Δt)において離れた位置に存在した低い雲Ｍの一部ＯＡが、現在の時間フレームtで高い雲Ｎの下に入り込んでいる。第１～第３の実施形態では直前の時間フレーム(t-Δt)の衛星画像から現在の時間フレームtにおける雲の位置を移流計算に基づいて推定するが、第４の実施形態では、この移流計算を利用して、上層の雲に隠された雲の有無を推定する。

図１８は、本実施形態に係る雲分解部２０の一例のブロック図である。雲分解部２０は不確定部分特定部２６と移流推定合成部２７とを備えている。

不確定部分特定部２６は、上層の雲に隠されて、下層の雲の有無が不明な画素（対象画素）を特定する。上層の雲に隠れている領域の画素が対象画素に該当するため、衛星画像において雲がある領域を不確定領域として指定する。不確定領域に含まれる画素を対象画素とし、対象画素における雲の高度（輝度温度）を取得する。

移流推定合成部２７は、前の時間フレームt-Δtの衛星画像に基づき推定された雲塊（雲の部分）ごとの移流計算の結果に基づき、現在の時間フレームtの衛星画像に基づき得られた複数の雲塊の下層において他の雲が存在するか否か（又は他の雲が存在する範囲）を推定する。t-Δtは一例として第１時刻、ｔは一例として、第１時刻より後の第２時刻に対応する。

より詳細には、移流推定合成部２７は、前の時間フレームt-Δtで得られた移流計算の結果の中から、現在の時間フレームtに対応する雲の情報（３次元位置又は立体構造）を取得する。移流推定合成部２７は、取得した情報に基づき、不確定領域に含まれる各画素（対象画素）の位置に、現在の時間フレームの衛星画像において見えている雲よりも低い高度（高さレベルが低い）の雲が存在するかを判断する。移流推定合成部２７は、低い高度の雲が存在すると判断した場合は、現在の時間フレームtにおいて当該対象画素の位置に、高さレベルの低い高度の雲が存在するとみなす。移流推定合成部２７は、当該存在すると見なされた雲を、時間フレームtの衛星画像において、高さレベルが高い高度の雲に隣接して見えている高さレベルの低い雲に統合する。

以上、本実施形態によれば、高層の雲に隠れた他の雲を検出することができる。

（第５の実施形態）
上昇気流などの理由により、積乱雲などの雲が発生し、雲が急速に発達することがある。一方、水蒸気が過飽和な領域で作られた雲は、水蒸気の飽和度が小さい場所に移動すると、減衰する。このような雲の変化は画像上では雲が消滅しているものととらえることができる。第１の実施形態では、雲の状態の時間変化として位置の時間変化を推定したが、第５の実施形態では、位置の時間変化に加えて、雲の発達及び衰退を推定する。

図１９は、本実施形態に係る雲時間変化推定部３０の一例のブロック図である。図１９の雲時間変化推定部３０は、第１の実施形態における雲時間変化推定部３０に、雲マッチング部３３と、発達・減衰係数推定部３４と、雲発達・減衰計算部３５とが追加されている。雲時間変化推定部３０は、一つ前の時間フレームt-Δtの衛星画像に基づく雲塊（雲の部分）ごとの発達状態と、現在の時間フレームtの衛星画像に基づく雲塊ごとの発達状態との差分に基づいて、現在の時間フレームtの衛星画像に基づく雲塊ごとの雲の発達状態の時間変化を推定する。つまり、対象時刻における発達状態を推定する。雲の発達状態の時間変化は、雲粒子の密度の増加あるいは減少を意味する。より詳細には以下の通りである。

雲マッチング部３３は、現在の時間フレームtで得られた雲塊と、一つ前の時間フレームt-Δtで得られた雲塊とを対応づける。対応付けは、雲塊の位置に基づいて行われる。雲塊の位置として、雲塊が占める領域の中心位置（重心位置）と当該中心位置に対応する高さとを用いることができる。雲マッチング部３３は、雲塊の位置間の距離を計算し、距離が最も近いもの同士を対応付ける。以下の説明において、現在の時間フレームの雲塊をインデックスnで表し、当該雲塊に対応する一つ前の時間フレームの雲塊も同じインデックスnで表す。

発達・減衰係数推定部３４は、現在の時間フレームにおける雲塊nの画素毎の反射強度の合計と、時間フレームが一つ前の雲塊nの画素毎の反射強度の合計との関係に基づいて、発達・減衰係数αとして計算する。一例として、これら２つの合計の比に基づいて発達・減衰係数αを計算する。２つの反射強度の比によって、発達・減衰係数αを計算する計算式の例を、以下の式に示す。反射強度は雲粒子の密度を表すことから、αは時間Δtの間に雲粒子の密度が増加あるいは減少する割合を表している。

雲発達・減衰計算部３５は、発達・減衰係数αに基づいて、移流推定した雲の反射係数を、下記の式(6)により修正する。修正後の反射係数を

とする。式(6)により移流推定した雲の反射係数を修正することで、雲の発達状態の時間変化（雲粒子の密度の増加あるいは減少）を推定することができる。

以上、本実施形態によれば、雲の状態の時間変化として位置の時間変化に加えて、雲の発達及び衰退を推定することができる。

（第６の実施形態）
本実施形態は第1の実施形態に対して、雲の分解とそれに伴う移流ベクトルの推定とにおける不確定性をアンサンブル予測により取り扱うものである。第1の実施形態で取り上げた輝度温度を用いた雲の分解において、閾値の設定に不確定性があることを例として説明する。

すなわち、第１の実施形態では、与えられた輝度温度の閾値T₁に基づいて、閾値T₁より低い輝度温度を有する高層にある雲と、閾値T₁以上の輝度温度を有する下層の雲とに雲（雲画像）を分解することを示した。T₁を摂氏0℃(273K)に設定することで、水と氷の雲に分離できる。しかし、実際には雲は必ずしも摂氏0℃以下なら氷というわけではない。また、衛星画像上で観測される輝度温度が実際の値と正確には一致しない可能性もある。そこで、本実施形態では、輝度温度の閾値T₁の値が取りうる候補を生成し、候補ごとに雲分解部２０、雲時間変化推定部３０及び雲合成部４０の処理を行う。候補ごとに行った処理の結果を合成する。これにより不確実性を補完し、移流推定の精度を向上させることができる。以下、本実施形態について詳細に説明する。

図２０は、本実施形態に係る情報処理装置として雲移動推定装置１００Ａの一例のブロック図である。

雲分割パラメータ推定部５０は、輝度温度の閾値T₁の値が取りうる候補
T₁ ⁽¹⁾, T₁ ⁽²⁾, …, T₁ ^(m),
を生成する。例えば、輝度温度の閾値が273Kの近くで複数な値を取り得る可能性を考慮して
T₁ ⁽¹⁾= 273-(m/2)ΔT, T₁ ⁽²⁾ = 273-(m/2-1)ΔT, …, T₁ ^(m) = 273+(m/2-1)ΔT
を候補として生成する。ただし、ΔTは任意の値であり、例えばΔT=1である。

図２１は、分割パラメータT₁の違いにより雲の分解の相違を示す。T₁が下がるほど下層の雲の占める割合が高くなる。また、それにより雲の移流ベクトルも変わる可能性がある。

T₁の各候補に対して、雲分解部２０、雲時間変化推定部３０及び雲合成部４０の処理がそれぞれ独立にm個行われる。

推定合成部６０は、独立に行われた雲位置の推定結果を合成する。合成の方法として以下に２つの例を示すが、他の方法も可能である。

（平均値の導出）
推定したい量をXとする。Xは例えば雲の中心位置、又は、ある地点の日射量などである。この時、平均値は、

で計算される。ただし、X⁽ⁱ⁾は、i番目の雲合成部４０から推定したXの値である。

（確率的な推定）
例えばXの確率分布P(X)を、カーネル密度推定を用いて以下の式で算出する。

K(X, X⁽ⁱ⁾)はカーネル関数と呼ばれるもので、例えば以下の式で与えられる。

ここでσはカーネル関数の幅を表すパラメータである。

カーネル密度推定を用いた場合、例えばカーネル関数を累積確率分布に変換し、所定の確率（例えば0.8）に対応するＸを最終的な推定結果としてもよい。また、０～１の範囲で乱数を生成し、乱数の値に対応するXを最終的な推定結果としてもよい。乱数を複数回生成して複数のXを算出し、複数のXの平均を最終的な推定結果としてもよい。その他の方法でカーンル関数に基づき、最終的な推定結果を算出してもよい。

以上、本実施形態によれば、閾値となる複数の輝度候補ごとに移流推定を行い、温度候補ごとの移流推定を統合することにより、単一推定で発生しやすい大きな誤差を回避し、移流推定の精度を向上させることができる。

（第７の実施形態）
本実施形態は、第1の実施形態に対して数値気象計算の結果を利用する機能を追加したものである。

図２２は、本実施形態の概要を説明する図である。雲時間変化推定部３０で推定される移流ベクトルＶ１は、推定時点（τ＝０）では雲の動く方向を高精度に表していると考えられる。第１の実施形態では基本的にその後も移流ベクトルが変わらないものと仮定して雲の将来（τ＝１，２，・・・）の位置を計算する例を記載した。しかしながら、風が時間とともに変化する場合、雲の移動方向も、リードタイムτの増加とともに最初に推定した移流ベクトルから外れると予想される。本実施形態では数値気象計算の結果を利用して、移流ベクトルの誤差を修正する。

数値気象計算では大気のダイナミクスの基づく気象の変化を計算する。したがって、予測リードタイムの増加とともに、数値気象計算により得られる風から導かれる移流ベクトルの方が、移流ベクトルＶ１より正確な推定値となる可能性が高い。移流ベクトルＶ２＿０、Ｖ２＿１、・・・、Ｖ２＿６は、τ＝０，２，・・・、６において数値気象計算により得られる風から導かれる移流ベクトルであり、時間の経過に応じて、移流ベクトルＶ１より信頼性が高くなる。

そこで、移流ベクトルＶ１と移流ベクトルＶ２（Ｖ２＿０～Ｖ２＿６）とにリードタイムに応じた重みを設定し、両ベクトルを合成することで、τ＝０，２，・・・、６において、移流計算に用いる移流ベクトルＶ３＿０～Ｖ３＿６を算出する。重みの設定では、移流計算の開始時（τ＝０）は移流ベクトルＶ１の重みの方を大きくする。例えば移流ベクトルＶ１の重みを１、移流ベクトルＶ２＿０の重みを０とする。リードタイムの増加に応じて徐々に移流ベクトルＶ１の重みを小さくし、移流ベクトルＶ２の重みの方を大きくしていく。

図２３は、本実施形態に係る情報処理装置として雲移動推定装置１００Ｂの一例のブロック図である。

数値気象計算部２００は、大気（風等）の物理的なダイナミクスを数値的にモデル化し、その時間的な発展を計算する。数値気象計算部２００は、モデルを実装したシミュレーションソフトを用いて計算してもよい。任意地点及び任意時刻の気温、風向・風速、気圧、湿度の推定値を計算することができる。数値気象計算部２００により計算される時刻t、位置x、高度（高さ）zにおける風をW(t,x,z)と表す。W(t,x,z)は少なくとも風向を含み、その他、気温、風速等を含んでもよい。なお、本装置１００Ｂは、数値気象計算の結果を、機関(気象庁など)のサーバから受信しても良い。

風移流ベクトル算出部３６は、数値気象計算の結果における風の場（例えば、ある時刻で、どの位置・高さで、どの方向に、どのような風速の風が吹いているか）から、風移流ベクトルV_s (t,x)を算出する。風移流ベクトルは一例として風向及び風速の情報を含む。風が雲の移流をもたらすため、雲塊nの高さH_nを使って、風移流ベクトルを以下の式で表す。

以下では、このように数値気象計算の結果から計算した移流ベクトルを風移流ベクトルと呼ぶ。

移流ベクトル合成部３７は、移流ベクトル推定部３１により推定された時間フレームtにおける移流ベクトルV(t,x)と、風移流ベクトル算出部３６により算出された風移流ベクトルV_s(t+τ,x)とを合成して、リードタイムτにおける合成移流ベクトル（移流ベクトルV'(t,τ,x)）を算出する。例えば、リードタイムτにおけるV(t,x)とV_s(t+τ,x)との重みを、それぞれ第１重みに対応するw_τ(x)、第２重みに対応する1-w_τ(x)として、以下の式により、合成移流ベクトルを計算する。重みw_τ(x)、1-w_τ(x)は、後述する重み計算部２８０により計算される。

衛星画像データベース２２０は、ある一定期間（例えば数か月、数年）の過去の衛星画像データを記憶したデータベースである。
図２４は、衛星画像データベース２２０の一例を示す。時間フレーム毎に可視バンドのデータ（可視画像）及び赤外バンドのデータ（赤外画像）が格納されている。可視バンドのデータは画素毎に反射強度の値を含む可視画像である。赤外バンドのデータは、画素毎に輝度温度の値を含む赤外画像である。

データベース解析部２４０は、第１の実施形態における同一名称の要素と同様の機能を有する雲構造推定部２４１、雲分解部２４２及び移流ベクトル推定部２４３を備える。データベース解析部２４０は、衛星画像データベース２２０に記憶されている過去の衛星画像データの各時間フレームに対して、第１の実施形態と同様の処理により、移流ベクトルの推定を行う。

移流ベクトル誤差履歴算出部２５０、過去の衛星画像データに対して算出された移流ベクトルの誤差を、下記の式に基づいて計算する。

V(h,x)は、衛星画像データベース２２０における過去の時刻hのデータから計算された移流ベクトルを表す。V(h,x)は、時刻hにおける移流ベクトルに対応する。V(h,x)は、第３時刻の衛星画像に基づき推定した、位置ｘにおける移動の方向及び速さを表すベクトルに対応する。
V(h+τ,x)は、衛星画像データベース２２０における過去の時刻h+τのデータから計算された移流ベクトルを表す。時刻h+τにおける移流ベクトルを表す。V(h+τ,x)は、第３時刻より後の第４時刻の衛星画像に基づき推定した、位置ｘにおける移動の方向及び速さを表すベクトルに対応する。
一例として時間フレーム間隔Δt＝Δτであるが、これに限定されない。
e(h,τ,x)は、上述の２つの移流ベクトルの誤差（差分）である。

数値気象計算履歴データベース２３０は、数値気象計算の結果を記憶したデータベースである。

図２５は、数値気象計算履歴データベース２３０の一例を示す。衛星画像データベース２２０でカバーされている期間及び領域と共通する期間及び領域の数値気象計算結果が格納されている。時刻ごとに気象変数が格納されている。図の例では、反射強度、気温、風などの気象変数が格納されている。

時刻毎の各気象変数は、表形式を有している。気象変数が反射強度の場合、反射強度の情報は、二次元の位置(緯度・経度)毎に反射強度に関する値を格納した表である。気象変数が気温の場合、気温の情報は、三次元の位置(緯度・経度・高度)毎に気温に関する値を格納した表である。気象変数が風の場合、風の情報は、三次元の位置(緯度・経度・高度)毎に、風場に関する値を格納した表である。

図２５の数値気象計算履歴データベース２３０の時刻は必ずしも衛星画像の時間フレームと一致する必要は無い。この場合、補間により時間フレームに対応した時刻における気象変数の値を得ればよい。

風移流ベクトル履歴算出部２６０は、数値気象計算履歴データベース２３０に基づき、データベース解析部２４０の雲分解部２０により取得された各雲塊に対して、風移流ベクトルを算出し、風移流ベクトルの履歴を算出する。風移流ベクトルの算出方法は、風移流ベクトル算出部３６と同様である。

風移流ベクトル誤差履歴算出部２７０は、風移流ベクトル履歴算出部２６０により算出された風移流ベクトルの履歴と、データベース解析部２４０により得られた移流ベクトルの履歴とに基づいて、風移流ベクトルと移流ベクトルとの誤差の履歴を得る。風移流ベクトルの誤差は、式(12)に基づいて算出する。

V_s(h+τ,x)は、リードタイムτにおける風移流ベクトルである。V_s(h+τ,x)は、第４時刻に対する数値気象計算の風向及び風速を表すベクトルに対応する。
e_s(h+τ,x)は、時刻h+τにおける風移流ベクトルV_s(h+τ,x)と、時刻h+τにおける移流ベクトルV(h+τ,x)との誤差（差分）である。

重み計算部２８０は、移流ベクトルV(t,x)に与える重みw_τ(x)と、風移流ベクトルV_s(t+τ,x)に与える重み1-w_τ(x)とを計算する。

重み計算部２８０は、重み

を、以下の式で計算する。したがって、重みは、過去の履歴に基づいて計算される。

重みは位置（地点）x毎に計算される。地点ごとに重みが大きく変化することは無いと仮定される場合、以下の式のように、全地点で平均した重みを用いてもよい。

重みの計算のために用いる過去の履歴の時刻（第３時刻及び第４時刻）は、現在の時間フレームより前の時刻であれば、どの時刻でもかまわないが、現在の時間フレームに近い場合は、より高精度な移流推定が期待できる。

以上、本実施形態によれば、数値気象計算の結果を用いることで、雲の移流推定（特に雲の位置の推定）を高精度に行うことができる。

（第８の実施形態）
第８の実施形態は第1実施形態における雲位置の推定の結果と、数値気象計算の結果とを合成して、推定された雲（雲塊）の反射強度を推定する。第７の実施形態では、雲位置の推定の途中で得られる移流ベクトルと、数値気象計算の結果に基づく風移流ベクトルとを合成した。第８実施形態では移流推定の結果として得られる雲（雲塊）の反射強度と、数値気象計算の結果に基づく反射強度との合成を行う。短時間予測に優れている第１の実施形態における雲位置の推定と、長時間予測に優れている数値気象計算の結果とを合成することで、予測時間の長さに拘わらず、高精度な予測が可能になる。

図２６は、本実施形態に係る情報処理装置として雲移動推定システムの一例のブロック図である。

雲移動推定装置１００は、第１の実施形態に係る雲移動推定装置（図１参照）と同じである。数値気象計算部２００、衛星画像データベース２２０及び数値気象計算履歴データベース２３０は、第７の実施形態（図２３参照）と同じである。

過去履歴用雲移動推定装置３１０は、衛星画像データベース２２０における衛星画像の過去履歴に基づき、過去の時間フレームの時刻h、リードタイムτ、位置xの反射強度の推定値A'(h,τ,x)を算出する。過去履歴用雲移動推定装置３１０は、雲移動推定装置１００と同じ構成を有する。雲移動推定装置１００が、過去履歴用雲移動推定装置３１０を兼ねてもよい。この場合、過去履歴用雲移動推定装置３１０を省略してよい。

反射強度推定誤差履歴算出部３２０は、衛星画像データベース２２０に基づき、過去履歴用雲移動推定装置３１０により推定された反射強度A'(h,τ,x)の誤差を、以下の式により計算する。

反射強度A'(h,τ,x)は、推定対象となる時刻（h+τ）における雲を含む領域（衛星画像の撮像領域）の各位置の反射強度を表す。
A(h+τ,x)は、衛星画像データベース２２０に記憶されている時刻h+τにおける位置xの反射強度である。つまりA(h+τ,x)は、時刻h+τにおける衛星画像（撮像画像）の画素毎の反射強度を表す。
e(h,τ,x)は、第３時刻（h）の衛星画像（撮像画像）に基づき第４時刻（h+τ）に対して推定された位置ｘ毎の反射強度と、第４時刻の衛星画像に基づく画素毎の反射強度との差分を表す。
Δτは、衛星画像が取得される時間間隔（時間フレームの時間間隔）のΔtと同じとするが、これに限定されない。

反射強度計算履歴算出部３３０は、数値気象計算過去履歴データベースから過去の数値気象計算の結果として、時刻h+τ、位置xにおける反射強度A_s(h+τ,x)を取得する。

反射強度計算誤差履歴算出部３４０は、以下の式により数値気象計算の反射強度の誤差を算出する。

e_s(h+τ,x)は、第４時刻（h+τ）に対する数値気象計算の前位置の反射強度と、第４時刻の衛星画像に基づく画素毎の反射強度との差分に対応する。

重み計算部３５０は、第７の実施形態と同様、以下の式により、リードタイムτ、位置xにおける重みw_τ(x)を計算する。

地点ごとに重みが大きく変化することは無いと仮定される場合、以下の式のように、全地点で平均した重みを用いてもよい。

合成部３６０は、雲移動推定装置１００によって推定された時間フレームt、リードタイムτ、位置xにおける反射強度A'(t,τ,x)を、第３重みに対応するw_τによって重み付けする。合成部３６０は、数値気象計算部２００によって計算された時刻t+τ、位置xにおける反射強度

を、第４重みに対応する1-w_τ(x)によって重み付けする。合成部３６０は、重み付けしたこれらの反射強度を合計（合成）することで、反射強度の新たな推定値A"(t,τ,x)を算出する。A"(t,τ,x)を算出する式を以下に示す。

つまり合成部３６０は、推定された各位置の反射強度に第３重みを設定し、数値気象計算の各位置の反射強度に第４重みを設定し、第３重み及び第４重みに基づき、反射強度を合成することにより、推定対象となる時刻（t+τ）における各位置の反射強度を推定する。

出力部は、合成部３６０により算出された新たな反射強度に基づく情報を出力（例えば画面に表示）する。

以上、本実施形態によれば、数値気象計算の結果を用いることで、雲の移流推定（特に雲の発達状態の推定）を高精度に行うことができる。

（第９の実施形態）
雲移動推定の重要な応用の一つとして、地表面への日射量の推定がある。近年、太陽光発電の普及が進んでいる。太陽光発電の場合、発電量が日射量に左右される。電力を有効又は効率的に活用するためには予め発電量を把握しておく必要があり、発電量の把握のために日射量の予測が必要となる。本実施形態では、第１の実施形態に係る雲移動推定装置１００に対して日射量推定部を追加することにより、日射量の推定を可能とする。

図２７は、本実施形態に係る情報処理装置として日射両推定装置の一例のブロック図である。

雲移動推定装置１００は、第１の実施形態に係る雲移動推定装置（図１参照）と同じである。

日射量計算部３８０は、雲移動推定装置１００により推定された時間フレームt、リードタイムτ、位置（地点）xにおける反射強度R(t,τ,x)に基づき、地点xに対応する地表面日射量（地点xの真下の地点における地表面日射量）を計算する。

雲の反射強度Rに基づき地表面日射量Sを計算する方法としていくつかの方法が知られている。例えば、比較的単純なモデル化を行っているが、以下の式を用いてもよい（参考文献８参照）。

Iは地球の大気の上端において、太陽からの放射線に対して垂直な面における単位面積・単位時間当たりに受けるエネルギー（仕事率）を表す。季節での変動はあるが、平均値は約1.37×10³(W/m²)である。

θは太陽の天頂角を表す。天頂角は、日時から計算できる。

A_Gは地表面アルベドである。地表面アルベドは、リモートセンシングデータから計算できる。

以上、本実施形態によれば、地表面への日射量を高精度に推定することができる。

（ハードウェア構成）
図２８は、各実施形態に係る雲移動推定装置（情報処理装置）のハードウェア構成を示す。雲移動推定装置は、コンピュータ装置６００により構成される。コンピュータ装置６００は、ＣＰＵ６０１と、入力インタフェース６０２と、表示装置６０３と、通信装置６０４と、主記憶装置６０５と、外部記憶装置６０６とを備え、これらはバス６０７により相互に接続されている。

ＣＰＵ（中央演算装置）６０１は、主記憶装置６０５上で、コンピュータプログラムである情報処理プログラムを実行する。情報処理プログラムは、雲移動推定装置の上述の各機能構成を実現するプログラムのことである。情報処理プログラムは、１つのプログラムではなく、複数のプログラムやスクリプトの組み合わせにより実現されていてもよい。ＣＰＵ６０１が、情報処理プログラムを実行することにより、各機能構成は実現される。

入力インタフェース６０２は、キーボード、マウス、及びタッチパネルなどの入力装置からの操作信号を、雲移動推定装置に入力するための回路である。入力インタフェース６０２は入力装置１２０に対応する。

表示装置６０３は、雲移動推定装置から出力されるデータを表示する。表示装置６０３は、例えば、ＬＣＤ（液晶ディスプレイ）、有機エレクトロルミネッセンスディスプレイ、ＣＲＴ（ブラウン管）、又はＰＤＰ（プラズマディスプレイ）であるが、これに限られない。コンピュータ装置６００から出力されたデータは、この表示装置６０３に表示することができる。

通信装置６０４は、雲移動推定装置が外部装置と無線又は有線で通信するための回路である。データは、通信装置６０４を介して外部装置から入力することができる。外部装置から入力したデータを、主記憶装置６０５や外部記憶装置６０６に格納することができる。

主記憶装置６０５は、情報処理プログラム、情報処理プログラムの実行に必要なデータ、及び情報処理プログラムの実行により生成されたデータなどを記憶する。情報処理プログラムは、主記憶装置６０５上で展開され、実行される。主記憶装置６０５は、例えば、ＲＡＭ、ＤＲＡＭ、ＳＲＡＭであるが、これに限られない。雲移動推定装置の各記憶部又はデータベースは、主記憶装置６０５上に構築されてもよい。

外部記憶装置６０６は、情報処理プログラム、情報処理プログラムの実行に必要なデータ、及び情報処理プログラムの実行により生成されたデータなどを記憶する。これらの情報処理プログラムやデータは、情報処理プログラムの実行の際に、主記憶装置６０５に読み出される。外部記憶装置６０６は、例えば、ハードディスク、光ディスク、フラッシュメモリ、及び磁気テープであるが、これに限られない。雲移動推定装置の各記憶部又はデータベースは、外部記憶装置６０６上に構築されてもよい。

なお、情報処理プログラムは、コンピュータ装置６００に予めインストールされていてもよいし、ＣＤ－ＲＯＭなどの記憶媒体に記憶されていてもよい。また、情報処理プログラムは、インターネット上にアップロードされていてもよい。

また、雲移動推定装置は、単一のコンピュータ装置６００により構成されてもよいし、相互に接続された複数のコンピュータ装置６００からなるシステムとして構成されてもよい。

なお、本発明は上記各実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。また、上記各実施形態に開示されている複数の構成要素を適宜組み合わせることによって種々の発明を形成できる。また例えば、各実施形態に示される全構成要素からいくつかの構成要素を削除した構成も考えられる。さらに、異なる実施形態に記載した構成要素を適宜組み合わせてもよい。

［参考文献リスト］
・参考文献１
尾関一頼, 佐々木幸男. 静止気象衛星ひまわり 8 号及び 9 号の概要 (ひまわり 8 号及び 9 号の地上シス
テム総合報告). 気象衛星センター技術報告, pp. 3-16, 2016.
・参考文献２
Maarten Reyniers. Quantitative precipitation forecasts based on radar observations: Principles, algo-
rithms and operational systems. Institut Royal M´et´eorologique de Belgique, 2008.
・参考文献３
Jerome Friedman, Trevor Hastie, Robert Tibshirani, et al. The elements of statistical learning, Vol. 1.
Springer series in statistics New York, 2001.
・参考文献４
Mryka Hall-Beyer. GLCM texture: a tutorial. National Council on Geographic Information and Analysis Remote Sensing Core Curriculum, Vol. 3, p. 75, 2000.
・参考文献５
M Lalitha, M Kiruthiga, and C Loganathan. A survey on image segmentation through clustering algorithm. International Journal of Science and Research, Vol. 2, No. 2, pp. 348-358, 2013.
・参考文献６
Benoit Cushman-Roisin and Jean-Marie Beckers. Introduction to geophysical fluid dynamics: physical and numerical aspects. Academic press, 2011.
・参考文献７
Gerhard Winkler. Image analysis, random fields and Markov chain Monte Carlo methods: a mathe-matical introduction, Vol. 27. Springer Science & Business Media, 2012.
・参考文献８
G Dedieu, PY Deschamps, and YH Kerr. Satellite estimation of solar irradiance at the surface of the earth and of surface albedo using a physical model applied to metcosat data. Journal of Applied Meteorology and Climatology, Vol. 26, No. 1, pp. 79-87, 1987.

１ 気象衛星
１０ 雲構造推定部
１１ 雲水平位置推定部
１２ 雲立体構造推定部
２０ 雲分解部
２１ 画像特徴量計算部
２２ 画像セグメンテーション部
２３ 統計モデル部
２４ 推定部
２６ 不確定部分特定部
２７ 移流推定合成部
３０ 推定処理部
３０＿１～３０＿Ｎ 雲時間変化推定部
３１ 移流ベクトル推定部
３２ 移流計算部
３３ 雲マッチング部
３４ 発達・減衰係数推定部
３５ 雲発達・減衰計算部
３６ 風移流ベクトル算出部
３７ 移流ベクトル合成部
４０ 雲合成部
５０ 雲分割パラメータ推定部
６０ 推定合成部
１００ 雲移動推定装置
１００Ａ 雲移動推定装置
１００Ｂ 雲移動推定装置
１１０ 入力部
１２０ 入力装置
１３０ 記憶部
１３０ 出力装置
２００ 数値気象計算部
２２０ 衛星画像データベース
２３０ 数値気象計算履歴データベース
２４０ データベース解析部
２４１ 雲構造推定部
２４２ 雲分解部
２４３ 移流ベクトル推定部
２５０ 移流ベクトル誤差履歴算出部
２６０ 風移流ベクトル履歴算出部
２７０ 風移流ベクトル誤差履歴算出部
２８０ 重み計算部
３１０ 過去履歴用雲移動推定装置
３２０ 反射強度推定誤差履歴算出部
３３０ 反射強度計算履歴算出部
３４０ 反射強度計算誤差履歴算出部
３５０ 重み計算部
３６０ 合成部
３８０ 日射量計算部
６００ コンピュータ装置
６０２ 入力インタフェース
６０３ 表示装置
６０４ 通信装置
６０５ 主記憶装置
６０６ 外部記憶装置
６０７ バス

Claims (24)

1. 雲をセンシングしたデータに基づき、前記雲の高さ情報を推定する第１推定部と、
   前記高さ情報に基づき、前記雲を複数の部分に分解する分解部と、
   前記部分ごとに対象時刻における前記雲の状態を推定する第２推定部と、
   を備えた情報処理装置。
2. 前記データは、前記雲の温度分布を含み、
   前記第１推定部は、前記雲の温度分布に基づき、前記高さ情報を推定する
   請求項１に記載の情報処理装置。
3. 前記温度分布は、前記雲の赤外画像である
   請求項２に記載の情報処理装置。
4. 前記分解部は、前記温度分布に基づき複数の温度範囲ごとに前記雲を分解することにより前記複数の部分を取得する
   請求項２又は３に記載の情報処理装置。
5. 前記第２推定部は、前記部分ごとに前記雲の移動の方向及び速さを推定し、前記移動の方向及び速さに基づいて、前記雲の位置及び形状の少なくとも一方を推定する
   請求項１～３のいずれか一項に記載の情報処理装置。
6. 前記部分ごとに推定された前記雲を合成して画像を生成する合成部
   を備えた請求項１～５のいずれか一項に記載の情報処理装置。
7. 前記データは、前記雲の撮像画像を含み、
   前記分解部は、前記雲の撮像画像に含まれる画素の特徴量に基づき前記撮像画像をセグメンテーションすることにより、前記雲を前記複数の部分に分解する
   請求項１～６のいずれか一項に記載の情報処理装置。
8. 前記分解部は、前記複数の部分のうちの第１部分の境界において前記第１部分に隣接する他の前記部分と隣接する範囲に基づいて、前記第１部分の下層において他の雲が存在するか否か又は他の雲が存在する範囲を推定する
   請求項１～７のいずれか一項に記載の情報処理装置。
9. 前記分解部は、第１時刻の前記データに基づき推定された第２時刻における前記部分ごとの前記雲の位置に基づき、前記第２時刻の前記データに基づき得られた前記雲の前記複数の部分の下層において、他の雲が存在するか否か又は他の雲が存在する範囲を推定する
   請求項１～７のいずれか一項に記載の情報処理装置。
10. 前記第２推定部は、第１時刻の前記データに基づく前記部分ごとの前記雲の発達状態と、第１時刻より後の第２時刻の前記データに基づく前記部分ごとの前記雲の発達状態との差異に基づいて、前記第２時刻より後の前記対象時刻における前記部分ごとの前記雲の発達状態を推定する
    請求項１～９のいずれか一項に記載の情報処理装置。
11. 前記分解部は、複数の候補パラメータに基づき、前記雲を複数の部分に分解し、
    前記候補パラメータごとに推定される前記雲の状態を統合する推定統合部をさらに備えた
    請求項１～１０のいずれか一項に記載の情報処理装置。
12. 前記高さ情報は、前記雲の温度分布を含み、
    前記分解部は、前記温度分布に基づき複数の温度範囲ごとに前記雲を分解することにより前記複数の部分を取得し、
    前記複数の候補パラメータは、前記複数の温度範囲を規定するパラメータである
    請求項１１に記載の情報処理装置。
13. 前記第２推定部は、前記部分ごとに前記雲の移動の方向及び速さを推定し、推定した前記移動の方向及び速さと、数値気象計算の風向及び風速とに基づき、前記雲の位置を推定する
    請求項１～１２のいずれか一項に記載の情報処理装置。
14. 前記第２推定部は、前記雲の前記部分ごとに推定した前記移動の方向及び速さの第１ベクトルに第１重みを設定し、数値気象計算による前記部分の位置に対応する風向及び風速の第２ベクトルに第２重みを設定し、前記第１重み及び前記第２重みに基づき、前記第１ベクトル及び前記第２ベクトルを合成することにより、前記部分ごとに前記雲の位置を推定する
    請求項１３に記載の情報処理装置。
15. 第３時刻の前記データに基づき推定した前記移動の方向及び速さと、前記第３時刻より後の第４時刻の前記データに基づき推定した前記移動の方向及び速さとの差分と、第４時刻の前記データに基づき推定した前記移動の方向及び速さと、前記第４時刻に対する前記数値気象計算の風向及び風速との差分とに基づいて、前記第１重み及び前記第２重みを計算する重み計算部
    を備えた請求項１４に記載の情報処理装置。
16. 前記雲をセンシングしたデータは前記雲を含む第１領域の撮像画像を含み、前記撮像画像の画素は前記画素に対応する位置の反射強度を表し、
    前記第２推定部は、前記対象時刻における前記第１領域の各位置の反射強度を推定し、
    推定した前記各位置の反射強度と、数値気象計算の前記第１領域の各位置の反射強度とを合成する合成部を備えた
    請求項１～１５のいずれか一項に記載の情報処理装置。
17. 前記合成部は、推定された前記各位置の反射強度に第３重みを設定し、前記数値気象計算の前記第１領域の各位置の反射強度に第４重みを設定し、前記第３重み及び前記第４重みに基づき、前記反射強度を合成する
    請求項１６に記載の情報処理装置。
18. 第３時刻の撮像画像に基づき第４時刻に対して推定された前記位置毎の反射強度と、前記第４時刻の撮像画像に基づく前記画素毎の反射強度との差分と、前記第４時刻の撮像画像に基づく前記画素毎の反射強度と、前記第４時刻に対する前記数値気象計算の前記各位置の反射強度との差分とに基づいて、前記第３重み及び前記第４重みを計算する重み計算部
    を備えた請求項１７に記載の情報処理装置。
19. 前記雲をセンシングしたデータは前記雲を含む第１領域の撮像画像を含み、前記撮像画像の画素は前記画素に対応する位置の反射強度を表し、
    前記第２推定部は、前記対象時刻における前記第１領域の各位置の反射強度を推定し、
    前記第２推定部により推定された前記各位置の反射強度に基づき、地表面における日射量を計算する日射量計算部
    を備えた請求項１～１８のいずれか一項に記載の情報処理装置。
20. 前記対象時刻は、前記雲がセンシングされた時刻よりも後の時刻である
    請求項１に記載の情報処理装置。
21. 前記第１推定部は、前記雲をセンシングした前記データに基づき前記雲の位置情報をさらに推定する
    請求項１～２０のいずれか一項に記載の情報処理装置。
22. 前記高さ情報は、前記雲の水平方向の位置毎の高さを含む
    請求項１～２１のいずれか一項に記載の情報処理装置。
23. 雲をセンシングしたデータに基づき、前記雲の高さ情報を推定し、
    前記高さ情報に基づき、前記雲を複数の部分に分解し、
    前記部分ごとに対象時刻における前記雲の状態を推定する
    情報処理方法。
24. 雲をセンシングしたデータに基づき、前記雲の高さの分布を推定するステップと、
    前記高さ情報に基づき、前記雲を複数の部分に分解するステップと、
    前記部分ごとに対象時刻における前記雲の状態を推定するステップと、
    をコンピュータに実行させるためのコンピュータプログラム。

Images