卒業研究テーマ 静止気象衛星画像による 風ベクトルの推定

卒業研究テーマ
静止気象衛星画像による
風ベクトルの推定
教育学部情報科学教育課程
1
静止気象衛星画像による風ベクトルの推定
1
はじめに
まずはじめに、現在の気象観測の概要を述べる。
最初の静止気象衛星は 1966 年の末にアメリカで打ち上げられた ATS-1(Applications Technology Satellite)
である。日本の静止気象衛星 GMS（Geostationary Meteorological Satellite) 通称 ”ひまわり ”の第１号は、
1978 年に米国から打ち上げられ、続く GMS-2 からは種子島宇宙センターからの打ち上げとなり３号、４号
を経て現在では５号が運用されている。なお、国際的にも気象界でも ”ひまわり ”というよりは GMS とい
う名称のほうが一般的である。GMS は、東経 140 °の赤道上、約 36,000km の静止衛星軌道に位置し、常
時宇宙から広い範囲の雲の分布などの気象観測をおこなうとともに、雲画像の配信を定常的におこなって
いる。GMS は現在、気象業務になくてはならない存在となり、特に台風、豪雨、豪雪などの大きな災害に
つながる異常気象時に、より迅速・的確な予・警報が出せるようになり、災害の防止、軽減に大きな役割を
果たしている。GMS の雲画像は、我が国のみならずアジア、オセアニアの 20 数カ国・領域でも受信され、
それぞれの国における気象業務の改善にも大きく寄与している。
GMS は、世界気象機関（WMO) の世界気象監視（WWW) 計画の一環として、WMO と国際学術連合
（ICSU) が共同して行う地球大気開発計画（GARP) の推進を図り、併せて我が国の気象業務の改善に資す
るために開発が開始されたものである。
図 1: 世界の気象衛星観測網
国際的に見ると、GMS(日本) や GOES ×２（米国）、METEOSAT（欧州気象衛星機構）、GOMS（ロ
シア）の５つの静止気象衛星と、地球を南北に周りながら気象観測をおこなういくつかの軌道衛星 (ロシア
の METEOR や米国の NOAA) によって、地球を取り巻くように位置し、世界気象衛星観測網を展開して
いる (図１）。
2
研究目的
GMS から発信される静止気象衛星画像をもとに、台風や雲の動きを推定、また、台風の軌道を予測する
ことにより、それらの情報から防災へと応用させることができる。さらに雲の動きを予測することによっ
て、短期気象予報に応用できる。そこで本研究では 2 種類の風ベクトルの推定手法を実現することを目的
とする。実現する手法のうち１つは相互相関関数による対応付けを用いた相関法で、もうひとつは時空間
勾配＋変分法を用いた勾配法である。
2
静止気象衛星画像による風ベクトルの推定
3
相関法による風ベクトルの推定
3.1
相関法
まずはじめに相関法について述べる。相関法とは、連続する２枚の画像間において、濃淡パターンの対応
付けをおこない、２次元的な雲の移動ベクトル場、を推定する方法である。濃淡パターンの画像平面上での
移動によって、2 枚の連続する画像間で濃淡の時間変化が生じる。つまりこのとき相関関数を用いることに
より、連続する２枚の画像間の類似度を測ることができる。そこで、次に相関関数についてまとめてみた。
3.2
相関関数
相関関数は、自己相関関数と相互相関関数に分けられるが、まずはじめに自己相関関数の定義と意味につ
いてまとめた。
一般的に二つの変量 x と y との相互の関連の度合いを定量的に調べるには、x − y 面上に多くのサンプルか
ら得られる (x, y) で定まる点を打てばよい。x と y が時間の関数ならば x − y 面上の連続した曲線となる。
もし、x と y とに関連があれば、一つの直線ないしは (曲線）のまわりに点は分布するし、逆に関連がなけ
れば原点のまわりに一様に分布するであろう。
二つの変量の相関度は、
· ¸
E[xy]
x
E[xy]
rx = E
=
; ry = E
=
2
x
E[x ]
y
E[y 2 ]
hyi
x と y とを同じように扱うためには、
E[xy]
r= p
E[x2 ]E[y 2 ]
あるいは C = E[xy]
を測ればよい。ここに、E はアンサンブル平均 (母集団平均）を意味する。もし、x と y が無関係ならば
r=0
であり、x と y とが α 倍の違いで完全に一致 (x ≡ αy) すれば明らかに
r = 1(または − 1, α < 0)
となる。
一方、x(t) が周期 T の周期変動であるとすれば、
x(t) = x(t ± nT ) (n = 0, 1, 2, · · ·)
つまり、周期の整数倍だけ時間をずらすと元の波形と重なってしまう。不規則変動 x(t) が周期性の強い
ものならば、周期の整数倍だけ時間軸をずらすと元の波形とかなり似ているであろう。したがって、ある時
間 τ だけずらした波形が元の波形とどれだけ似ているかを調べ変動中の周期成分を判別するには、上の場
合と同じく x = x(t) と y = x(t + τ ) の相関を求めればよい。
時間に関する不規則変量を x(t) とするとき、τ 時間隔たった二つの変動の積の平均値
C = E[xy] から x = x(t), y = x(t + τ ) を代入して
3
静止気象衛星画像による風ベクトルの推定
図 2: 自己相関関数の説明図
C(t, τ ) = E[x(t)x(t + τ )]
(1)
で定義される統計関数を自己相関関数と呼ぶ。
また、隔たり時間 τ をラグ (lag) という。ここでの平均操作は、原義的には ”アンサンブル平均 ”である
が、定常確立過程では多くの場合これを式 (2) のように時間平均でおきかえることができ、また時刻 t には
無関係でラグ τ のみの関数である。
Cx (τ ) =
=
x(t)x(t + τ )
Z T2
1
lim
x(t)x(t + τ )dt
T →∞ T − T
2
（ここで———は時間平均を意味する。)
Cx (τ ) を τ = 0 の値 Cx (0) で割って正規化したものを、自己相関係数と呼ぶ。
Cx (τ )
x(t)x(t + τ )
=
Cx (0)
x2 (t)
ÁZ
Z T2
=
x(t)x(t + τ )dt
Rx (τ ) =
− T2
Rx (0)
T
2
− T2
= 1
最も単純な例として正弦波の自己相関を求める。いま、x(t) を、
x(t) = a sin 2πf t
とすれば、x(t + τ ) は、
x(t + τ ) = a sin 2πf (t + τ )
4
x2 (t)dt
(2)
静止気象衛星画像による風ベクトルの推定
である。このとき、式 (2) の右辺の被積分関数は、
x(t)x(t + τ ) =
a sin 2πf t · a sin 2πf (t + τ )
=
a2 {sin 2πf t · sin 2πf (t + τ )}
=
a2 {sin 2πf t · (sin 2πf t · cos 2πf τ + cos 2πf t · sin 2πf τ )}
=
a2 {sin2 2πf t · cos 2πf τ + sin 2πf t · cos 2πf t · sin 2πf τ }
=
a2 {(1 − cos2 2πf t) · cos 2πf τ + cos 2πf t · (cos 2πf t · cos 2πf τ − cos 2πf (t + τ ))}
=
a2 {cos 2πf τ − cos2 2πf t · cos 2πf τ + cos2 2πf t · cos 2πf τ − cos 2πf t cos 2πf (t + τ )}
1
a2 {cos 2πf τ − {cos 2πf (t + t + τ ) + cos 2πf (t − t − τ )}}
2
1
1
2
a {cos 2πf τ − cos 2πf (2t + τ ) − cos 2πf (−τ )}
2
2
a2
{cos 2πf τ − cos(2πf (2t + τ ))}
2
=
=
=
である。上式を t = 0 から n 波分 (周期 T =
Z
Z
nT
x(t)x(t + τ )dt
nT
=
0
0
Z
nT
=
0
1
f)
にわたって積分すれば、
a2
{cos 2πf τ − cos(2πf (2t + τ ))}dt
2
Z nT 2
a2
a
(cos 2πf τ )dt −
cos(2πf (2t + τ ))dt
2
2
0
ここで右辺第二項に着目すると、
Z nT 2
a
−
cos(2πf (2t + τ ))dt
2
0
Z
a2 nT
= −
cos(2πf (2t + τ ))dt
2 0
1
nt
=
[sin(2πf (2t + τ ))]0
4πf
T = f1 、n = 0, 1, 2, · · · であることより sin 4πn = 0 となるから、
Z
nT
x(t)x(t + τ )dt =
0
nT a2
(cos 2πf τ )
2
したがって、自己相関係数 Rx (τ ) は、
Á Z nT
x2 (t)dt
x(t)x(t + τ )dt
0
0
Á
nT a2
nT a2
=
(cos 2πf τ )
2
2
= cos 2πf τ
Z
nT
Rx (τ ) =
平均操作を時間平均ではなく、式 (1) のようにアンサンブル平均をとる場合は次のようになる。
すなわち、第 i 番目の標本を、
xi (t) = a sin(2πf t + θi )
とする。正弦波の波形は同じであるが、位相は標本ごとに θi だけずれている。このとき、
5
静止気象衛星画像による風ベクトルの推定
Cx (τ ) =
=
=
E[xi (t)xi (t + τ )]
N
hX
± i
lim
xi (t)xi (t + τ ) N
N →∞
lim
i=1
N
hX
N →∞
i=1
h
=
=
lim
N →∞
± i
a sin(2πf t + θi ) · a sin(2πf (t + τ ) + θi N
a2
N
X
± i
cos(2πf t + θi ) · cos(2πf (t + τ ) + θi ) − cos(2πf (2t + τ ) + 2θi ) N
i=1
i
a2
a2 h
cos 2πf τ − E cos(2πf (2t + τ ) + 2θi
2
2
いま、位相角 θi は、(0, 2π) 間に確立 p(θ) で不規則に分布するとすれば、
h
i Z 2π
E cos(2πf (2t + τ ) + 2θi =
cos(2πf (2t + τ ) + 2θ)p(θ)dθ
0
である。確率分布が一様である (p(θ) = 1/2π) とすると、上式は 0 となる。
したがって、正弦波の自己相関関数と自己相関係数はそれぞれ次のようになる。
a2
cos 2πf τ
2
Cx (τ ) =
Rx (τ ) = cos 2πf τ
自己相関関数の一般的性質として、自己相関関数は偶関数である。
式 (2) において τ = −τ1 とすれば、
1
Cx (−τ1 ) = lim
T →∞ T
Z
T
2
− T2
x(t)x(t − τ1 )dt
上式において積分変数を t1 = t − τ1 とおけば
1
Cx (−τ1 ) = lim
T →∞ T
Z
T
2−τ1
T
− 2−τ
x(t1 )x(t1 + τ1 )dt1
1
区間 (−T /2, T /2) と区間 (−T /2 − τ1 , T /2 − τ1 ) とは lim の極限では一致するから上式の右辺の積分区間は
(−T /2, T /2) に置き換えることができる。したがって、Cx (−τ ) は偶関数である。
Cx (τ ) = Cx (−τ )
次に、相互相関関数については、自己相関関数を 2 次元的に考えたものである。つまり、時間に関するあ
る変動量 x(t) に対し、τ 時間隔たった時刻の変動量 y(t + τ ) との積の平均値で定義される関数を相互相関
関数という。
相互相関関数を次のように定義する。
Cxy (τ ) =
=
x(t)y(t + τ )
Z T2
1
lim
x(t)y(t + τ )dt
T →∞ T − T
2
6
(3)
静止気象衛星画像による風ベクトルの推定
図 3: 相互相関関数の説明図
（ここで———は時間平均を意味する。)
Cxy (τ ) を τ = 0 の値、つまり Cxy (0) で割って、正規化したものを次の Rxy (τ ) で定義する。
=
Áp q
x(t)y(t + τ )
x2 y 2
Áq
Cxy (τ )
Cx (0)Cy (0)
=
1
Rxy (τ ) =
Rxy (0)
この場合、画像と画像の対応づけという処理の中で、この Rxy (τ ) の値が 1 に近ければ近いほど、画像の
濃淡パターンが似かよっているということになるため、台風や雲の動きを追うことができ、台風や雲の移動
ベクトルを推定することができるというわけである。
一方、式 (3) のような x(t) と y(t) の積による相関ではなく、差の絶対値を用いた次の相関関数によって
も相似度を測ることができる。
Z
Cxy (τ ) =
T
2
|x(t) − y(t + τ )|dt
(4)
− T2
我々の研究では、この差の絶対値を用いた相関関数を利用して、C 言語によるプログラムを実現した。付録
1 と付録 2 でその基礎となるプログラムを挙げ、付録 3 に風ベクトルを推定するプログラムを挙げる。
3.3
相関関数を用いた風ベクトル推定法
では C 言語のプログラムによる、台風や雲の移動ベクトルの推定とは具体的にどういう考え方で行われ
るものかを図 4 を用いて解説する。まず連続する 2 枚の画像を仮定する。左側を時刻 t = 0、右側を t = 1
とする。図中の２重丸が台風であるとする。t = 0 における台風が t = 1 のとき画像中のどこへ移動したか
を推定する。まず、t = 0 のときの台風を囲む部分を、t = 1 の画像中全てに対応付けていく。このとき差
の絶対値による相関関数を利用するため、ある部分の値が最も 0 に近いはずである。つまり濃淡パターン
が最も似かよっているということになり、t = 0 における台風が、その場所へ移動したということがわかる
のである。ここで、台風のまわりをどれくらいの大きさで囲み、対応付けていくかということは、観測す
る台風や雲の大きさがその時々によって異なるため、あらかじめプログラム中で指定しておく必要がある。
今回の研究では、画像全体のサイズが 2291 × 2291 に対し、囲む範囲を台風の大きさに合わせ、150 × 120
7
静止気象衛星画像による風ベクトルの推定
図 4: 相関関数による雲の移動速度の推定
とした。これと同じ処理を t = 1 と t = 2 で、さらに t = 2 と t = 3、というふうに繰り返すことにより雲
や台風の移動ベクトルを推定していくことができるのである。
3.4
実行結果
図５に示す左側の気象衛星は時刻 t = 0 のもである。今回の研究では、この時刻 t = 0 の画像から時刻
t = 4 までの合計 5 枚の静止気象衛星画像を用いて風ベクトルを推定した。
右側のグラフは、実現した先ほどのプログラムの実行結果から導かれた座標のプロットと、使用した連続す
る 5 枚の画像中の台風の目の座標を 1 つ 1 つ手動で確認し、プロットしたものを比較したものである。
多少の誤差はあるがほぼ同じ形の軌跡をたどっていることがわかる。
以上が相関法による風ベクトルの推定である。
なお、ここで使用した静止気象衛星画像は NASA ホームページよりダウンロードしたものである。
図 5: 実行結果
4
勾配法による風ベクトルの推定
次に、勾配法を用いた風ベクトルの推定法を示す。勾配法とは，パターンの特徴を表す画像関数（濃淡
分布）f (x, y, t) が，運動に際し不変に保たれるとの仮定 (f (x, y, t) = f (x + δ, y + δ, t + δ)) より，ある点
(x, y) における速度と，動画像の濃淡分布の空間勾配および時間勾配とを関係付ける式（拘束条件）をもと
にした解析手法である。解析の対象となる気象衛星画像は，ほぼ一定時間間隔で１枚ずつサンプリングさ
れた濃淡画像である。勾配法を濃淡動画像に適用する場合，基本的に各画素で独立にオプティカルフローを
決定することは困難である。そこで，何らかの第２の拘束条件式が必要となる。
8
静止気象衛星画像による風ベクトルの推定
4.1
基本式の導出
図 6 のように，時刻 t において座標 (x, y) にある濃淡パターンが，δt 時刻経過した後に座標 (x + δx, y + δy)
に，その濃淡値分布を一定に保ったまま移動したとき，その対応付けから次式が成り立つ。
ρ(x, y, t) = ρ(x + δx, y + δy, t + δt)
(5)
右辺をテイラー展開し整理すると，速度ベクトル (u, v) を用いて，次式となる。
ρx u + ρy v + ρt = 0
(6)
ここで、
ρx = ∂ρ/∂x, ρy = ∂ρ/∂y, ρt = ∂ρ/∂t,
u = δx/δt, v = δy/δt
とした。
これが，勾配法で一般に用いられる基本式である。これに，速度ベクトルの空間的変化を最小にする次の
条件式を付け加える。
(ux )2 + (uy )2 + (vx )2 + (vy )2 → min
(7)
式（6）と式（7）を用いて，誤差 E を次のように定義する。ここで、
ux
= ∂u/∂x, uy = ∂u/∂y,
vx
= ∂v/∂x, vy = ∂v/∂y
とした。
Z
{(ρx u +ρy v + ρt )2
E=
S
+α2 (u2x + u2y + vx2 + vy2 )}ds
(8)
α は式（6）と式（7）との相対的な重みを決定する。この評価関数 E を変分法を用いて，解を推定する。
図 6: 濃淡パターンの追跡の移動と速度ベクトル
9
静止気象衛星画像による風ベクトルの推定
4.2
変分法（１次元の場合）
図 7 のように最適曲線 x の変分を考える。
x̄(t) = x(t) + εη(t)
(9)
ここで ε は微小な定数， η は任意の微分可能な関数で，次式のようにその端点は固定されている。
η(t0 ) = η(tf ) = 0
(10)
図 7: 最適曲線 x とその近傍 x̄
関数 x̄ に対する誤差を J(ε) と書く。
Z
tf
J(ε) =
¯ t)dt
f0 (x, ẋ,
(11)
t0
式 (11) の停留条件は、次式をとる。
¯
µ
¶¸
Z tf ·
∂J ¯¯
∂f0
d
∂f0
=
−
η(t)dt = 0
∂ε ¯ε=0
∂x
dt
∂ ẋ
t0
(12)
ここで η(t) は任意であるので，上式の積分値が零になるためには [・] 内が零にならなければならない。す
なわち，次のオイラー・ラグランジュの方程式を得る。
µ
¶
∂f0
d
∂f0
−
=0
∂x
dt
∂ ẋ
(13)
式 (13) を解けば最適曲線 x が決定できる。
4.3
変分法による勾配法
Z Z
F (x, y, u, v, ux , uy , vx , vy )dxdy
G
10
(14)
静止気象衛星画像による風ベクトルの推定
が与えられたとする。これに、
½
u = u0 + tη
v = v0 + tξ
を代入し，１次元のときと同様に停留条件を適用すると、
Z Z
d
[
F (x, y, u0 + tη, v0 + tξ, u0x + tηx , v0x + tξx , u0y + tηy , v0y + tξy )dxdy]t=0 = 0
dt
G
(15)
ここで，Fux , Fvx は F (x, y, u, v, ux , vx , uy , vy ) の８個の独立な変数 x, y, u, v, ux , vx , uy , vy の関数とみて，
これを ux , vx について偏微分したものを表す。Fuy , Fvy についても同様に、
Z
Z
Z
2
(2α2 uy )ds Z
Z
Z
+ 2ρy (ρx u + ρy v + ρt) + (2α2 vx )ds + (2α2 vy )ds = 0
2ρx (ρx u + ρy v + ρt)ds +
S
(2α ux )ds +
S
S
S
S
(16)
S
第２項にグリーンの公式を適用する。３，５，６項についても同様である。
R
S
=
=
(2α2 ux )ds Z
Z Z
∂
G(2α2 ux )dy −
(2α2 ux )dxdy
∂x
∂
G
Z Z
∂
2
−
(2α ux )dxdy
G ∂x
同様に，３，５，６項も計算すると，
Z Z
{
G
+
∂
∂
(2α2 ux ) −
(2α2 uy )}η
∂x
∂x
∂
∂
{2ρy (ρy u + ρy v + ρt ) −
(2α2 vx ) −
(2α2 vy )}ξdxdy = 0
∂x
∂x
2ρx (ρx u + ρy v + ρt) −
となる。
ここで，η, ξ は任意の関数であるから，これが成立するためにはそれぞれの｛｝内が，＝０でなければな
らない。 すなわち，u = u0 , v = v0 のとき，
ρ2x u + ρx ρy v + ρx ρt − α2 ∇2 u = 0
したがって
ρ2x u + ρx ρy v = α2 ∇2 u − ρx ρt
これを解くと，風ベクトルが求まる。
11
静止気象衛星画像による風ベクトルの推定
4.4
グリーンの定理
自分自身に交わらない閉曲線を単一閉曲線という。D が，単一閉曲線 C で囲まれた xy 面上の領域とし，
関数 f (x, y)，g(x, y)，並びに，それらの１階偏導関数が，D および境界 C 上で連続であるとき，次式が成
り立つ。
Z Z
Z
(f dx + gdy) =
(
C
D
∂g
∂f
−
)dxdy.
∂x ∂y
(17)
ただし，C 上の線積分は，それに沿って積分するとき，領域 D を左側に見るように進むものとする。C
を回るこの向きを，正の向きという。これを グリーンの定理 という。
図 8: 単一閉曲線 C
証明
図 8 に示されているように，x 軸に平行な直線 M と C との交点は２個だけとし，その x 座標を x1 (y), x2 (y)
とする。また，M と C の接点を P1 , P3 とし，その y 座標を c, d とする。このとき領域 D の点は，c ≤ y ≤
d; x1 (y) ≤ x ≤ x2 (y) の範囲を動く。さて，∂g/∂x の D 上での積分は，最初に x について，それから y に
ついて行えば，つぎのようになる。
Z Z
D
∂g
dxdy
∂x
Z
Z
d
=
x2 (y)
dy
Z
c
=
x1 (y)
∂g
dx =
∂x
Z
d
{g[x2 (y), y] − g[x1 (y), y]}dy
c
Z
g[x2 (y), y]dy −
Z P1 P2 P3
=
g[x1 (y), y]dy
Z P1 P 2 P3
g[x2 (y), y]dy +
P1 P2 P3
g[x1 (y), y]dy.
P1 P 2 P3
ここで，最後の式をまとめると，これは関数 g[x(y), y] の値を C に沿って一周積分したものになる。曲
線 C を、
½
x = x(y) =
x1 (y) : P3 P4 P1 上
x2 (y) : P1 P2 P3 上
12
静止気象衛星画像による風ベクトルの推定
と表せば，結局、
Z Z
D
∂g
dxdy =
∂x
Z
Z
g[x(y), y]dy =
P1 P2 P3 P4 P1
g(x, y)dy
(18)
C
となる。最後の式は，簡単のため，x(y) は C 上の値をとるという約束のもとに，x と書いてある。
つぎに，y 軸に平行な直線 L と C の交点も２個だけとする。このとき，∂f /∂y の積分は，最初 y につい
て，ついで x について行えば，上と全く同じようにして，結局
Z Z
Z
∂f
dxdy = −
f (x, y)dx
D ∂y
C
(19)
となる。式 (18) と (19) の差をつくれば，グリーンの定理が得られる。
4.5
実行結果
左側の画像が基にした気象衛星画像である。なお，この気象衛星画像は NASA のホームページからダウ
ンロードしたものである。
この画像を勾配法を用いたプログラムによって風ベクトルの向きと大きさを推定したものが右の図である。
そのプログラムを付録 4 に示す。
気象衛星画像より，台風は上へと移動していた。プログラムの実行結果を見ても，風ベクトルを表す矢印が
上方向を向いている。よって，本研究で作成した勾配法を用いたプログラムにより，風ベクトルをほぼ正し
く推定できたこととなる。
以上が，勾配法による風ベクトルの推定法である。
5
まとめ
最後にまとめとして、私たちは「静止気衛星画像からの風ベクトルの推定」というテーマで 1 年間研究
をしてきた。そして本研究では気象衛星画像からの風ベクトルの推定をおこなうため、相関関数や変分法
13
静止気象衛星画像による風ベクトルの推定
を取り入れた解析法を実現した。現在の気象予報はコンピュータによる数値予報が主流である。数値予報と
は大気に関する数理モデルを計算機を用いて数値的に解き、一定時間後の気象を予報するものである。そ
の数値予報のためには気温や湿度、そして本研究の主題でもある風ベクトルの方向や大きさなどの情報を
初期値として入力する。そして気象予報には正確性が求められる。正確な気象予報を行うためには、正確な
先に述べた初期値が必要とされる。
本研究の今後の課題としては、より正確な風ベクトルの推定のために、さらに精度のよい推定法の研究を
進めることが挙げられる。それにより、気象予報の情報が正確で信頼性のあるものとなる。
謝辞
本研究では、GMS の画像をアクティブアーカイブセンターよりダウンロードした。
我々は、ゴッダードスペースフライトセンターのアクティブアーカイブセンターのデータ配布に対して感謝
の意を表します。
参考・引用文献
・
「気象研究ノート」第１４８号 (１９８３)
日本気象学会
・
「パソコンによる動画像処理」
三池 秀敏 他
・
「宇宙開発事業団（NASDA) ホームページ」
(http://www.nasda.go.jp/)
・
「NASA ホームページ」
(http://www.nasa.gov/)
14
静止気象衛星画像による風ベクトルの推定
付録 1
#include <stdio.h>
#include <stdlib.h>
int main()
{
FILE *fp_in;
FILE *fp_out;
int i,xs,ys;
char sfn[512];
char gms_data[763*763];
printf("file name1?");
fgets(sfn,512,stdin);
for(i=0; i<512;++i) if(sfn[i]==’\n’) sfn[i]=0;
fp_in=fopen(sfn,"rb");
fclose(fp_in);
printf("output file?");
fgets(sfn,512,stdin);
for(i=0; i<512;++i) if(sfn[i]==’\n’) sfn[i]=0;
fp_out=fopen(sfn,"wb");
fwrite(gms_data,1,xs*ys,fp_out);
printf("file name2?");
fgets(sfn,512,stdin);
for(i=0; i<512;++i) if(sfn[i]==’\n’) sfn[i]=0;
fp_in=fopen(sfn,"rb");
fread(gms_data,1,xs*ys,fp_in);
fclose(fp_in);
fwrite(gms_data,1,xs*ys,fp_out);
fclose(fp_out);
}
15
静止気象衛星画像による風ベクトルの推定
付録 1 のプログラムを簡潔に説明すると、2 枚の静止気象衛星画像データを、1 つのファイルに格納する
というプログラムである。このプログラムの処理の流れとしては、まず 1 枚目の画像データを読み込む、次
に 2 枚の画像を格納するための新規ファイルをのファイル名を指定する。そのファイルに先ほど読み込ん
だ 1 枚目の画像データを書き込む。次に 2 枚目の画像データを読み込み、1 枚目を書き込んだファイルと同
じファイルに 2 枚目の画像データも書き込む。これがこのプログラムの処理の流れである。つまり、2 枚の
静止気象衛星画像が 1 つのファイルに格納されたのである。
さらにこのプログラムを基礎として、2 枚の気象衛星画像間において、ある点の移動ベクトル (vx, vy と
する ) を求めるプログラムを付録 2 に示す。
付録 2
#include <stdio.h>
#include <math.h>
#define WX_SIZE 10
#define WY_SIZE 10
#define MAX_YSPEED 10
#define MAX_XSPEED 10
#define STR_LEN 512
int main()
{
int x,y,t,i,j,x_size,y_size;
int vx,vy,*dim;
int x0,y0;
char str[STR_LEN];
unsigned char *f0,*f1;
FILE *fp;
puts("x size=?");
x_size=atoi(fgets(str,STR_LEN,stdin));
puts("y size=?");
y_size=atoi(fgets(str,STR_LEN,stdin));
printf("image size:x=%d,y=%d\n",x_size,y_size);
f0=(unsigned char *)malloc(x_size*y_size*sizeof*f0);
f1=(unsigned char *)malloc(x_size*y_size*sizeof*f1);
puts("file name=?");
fgets(str,STR_LEN,stdin);
16
静止気象衛星画像による風ベクトルの推定
for(i=0;i<STR_LEN;++i)
if(str[i]==’\n’|| str[i]==’\r’)
str[i]=0;
fp=fopen(str,"rb");
fread( f0, sizeof *f0,x_size*y_size, fp ) ;
fread( f1, sizeof *f1,x_size*y_size, fp ) ;
fclose(fp);
puts("x0=?");
x0=atoi(fgets(str,STR_LEN,stdin));
puts("y0=?");
y0=atoi(fgets(str,STR_LEN,stdin));
det(x0,y0, x_size, y_size, &vx,&vy,f0,f1) ;
printf("vx=%d, vy=%d\n",vx,vy);
return 0;
}
int cor(x0,y0,x,y,x_size,y_size,f0,f1)
int x_size,y_size;
int x0,y0,x,y;
unsigned char *f0,*f1;
{
if ( chk_pos(x0,y0,x_size,y_size)==-1 ||chk_pos(x,y,x_size,y_size)==-1 )
return -1;
else {
int i,j,p0,p1 ;
double c;
c=.0;
for (j=-WY_SIZE/2;j<=WY_SIZE/2;++j)
for(i=-WX_SIZE/2;i<=WX_SIZE/2;++i) {
p0=x_size*(y0+j)+x0+i;
p1=x_size*(y+j)+x+i;
c+=fabs((double)f0[p0]-(double)f1[p1]);
}
return (c/(double)(WX_SIZE*WY_SIZE));
}
}
int det(x0, y0, x_size, y_size,pvx, pvy,f0,f1)
17
静止気象衛星画像による風ベクトルの推定
int x0, y0;
int x_size, y_size;
int *pvx, *pvy;
unsigned char *f0,*f1;
{
int x,y,i,j;
int min_x,min_y ;
double c, min_c;
min_c=1.e+20 ;
for(y=y0-MAX_YSPEED; y<=y0+MAX_YSPEED;++y)
for(x=x0-MAX_XSPEED; x<=x0+MAX_XSPEED;++x) {
if((c=cor(x0,y0,x,y,x_size,y_size,f0,f1))<min_c){
min_c=c;
min_x=x;
min_y=y;
}
}
printf( "min_c=%e, min_x=%d, min_y=%d\n",min_c, min_x, min_y ) ;
*pvx=(min_x-x0);
*pvy=(min_y-y0);
return 0;
}
int chk_pos(x,y,x_size,y_size)
int x,y,x_size,y_size;
{
if(x-WX_SIZE/2<0 || x+WX_SIZE/2>=x_size || y-WY_SIZE/2<0 || y+WY_SIZE/2>=y_size)
return -1;
else
return 0;
}
付録 2 のプログラムの処理の流れを解説する。まず、使用する気象衛星画像のサイズを指定する。
ここでは ”x size” と ”y size” である。使用する画像ファイルは、先ほどのプログラムを利用して、2 枚の
画像をあらかじめ 1 つのファイルにまとめたものを使用する。そして x0, y0 という座標を入力することに
より、1 枚目と 2 枚目の画像間で、差の絶対値を利用した相関関数により濃淡パターンを対応付け、vx, vy
という移動ベクトルを導くのである。
また、プログラム中に、f 0, f 1 という配列を確保しておく必要がある。これがこのプログラムの処理の流れ
である。
以上の考えを元に、最終的にこの研究で、相関法による風ベクトルを推定するプログラムを実現させた。
18
静止気象衛星画像による風ベクトルの推定
付録 3
#include <stdio.h>
#include <math.h>
#define WX_SIZE 150
#define WY_SIZE 120
#define MAX_YSPEED 100
#define MAX_XSPEED 100
#define STR_LEN 512
int main()
{
int n,m,i,x_size,y_size;
int px,py;
int x0,y0;
char str[STR_LEN], fn[STR_LEN];
char f0[STR_LEN],f1[STR_LEN] ;
unsigned char *image0,*image1;
FILE *fp,*fpd;
puts("x size=?");
x_size=atoi(fgets(str,STR_LEN,stdin));
puts("y size=?");
y_size=atoi(fgets(str,STR_LEN,stdin));
printf("image size:x=%d,y=%d\n",x_size,y_size);
puts("x0=?");
x0=atoi(fgets(str,STR_LEN,stdin));
puts("y0=?");
y0=atoi(fgets(str,STR_LEN,stdin));
printf("point:x0=%d,y0=%d\n",x0,y0);
image0=(unsigned char *)malloc(x_size*y_size*sizeof*image0);
image1=(unsigned char *)malloc(x_size*y_size*sizeof*image1);
puts("file name?");
fgets(f0,512,stdin);
19
静止気象衛星画像による風ベクトルの推定
printf("file name=%s\n",f0);
for(i=0;i<STR_LEN;++i)
if(f0[i]==’\n’||f0[i]==’\r’)
f0[i]=0;
fp=fopen(f0,"r");
printf("file name=%s\n",f0);
fgets(fn,512,fp);
for(i=0;i<STR_LEN;++i)
if(fn[i]==’\n’||fn[i]==’\r’)
fn[i]=0;
printf("file name=%s\n",fn);
fpd=fopen(fn,"rb");
fread(image0,sizeof*image0,x_size*y_size,fpd);
fclose(fpd);
for(n=0;n<4;++n){
printf("n=%d\n",n);
fgets(fn,STR_LEN,fp);
for(i=0;i<STR_LEN;++i)
if(fn[i]==’\n’||fn[i]==’\r’)
fn[i]=0;
printf("file name=%s\n",fn);
fpd=fopen(fn,"rb");
fread(image1,sizeof*image1,x_size*y_size,fpd);
fclose(fpd);
det(x0,y0, x_size, y_size, &px,&py,image0,image1) ;
x0=px;
y0=py;
printf("px=%d, py=%d\n",px,py);
for(m=0;m<x_size*y_size;++m)image0[m] = image1[m] ;
}
return 0;
}
20
静止気象衛星画像による風ベクトルの推定
int cor(x0,y0,x,y,x_size,y_size,image0,image1)
int x_size,y_size;
int x0,y0,x,y;
unsigned char *image0,*image1;
{
if ( chk_pos(x0,y0,x_size,y_size)==-1 ||chk_pos(x,y,x_size,y_size)==-1 )
return -1;
else {
int i,j,p0,p1 ;
double c;
c=.0;
for (j=-WY_SIZE/2;j<=WY_SIZE/2;++j)
for(i=-WX_SIZE/2;i<=WX_SIZE/2;++i) {
p0=x_size*(y0+j)+x0+i;
p1=x_size*(y+j)+x+i;
c+=fabs((double)image0[p0]-(double)image1[p1]);
}
return (c/(double)(WX_SIZE*WY_SIZE));
}
}
int det(x0, y0, x_size, y_size, ppx,ppy,image0,image1)
int x0, y0;
int x_size, y_size;
int *ppx, *ppy;
unsigned char *image0,*image1;
{
int x,y;
int min_x,min_y;
double c, min_c;
min_c=1.e+20 ;
for(y=y0-MAX_YSPEED; y<=y0+MAX_YSPEED;++y)
for(x=x0-MAX_XSPEED; x<=x0+MAX_XSPEED;++x) {
if((c=cor(x0,y0,x,y,x_size,y_size,image0,image1))<min_c){
min_c=c;
min_x=x;
min_y=y;
}
}
printf( "min_c=%e, min_x=%d,min_y=%d\n", min_c, min_x, min_y ) ;
*ppx=min_x;
21
静止気象衛星画像による風ベクトルの推定
*ppy=min_y;
return 0;
}
int chk_pos(x,y,x_size,y_size)
int x,y,x_size,y_size;
{
if((x-WX_SIZE/2)<0 || (x+WX_SIZE/2)>=x_size ||(y-WY_SIZE/2)<0 || (y+WY_SIZE/2)>=y_size)
return -1;
else
return 0;
}
このプログラムの流れとしては、使用する連続した何枚かの画像（今回の研究では 5 枚）のファイル名の
みを、あらかじめ別のファイルに格納しておく。プログラムを実行した時に入力する初期値としては、使用
する気象衛星画像の縦、横のサイズ、移動ベクトルを推定するための中心座標（例えば台風の目の座標）、
画像ファイル名を格納したファイルのファイル名などである。
このプログラムを解説していく。
プログラム中で image0 と image1 という配列を用意した。処理の流れとしては、image0 に、格納して
いたファイル名から、先頭の 1 枚目の画像データを読み込み、image1 に 2 枚目の画像データを読み込む。
ここで image0 と image1 において、相関関数による対応付けの処理を行う。
その１つの処理が終わると、次に image1 にあるデータを image0 へ移し、image1 に次の画像データを読
み込む。image1 にあるデータを image0 へ移すことにより、読み込むデータが image1 の 1ヶ所で済むの
である。
そして再び相関関数による対応付け、というふうに繰り返していくことにより、最終的に台風の目が移動し
ていく座標を得ることができ、それにより移動ベクトルを推定することができるのである。
付録 4
#include <stdio.h>
#include <stdlib.h>
#define STR_LEN 512
/*
#define X_SIZE 763
#define Y_SIZE 763
*/
int read_data( image, x_size, y_size )
unsigned char *image ;
int
x_size, y_size;
{
int i ;
char fn[STR_LEN] ;
22
静止気象衛星画像による風ベクトルの推定
FILE *fp ;
puts( "gms image data = ?" ) ;
fgets( fn, STR_LEN, stdin ) ;
for ( i=0; i<STR_LEN; ++i ) if ( fn[i]==’\n’ ) fn[i]=0 ;
fp = fopen( fn, "rb" ) ;
fread( image, sizeof *image, y_size*x_size*2, fp ) ;
fclose( fp ) ;
return 0 ;
}
int write_data( vx, vy, x_size, y_size)
float *vx, *vy ;
int x_size, y_size;
{
int i ;
char fn[STR_LEN] ;
FILE *fp ;
puts( "velocity vector field = ?" ) ;
fgets( fn, STR_LEN, stdin ) ;
for ( i=0; i<STR_LEN; ++i ) if ( fn[i]==’\n’ ) fn[i]=0 ;
fp = fopen( fn , "wb" ) ;
for ( i=0; i<y_size*x_size; ++i ) {
fwrite( &(vx[i]), sizeof *vx, 1, fp ) ;
fwrite( &(vy[i]), sizeof *vy, 1, fp ) ;
}
fclose( fp ) ;
return 0 ;
}
double av( x, y, v, x_size )
int x, y;
float *v ;
int x_size;
{
double a ;
23
静止気象衛星画像による風ベクトルの推定
int pos ;
pos = y*x_size + x ;
a = ( v[pos-x_size] + v[pos+1]
+ v[pos+x_size] + v[pos-1] ) / 6.0
+ ( v[pos-x_size-1] + v[pos-x_size+1]
+ v[pos+x_size+1] + v[pos+x_size-1] ) / 12.0 ;
return a ;
}
int fxyt( image, fx, fy, ft, x_size, y_size )
unsigned char *image ;
double *fx, *fy, *ft ;
int x_size, y_size;
{
int x, y, pos0, pos1 ;
for ( y=1 ; y<y_size-1 ; ++y)
for ( x=1 ; x<x_size-1 ; ++x) {
pos0 = y*x_size+x ;
pos1 = pos0 + y_size*x_size ;
fx[pos0] = ( (int)image[pos0+1]
-(int)image[pos0]
+ (int)image[pos0+x_size+1]
-(int)image[pos0+x_size]
+ (int)image[pos1+1]
-(int)image[pos1]
+ (int)image[pos1+x_size+1]
-(int)image[pos1+x_size] ) / 4. ;
fy[pos0] = ( (int)image[pos0+x_size]
-(int)image[pos0]
+ (int)image[pos0+x_size+1]
-(int)image[pos0+1]
+ (int)image[pos1+x_size]
-(int)image[pos1]
+ (int)image[pos1+x_size+1]
-(int)image[pos1+1] ) / 4. ;
24
静止気象衛星画像による風ベクトルの推定
ft[pos0] = ( (int)image[pos1]
-(int)image[pos0]
+ (int)image[pos1+1]
-(int)image[pos0+1]
+ (int)image[pos1+x_size]
-(int)image[pos0+x_size]
+ (int)image[pos1+x_size+1]
-(int)image[pos0+x_size+1] ) / 4. ;
}
return 0 ;
}
int determine_of( image, vx, vy, x_size, y_size )
unsigned char *image ;
float *vx, *vy ;
int x_size, y_size;
{
char str[STR_LEN] ;
int i, in, n, pos, x, y ;
double ap, d, avx, avy, tt ;
double *fx, *fy, *ft ;
fx = (double *)malloc(x_size*y_size*sizeof *fx) ;
fy = (double *)malloc(x_size*y_size*sizeof *fy) ;
ft = (double *)malloc(x_size*y_size*sizeof *ft) ;
puts( "Iterative number = ?" ) ;
fgets( str, STR_LEN, stdin ) ;
in = atoi( str ) ;
puts( "Alpha = ?" ) ;
fgets( str, STR_LEN, stdin ) ;
ap = atof( str ) ;
for ( i=0; i<y_size*x_size; ++i ) {
vx[i] = .0 ; vy[i] = .0 ;
fx[i] = .0 ; fy[i] = .0 ; ft[i] = .0 ;
}
fxyt( image, fx, fy, ft, x_size, y_size ) ;
25
静止気象衛星画像による風ベクトルの推定
for ( n=0; n<in; ++n ) {
printf("n=%d\n", n ) ;
for ( y=1; y<y_size-1; ++y )
for ( x=1; x<x_size-1; ++x ) {
avx = av( x, y, vx, x_size ) ;
avy = av( x, y, vy, x_size ) ;
pos = y*x_size + x ;
d = fx[pos]*avx + fy[pos]*avy + ft[pos] ;
tt = ap*ap + fx[pos]*fx[pos] + fy[pos]*fy[pos] ;
vx[pos] = avx-fx[pos]*d/tt ;
vy[pos] = avy-fy[pos]*d/tt ;
}
}
puts("ok");
return 0 ;
}
int main( void )
{
float *vx, *vy ;
unsigned char *image ;
char str[STR_LEN];
int x_size, y_size;
puts("X_SIZE=?");
fgets( str, STR_LEN, stdin );
x_size = atoi( str );
puts("Y_SIZE=?");
fgets( str, STR_LEN, stdin );
y_size = atoi( str );
printf("x_size=%d, y_size=%d\n",x_size, y_size);
vx = (float *)malloc(x_size*y_size*sizeof *vx) ;
vy = (float *)malloc(x_size*y_size*sizeof *vy) ;
image = (unsigned char *)malloc(2*x_size*y_size*sizeof *image) ;
read_data( image, x_size, y_size ) ;
determine_of( image, vx, vy, x_size, y_size ) ;
write_data( vx, vy, x_size, y_size ) ;
26
静止気象衛星画像による風ベクトルの推定
return 0 ;
}
付録 4 のプログラムの流れを解説する。
まず，ファイルに格納されている２枚の衛星画像を C 言語の配列に読み込む。
読み込んだ衛星画像データに対して勾配法を適用する。このとき，勾配法のオイラー・ラグランジュの方
程式から導かれる大規模連立方程式を解くために，ガウス・ザイデル法と呼ばれる反復法を用いているが，
この反復回数をパラメータとして指定しておく。
また，条件式の中のパラメータ α の値も指定する。
これらのパラメータの下で，勾配法によって得られた風ベクトル場の情報をファイルに格納する。
27