2019 年度情報処理学会関西支部 支部大会

G-39

物理量の時空間分布推定と時空間特異点検知
Time and Spatial Distribution Estimation of Physical Parameters and Singularity Detection
武内 俊樹 †

岡留 剛 †

Toshiki Takeuchi

Takeshi Okadome

概要
西日本に起こった平成 30 年 7 月豪雨のような 100 年に 1 度と言われる異常気象（特異点）に対して，事前に予測することは一般に難しい．特異点とは，過去に経験した現象や観測したデータと比べ，大きく外れた現象・データのことを示す．ただし，毎年のように起こる現象はいかなるものでも特異点ではなく，正常なものとして扱われることに注意する．本研究では，気象において過去に前例のない（100 年に 1 度起こると言われるような）特異点の検知を目的とし，物理量の時空間分布推定と時空間特異点検知手法を提案する． 時空間分布推定では，高精度センサと低精度センサを利用し，個々のセンサ値を補正または欠損値を推定する．時空間分布推定により，高精度センサがない地
域においても特異点検知手法を適用できる利点がある．特異点検知手法では，ある地域の降水量時系列データを部分
時系列に変換し，部分時系列における各ベクトルを多次元空間上の 1 点に対応づける．多次元空間上で近傍距離を計算し，一定の基準距離（閾値）と比較することで部分時系列ごとに異常を判定する．予備実験として，降水量において過去に特異点があった 20 箇所の地域に対し，提案手法によって対象の特異点を早く確実に捉えられるか検証した．

For abnormal weather (singularity) said to be once in 100 years like heavy rain in July 2018 which happened towestern Japan, It is generally diﬃcult to predict in advance. Singularity refers to phenomena and data that deviatesignificantly from phenomena experienced in the past and observed data. However, it should be noted that any
phenomena that occur every year are treated as normal ones, not as singular points. In this study, we propose timeand spatial distribution estimation of physical parameters and singularity detection methods for the purpose of
detecting unprecedented singularity in the weather (which is said to occur once in 100 years). In time and spatial
distribution estimation, high accuracy sensors and low accuracy sensors are used to correct individual sensor values
or estimate missing values. By time and spatial distribution estimation, there is an advantage that the singularpoint detection method can be applied even in areas without high precision sensors. In the singularity detection
method, the precipitation time series data of a certain area is converted to a partial time series, and each vector
in the partial time series is associated with one point in the multidimensional space. The neighborhood distanceis calculated in multidimensional space, and the anomaly is judged for each partial time series by comparing the
fixed reference distance (threshold). As a preliminary experiment, we verified whether the proposed method could
capture the singularity of the object quickly and reliably by the proposed method, for the 20 regions that had
singularity in the past in the precipitation amount.

1.

はじめに

西日本に起こった平成 30 年 7 月豪雨のような 100 年に側面がある．気象においても過去に前例のない異常は，従来の特異点検知手法を適用できる可能性がある．

1 度と言われる異常気象（特異点）に対して，事前に予そこで本研究では，気象において過去に前例のない測し判断することは一般に難しい．特異点とは，過去に（100 年に 1 度起こると言われるような）特異点の検知経験した現象や観測したデータと比べ，大きく外れた現を目的とする．物理量は降水量を扱い，特異点は記録的

象・データのことを示す．主に特異点検知手法は，機器豪雨などに着目する．アメダス雨量値（高精度）とレーの故障やビジネスにおいて売り上げ変化を捉えるといっ

ダー雨量値（低精度）を利用する．ただし，毎年のようにた実用的な分野で活用されている．記録的豪雨に代表さ起こる降水パターンはいかなるものでも特異点ではなく，れる異常は過去に前例のない降水量を観測するため，機正常なものとして扱われることに注意する．また，観測器の故障時に機器が過去にはない動作をする傾向と同じする物理量における異常は時間変動として極少数である
ものとする．本研究は，物理量の時空間分布推定と，時†関 西 学 院 大 学 大 学 院, Kwansei Gakuin University
Graduate

空間特異点検知の大きく２段階に分け，手法を提案する．

時空間分布推定では，高精度センサと低精度センサを

過去のデータに前例がないパターンを発見することに着利用し，個々のセンサ値を補正または欠損値を推定する．

目し，気象予報士でも予測困難である特異点を予測する

推定したい空間を特定のメッシュで区切り，メッシュのことを目的としている．

格子点における物理量を，マルコフ確率場を利用して空間分布推定を行なう．実用上，環境物理量には不連続性を考慮しないといけないため，ラインプロセス手法により不連続領域を検出する．時空間分布推定を行なうために，連鎖グラフを導入し確率変数間の因果関係を持たせることで時系列データの扱いを可能とする．レーダー雨量値は日本全国において降水量を観測したものであり，ノイズが大きく欠損値があるため，高精度センサがない地域で特異点検知手法が使用できない問題がある．時空間分布推定を行なうことで，アメダスがない地域に対しても特異点検知手法を適用できる利点がある．
特異点検知手法では，ある地域の降水量時系列データ
を部分時系列に変換し，部分時系列における各ベクトルを多次元空間上の 1 点に対応づける．多次元空間上で近傍距離を計算し，一定の基準距離（閾値）を比較することで部分時系列ごとに異常を判定する．予備実験として，降水量において過去に特異点があった 20 箇所の地域に対して，提案手法によって対象の特異点を早く確実に捉

3.

時空間分布推定

本研究は，物理量の時空間分布推定と時空間特異点検
知の２段階に分けて説明する．まずは時空間分布推定手法を説明し，後に時空間特異点検知手法を説明する．
図 1 に示すように，時空間分布推定では，地上レー
ダー雨量（低精度）と地域気象観測システムから得られるアメダス雨量（高精度）を用いる．地上レーダー雨量
は，日本全国の雨量を観測することができるが，ノイズが大きく，時刻によっては欠損値が存在する．そのため，
アメダスがない地点において特異点検知を行なう場合は，
地上レーダー雨量のノイズを小さくし，欠損値を推定した後が望ましい．まず本研究では，少数の高精度センサと多数の低精度センサを用いてセンサフュージョンを行
ない，物理量の時空間分布推定を行なう．さらに，ラインプロセス手法により局所的な降水量の変化を捉え，不連続性を考慮する．

えられるか検証した．その結果，20 箇所で観測史上１位となる記録的豪雨を 24 時間以上前に検知でき，高い確率で特異点検知が可能であることを示した．しかし，降水量のみでは様々な特性を持つ特異点を捉える精度は低
く，降水量以外の物理量を利用した手法の検討が今後の重要な課題である．本論文の貢献は，レーダー雨量の観測可能な地点ならば，特異点検出手法を気象に適用し，
記録的豪雨といった特異点は検出可能としたことである．

図 1: アメダス雨量（左）と地上レーダー雨量（右）の

以下各節で，関連研究と提案手法・予備実験について詳例．気象庁が提供する気象観測データより (Copyright(c)

細に述べまとめを記す．

2006 by Japan Meteorological Agency)．

2.

関連研究

車両機器に関する異常検知手法として，近藤稔 [1] の手法がある．機器の振動や音の音色をデータ化し，時系

3.1

物理量の不連続性を考慮した空間分布推定

まず時空間分布推定の前に，マルコフ確率場 [3] を利

列データを部分時系列に変換した後，多次元空間上の 1用した空間分布推定を説明する．マルコフ確率場（無向点の座標に対応付ける．対応付けられた点と他の点とのグラフ）により，確率変数間の緩い依存関係を表現した

近傍距離から異常度を算出し，ある閾値と比較すること

例を図 2 に示す．推定したい空間を決め，地上レーダーで，人が判断するよりも早く故障を発見することが可能

雨量の観測点を通るようにメッシュで切る．メッシュのである．本研究では，降水量データを用いて，人よりも格子点に真の値を反映すると思われる i 番目の確率変数

早く特異点を発見することを目的とする．

を xi とする．低精度センサの観測値は図中における水色

また，機械学習と気象衛星ひまわり 8 号の雲画像によのノードに対応し，高精度センサの観測値は赤色のノーる台風検出手法として，金崎拓郎ら [2] の手法がある．過ドに対応する．例えば，図 2 右のモデルにおける中央の去の気象衛星雲画像を大量に用いて台風位置を学習する
ことで，発達前の台風や衰弱した台風を含めて検出できる．しかし，異常気象の中には突発的なものも含まれ，台風以外が原因で発生するものも存在する．本研究では，ノードは隣接する 4 方向のノードと、対応しているセン
サに依存関係がある．また，無向エッジにはそれぞれ対応するエッジ変数があるとする．

は，使用するセンサのダイナミックレンジで決定する．
∑
2
E (xC , uC ) =w0{li − xi }

エッジ変数

確率変数

i∈Dlow+ w1：真の値が反映する確率変数：高精度センサ値∑2{hi − xi }i∈Dhigh：低精度センサ値+ w2

図 2: 兵庫県をメッシュで切った例（左）と 3×3 格子の

∑

∑

(3)

f (xi , a)i∈N a∈bd(xi )

例（右）．

2

f (xi , a) = {xi − a} {1 − σ (ua )}このとき，マルコフ確率場の同時確率は，クリークの概念を導入すると便利に表現できる．クリークは全結合+ cσ (ua )の条件を満たしたグラフの部分集合であり，あと 1 つ+ σ (ua ) {1 − σ (ua )}(4)

ノードを加えると全結合ではなくなる状態のクリークを
極大クリークとする．クリークは C で表し，クリーク内σ=の変数集合を xC とし，クリーク内のエッジ変数集合を

1
1 + e−ua(5)

uC とする．Z は規格化定数である．同時分布 p(x, u) に

関数 f は，ラインプロセスを含んだ関数であり式 (4)対して，極大クリーク上のポテンシャル関数の積で表現で表現される．以上の式を定義し，同時確率 p(x, u) がしたものを以下の式とする．

観測値が与えられたもとで事後確率最大になるよう推定を行なう．実用上は物理量の不連続領域を考慮しなけれ1 ∏p(x, u) =ψC (xC , uC )Z

(1)

ψC (xC , uC ) = exp {−E (xC , uC )}(2)

ばならず，式 (3) の第 3 項がその役割を果たしている．
次節でラインプロセス手法による不連続領域を考慮する

C手法を示す．

ポテンシャル関数はエネルギー関数 E (xC , uC ) によって決まるものとする．このエネルギー関数を最小化することによって，与えられた観測値にフィッティングする [4]．
第 1 項と第 2 項はデータフィッティング項であり，第 1 項は低精度センサ，第 2 項は高精度センサと対応している．
第 3 項は拘束条件項であり，隣接する確率変数間の依存関係を表現している．さらに，拘束条件項にラインプロ
セス手法 [5] を導入し，不連続を考慮可能としている．な
お，ラインプロセス手法の詳細は次節説明する．格子点の総数を N とし，x = {x1 , x2 , ..., xN } とする．Dlow は
低精度センサ値を持つ格子点の集合であり，i 番目の格子点に対応する低精度センサ値を li とする．Dhigh は高精度センサ値を持つ格子点の集合であり，i 番目の格子点に対応する高精度センサ値を hi とする．ただし，高精度センサの数は低精度センサの数に比べ圧倒的少数とする．bd(xi ) は確率変数 xi に隣接する確率変数集合である．このとき，ua は確率変数 xi と確率変数 xi に隣接する確率変数 a を接続するエッジに対応したエッジ変数とする．σ はシグモイド関数であり，式 (5) で表される．

c はラインプロセス手法における閾値であり，不連続領域の検出しやすさを意味する．各項に関わるパラメータ3.2

ラインプロセス手法による不連続領域検出

環境物理量の空間分布推定において，実用的には不連続領域の検出が必要である．不連続領域とは，照度においては人の影の領域，雨量においては浮遊物によって雨が遮断される領域，温度においては建築物の壁などがある．不連続領域における観測は周囲とは全く異なる値を取るため，センサフュージョンの際には精度低下の原因になる．そこでラインプロセス手法により，不連続領域
を別の領域として分割し，領域ごとに空間分布推定を行なう．
ラインプロセス手法は，式 (4) の関数 f に既に反映されている．2 乗誤差は確率変数 xi と xi の隣接する確率
変数 a との依存関係を表しており，xi と a の値が近づくほど，エネルギー関数によって返される値が小さくな
る．式 (4) は閾値 c よりも 2 乗誤差が小さい場合は二乗誤差によって得られる値をエネルギー関数に加え，大きい場合は閾値 c をエネルギー関数に加えることを意味する．なぜなら，エネルギー関数の最小化を行なうため，最小化に有利となる値を選択していくからである．すな
わち，2 乗誤差と閾値 c を比較し，値が大きい方はエネルギー関数を大きくするために無視し，値が小さい方をエネルギー関数に加える．この流れは，ユーザーの手によらず式 (3) のエネルギー関数を最小化する過程で行なわれる．その仕組みを図 3 に示す．

では，このブロックと順序を用いて時系列データを扱う．
また，図 4 中央は有向グラフであり，ブロック b1 = {1}，b2 = {2}，b3 = {3}，b4 = {4}，b5 = {5} のブロック集
合を持つ連鎖グラフと言える．図 4 右は無向グラフであり，b1 = {1, 2, 3, 4, 5} のブロック集合を持つ連鎖グラフと言えるため，連鎖グラフは無向グラフと有向グラフを
内包する．
１

４

１

４

１

４

２

５

２

５

２

５

３

最小化

３

３

図 4: 連鎖グラフの例．

図 3: 関数 f が最小化で起こる動作の様子．横軸 ua ．
前章に紹介した空間分布推定の役割をする無向グラ
フに，t 時刻目の状態を表す潜在変数を z (t) を対応させ
関数 f が最小化で起こる動作として，図 3 の下図から

ua が ∞ か −∞ に徐々に近づいていくことを意味する．
近づく速さが遅い場合は，重み付けなどによる工夫が必要である．ua が ∞ に近づいた場合，関数 f において二乗誤差項が消えるので，エネルギー関数に含まれず依存関係がなくなる．その代わり閾値 c がエネルギー関数に加わるが，閾値は定数であるためエネルギー関数の最小化には何も影響を与えない．これが不連続領域と解釈さる．図 5 に格子を 2 × 2，時系列を 2 とした提案モデル{}の連鎖グラフを示す．図 5 のブロックは，b1 = z (1) ，
{
}
{
}
(1)
(1)
(1)
(1) (1) (1)(1)
b2 = x1 , x2 , x3 , x4 , l1 , l2 , h1 ，b3 = z (2) ，
{
}
(2)
(2)
(2)
(2) (2) (2)(2)
b1 = x 1 , x 2 , x 3 , x 4 , l 1 , l 2 , , h 1となる．潜在変数間に連結される有向エッジは時刻間の因果関係を表し，潜在変数 z と確率変数 x を連結する有向エッジは，状態と真の値と思われる確率変数との因果関係を表す．

れる．逆に，ua が −∞ に近づいた場合，関数 f において閾値 c が消えることになり，二乗誤差項がエネルギー関数に加えれる．二乗誤差項の値が，閾値 c よりも小さい場合は連続領域と解釈される．

3.3

連鎖グラフを用いた時空間分布推定

これまで解説した空間分布推定から，時系列データを扱う時空間分布推定を考える．無向グラフは確率変数間の緩やかな依存関係を表し，有向グラフは確率変数間の図 5: 2 × 2 格子，時系列の長さ 2 の連鎖グラフ．

因果関係を表す．時系列を扱うためには，時刻間の因果関係を表現する有向グラフが必要であるため，連鎖グラ

低精度センサ値集合 l ，高精度センサ値集合 h ，学習

フ [6] [7] を導入する．連鎖グラフは，無向グラフと有向済みパラメータ集合 θ が与えられたもとで，真値が反映グラフを融合したグラフであり，有向エッジと無向エッされる確率変数集合 x ，潜在変数集合 z，エッジ変数集

ジを両方持つグラフに拡張できる．

合 u の推定を行なう．θ = {π, A, ϕ} であり，遷移確率

連鎖グラフを扱う上で，新たに「ブロック」と「順序」

A，出力確率 ϕ，初期要素は π とする．潜在変数の状態の 2 つの概念を導入する．ブロックとは，無向グラフのみ

数は訓練データを分割し，最も精度の良い状態数を与えで形成される部分グラフのことであり，ブロック内のグるものとする．無向エッジのみで接続されたノード集合

ラフは有向エッジを持たないとする．例えば図４の左図であるブロックの概念を導入し，グラフのブロック数をにおいて，ブロック b1 = {1, 2}，b2 = {3}，b3 = {4, 5}

B としたとき，連鎖グラフの同時分布は，ブロックの積の 3 つのブロックに分けられる．ブロックの概念を導入の形として以下の式で因数分解の表現ができる．

したことにより，ブロック同士は有向エッジのみで連結され，ブロック間に順序の概念が生まれる．連鎖グラフ

p(x, z, u|l, h, θ) = p(b1 )∏B

i=1 p(bi |bi−1 )(6)ブロック単体は無向グラフであるため，ブロックごとにマ分時系列によって時系列データを任意の部分時系列デー

ルコフ確率場の枠組みを与える．ブロックごとに同時分タに変換し，多次元空間上に 1 点に対応付ける．対応付布を計算することで確率的解釈が可能となり，p(x, z, u)

けられた点に対して，近傍距離を利用し異常判定するたの同時分布が最大になるように推定を行う．なお図 5 にめに参考する異常度を定義する．異常を判定する閾値とおいて，大域的マルコフ性 [7] により，連鎖グラフにおして基準距離を決定し，基準距離より大きい場合は対応

ける条件付き独立性 b2 ⊥⊥ b3 | b1 があるためブロックの積

する部分時系列データは異常であると判定する．最後に，で表現できる．例として，図 5 の 2 時刻における連鎖グ

気象庁が提供する降水量データを用いて，過去の特異点

ラフの因数分解を考えると，以下の式になる．

を提案手法で検出可能か予備実験を行なった．

p(b1 , b2 , b3 , b4 ) = p(b1 )p(b2 | b1 )p(b3 | b1 )P (b4 | b3 ) (7)p(b1 ) = p(z

(1)

(1)

(1)

(1)

(1)

(2)

(2)

(2)

時系列データを部分時系列に変換

降水量の時系列データの中から異常部位を検出するため，１つの時系列データをスライド窓によって複数の部

)

(8)
(1)

(1)

(1)

p(b2 | b1 ) = p(x1 , x2 , x3 , x4 , l1 , l2 , h1 | z (1) )
(9)
(2)
(1)
p(b3 | b1 ) = p(z | z )
(10)
(2)

4.1

(2)

(2)

(2)

p(b4 | b3 ) = p(x1 , x2 , x3 , x4 , l1 , l2 , h1 | z (2) )(11)
このとき，ある時刻 t における p(bi+1 | bi ) は前節の説

分時系列に表現する [8] [9]．スライド窓の窓幅を 3 時間に設定した場合，3 時間ごとの部分時系列に分割される．
例として，i 番目の地点における長さ T の時系列データを以下の式で表現する．
{
}
(1)
(2)
(T )
D = xi , x i , . . . , x i

式 (13) を M 次元の部分時系列に変換する際，式 (13) の時系列データは以下の式で変換される．
{
}
(1)
(2)
(N )
D = mi , mi , . . . , mi明と同様，クリークの概念を用いて因数分解を行なう．

N =T −M +1

なお，連鎖グラフのモラルグラフとして，有向エッジは
(t)

1
(t)
Zi

(t)

(t)

(t)

(15)
(t)

本研究で扱う降水量の時系列データを部分時系列に変換ψ1 (x1 , x2 )
(t)

(14)

このとき，mi は M 次元ベクトルであり，mi
=
{ 1 2
}
{ 2 3
}
(t+1)
M +1
M
xi , xi , ..., xi ，mi= xi , xi , ..., xiとなる．

無向エッジへと変換している．

p(bi+1 | bi ) =

(13)

(t)

する様子を図 6 に示す．図中のグラフは，岡山県倉敷市(t)

ψ1 (x1 , x3 )ψ1 (x2 , x4 )(t)(t)(t) (t)ψ1 (x3 , x4 )ψ1 (x1 , l1 )(t)(t)(t) (t)ψ1 (x1 , h1 )ψ1 (x4 , l2 )(t)(t)
ψ2 (z (t) , x1 )ψ2 (z (t) , x2 )(t)(t)
ψ2 (z (t) , x3 )ψ2 (z (t) , x4 )の 2018 年 1 月 1 日 1 時から 2018 年 12 月 31 日 24 時までにおける 1 時間ごとの降水量の推移を緑線で表してい(12)る．なお，部分時系列は 1 時刻ずつスライドして分割していく．
部分時系列ごとに，異常度（大きいほど異常であることを示す）を求め，異常部位がどこか判定する．部分時
系列ごとに異常度を計算することで，継続した異常を検

ポテンシャル関数 ψ1 を決定するエネルギー関数は前章

知でき，観測した際のノイズによる誤検知を減らすことで設定した式 (3) を利用する．ただし，状態の潜在変数ができる．異常度の計算と，異常部位の判定方法は次のを含むポテンシャル関数 ψ2 はガウス分布の式を利用す

章で詳細を記す．窓幅が小さいほど集中豪雨などの短期

る．レーダー雨量値における欠損値などを推定すること的な異常気象を捉えることに便利であり，窓幅が大きいが可能となり，アメダスがない地域において時空間特異と梅雨前線による多雨など長期的な異常気象を捉えるこ

点検知が可能となる．

とに便利である．

4.2

4.

時空間特異点検知

近傍法による異常度の決定

部分時系列の変換によって得られたデータ D を訓練

時系列データを任意の長さで分割した部分時系列ごとデータとし，新たに与えられたデータが異常か判定する．

に異常を判定する．過去の前例のない異常を捉えることこのとき部分時系列の各要素を，多次元空間上の 1 点のが目的であり，例年のように頻発している事象は特異点座標に対応付ける．2018 年 6 月と 7 月における倉敷市のには含めないものとする．手法の流れとしては，まず部

降水量の推移を，スライド窓の窓幅は 2 として，多次元2018年7月倉敷市の降水量

・・・

(mm)

降水量 (mm)

・・・

降水量

降水量の
時系列データ2018年6月倉敷市の降水量

時間

部分時系列

30mm

図

6:
時系列データを部分時系列に変換する
図 ．ア メ ダ ス 雨 量 値 は 気 象 庁 ホ ー ム ペ ー ジ（https://www.data.jma.go.jp/gmd/risk/obsdl/index.php）
より参照．

空間に対応づける様子を図に示す．窓幅は 2 であるため，

2 次元空間に対応付けられる．図に示した 2 次元空間の横軸は，2 時刻で分割されるうちの 1 つ目の時刻におけ

30mm

0mm

図 7: 部分時系列を多次元空間上に対応づける様子．

る値を表しており，縦軸は 2 つ目の時刻における値を表
す．青で示した 2 時刻は，多次元空間に対応づけると，の近傍を考慮した方が精度が向上する．本研究では，訓１箇所に固まることがわかる．赤で示した 2018 年 7 月

練データを複数に分割し，最も精度が良かった近傍デー

豪雨における 2 時刻を多次元に対応づけると，他のデー

タ数 kN N に決定する．

タと離れた位置にプロットされる．ただし，訓練データ

異常度の式 (16) で，図 7 における倉敷市の降水量デー

は正常標本がほとんどであり，異常標本は限りなく少なタを，異常度に変換した様子を図 8 に示す．緑線は降水

いものとする．

量の推移を表し，青線は異常度の推移を表している．異常度は 0 から 1 の間に収まるよう正規化しており，大き異常度を k 近傍法によって，他のデータ間とのユーク

いほど異常であることを意味する．スライド窓の窓幅をリッド距離で決定する．新しく与えられるデータベクト

10 時間と設定している．スライド窓の幅を調節すること

ル m の異常度を，a(x) としたとき，異常度は以下の式

で，捉えたい特異点の性質が異なる．本研究では，複数で近藤稔 [1] より定義されている．

の窓幅を用いて，別々に異常度を計算する．

a(x) =

1
kN N

k∑
NN
k=1

|x − N Nk (x)|
−1
dk

(16)

異常度

降水量

降水量

番目に近いデータを表す．dk は k 番目に近いデータに対異常度

kN N は近傍データ数を表し，N Nk は入力データ m と k

2018年7月倉敷市

2018年6月倉敷市

応する基準距離（閾値）であり，基準距離は k 番目ごとに別々に値を求める必要がある．入力データ m と他のデータ点との距離を基準距離で割り 1 を引くと，負のときは入力データ m は正常であると判断し，正のときは入図 8:

力データ m は異常であることを表現する．異常度 a(x)ス ラ イ ド 窓 の 窓 幅 を 10 時 間 と し て 異 常 度 を 計 算

は，小さい値であるほど正常であり，大きい値であるほしている．アメダス雨量値は気象庁ホームページ

ど異常である．1 近傍のみだと訓練データに偶然似た異（https://www.data.jma.go.jp/gmd/risk/obsdl/index.php）

常度があった場合，入力データとの距離が近く正常と判より参照．

断される可能性がある．また時刻が隣り合うデータが両方とも異常であった場合，先の時刻は異常と判断されるが後の時刻は正常と判断される可能性があるため，複数

図 7 の降水量データを異常度に変換した図．

4.3

基準距離の決定では報告レポートに記載がある最大 3 時間降水量から最新しい入力データが異常か正常かの判断のためには，

大 72 時間降水量まで扱う．気象庁の報告レポートに記基準距離を決定する必要がある．基準距離を大きくする載があった 20 箇所において特異点検出可能か，また最と，正常標本を間違えて異常であると検出しないように

大降水量に至るどのくらい前に検出できるか検証した．

なるが，逆に異常を見逃す問題がある．また，基準距離を

以下表 1 に，今回の実験の対象となった特異点を記す．

小さくすると，確実に異常を発見するが，逆に正常標本を間違えて異常と判断する問題がある．そこで基準距離

表 1: 実験対象の特異点．出典：気象庁ホームページの決定には，この 2 つのトレードオフを考えるため，正（https://www.data.jma.go.jp/gmd/risk/obsdl/index.php）

常標本精度 r0 と異常標本精度 r1 を用いる．井出剛ら [8]をもとに作成．地点名とは，アメダス雨量値を観測可能は以下のように説明している．正常標本精度とは，実際である箇所を示す．

に正常である標本の中で正しく正常と判定できた数を実際に正常である標本の総数で割った値である．異常標本

特異点

都道府県

市町村

地点名

降水量の時間幅

No.1

高知県

宿毛市

宿毛

6時間降水量

No.2

広島県

呉市

倉橋

12時間降水量

定できた数を実際に異常である標本の総数で割った値で

No.3

岡山県

笠岡市

笠岡

24時間降水量ある．異常度の基準距離により，正常標本精度と異常標

No.4

兵庫県

養父市

八鹿

48時間降水量

No.5

愛媛県

今治市

今治

72時間降水量

岐阜県

大野郡白川村御母衣

72時間降水量

No.7

山口県

岩国市

岩国

72時間降水量

No.8

大阪府

豊中市

豊中

72時間降水量No.9鳥取県

八頭群若桜町若桜

72時間降水量

No.10

福岡県

田川郡添田町添田

72時間降水量

No.11

茨城県

古河市

古河

3時間降水量

No.12

宮城県

栗原市

鶯沢

24時間降水量No.13 平成27年台風18号 埼玉県

越谷市

越谷

24時間降水量

No.14

栃木県

日光市

今市

24時間降水量

No.15

福島県

南会津群南会津町 舘岩

24時間降水量

No.16

三重県

津市

笠取山

3時間降水量

No.17

北海道

下川町

下川

3時間降水量No.18 平成26年8月豪雨 京都府

福知山市

福知山

24時間降水量

No.19

石川県

羽咋市

羽咋

24時間降水量

No.20

徳島県

阿南市

蒲生田

48時間降水量

精度とは，実際に異常である標本の中で正しく異常と判

本精度は大きく変化するため，正常標本精度 r0 と異常
標本精度 r1 を基準距離の関数とする．
正常標本精度と異常標本精度で異常の見逃しの指標としている．正常標本精度 r0 = 1 のとき，正常標本に対しては精度が良く判定できるが，異常標本を発見できないことを意味する．異常標本精度 r1 = 1 のとき，異常標本に対しては精度が良く判定できるが，正常標本を異常と判定することになる．正常標本精度と異常標本精度の 2 つの観点から見て，丁度良い箇所を分岐点とし，基準距離に設定する．実際には，正常標本精度と異常標本精度をグラフにプロットし，交差する点を調べる必要がある．

2r0 r1(17)
r0 + r1
実用上は式 (17) のように F 値を求めると便利である．

No.6

平成30年7月豪雨

f≡

異なる基準距離で繰り返し計算し，F 値が最大となる基準距離を探すことになる．近傍距離より異常度を定義しているため，各近傍ごとに基準距離を求める．

5.

予備実験

提案した時空間特異点検知手法を用いて，特異点を検出できるか実験を行なった．気象庁が提供するアメダス

2008 年から 2018 年の 10 年間の降水量（1 時間ごと）
を利用し，降水量の時系列データを部分時系列に変換し多次元空間上の 1 点に対応づける．式（16）より異常度を計算し，異常度が基準距離より大きい場合は特異点と判断する．スライド窓の窓幅は，3 時間，6 時間，12 時間，24 時間，48 時間，72 時間の窓幅を設定し実験を行
なった．

5.1

実験結果

雨量値を利用する．対象とする特異点は，気象庁もしく実験結果を表 2 に示す．表 2 の No.1 の記載は，表 1 のは地方自治体の報告レポートなどに，異常気象であり，No.1 と対応している．最大降水量の観測時刻とは，対応

過去に事例がない事象かどうかはっきりと明記されていする地点において観測史上 1 位となる降水量を観測したるものとする．主に気象庁ホームページで公開されてい時刻である．なお一定期間の降水量の総量であり，表 1

る「異常気象の特徴と要因に関する情報」の報告レポーの降水量の時間幅を参考にしている．特異点の検知時刻

トを参照した．報告レポートで最大降水量の観測時刻のとは，特異点検知手法によって検知した時刻である．表

記載がある地点に着目した．またある一定期間の雨量に2 に示した時刻はスライド窓の窓幅を考慮した時刻であおける総量で，降水量の値を見ていくものとする．一定

り，最も早く検知可能な時刻である．結果として，20 箇

期間を 3 時間とした場合，3 時間の降水量の総量が観測所のうち全ての箇所において異常検知可能であった．ま

史上最大であることを最大 3 時間降水量という．本実験

た，最大降水量の観測時刻より，特異点検知した時刻の方が 3 日ほど早い結果となった．異常気象で大きな被害における時空間分布推定と時空間特異点検知を併用し，をもたらすほどの過去に前例のない豪雨が来る 3 日ほど

レーダー雨量で観測可能な場所は多種多様な特異点検知

前には，検知可能であることがわかる．

が可能になることを目指す．

表 2:

実 験 結 果 を 示 す．最 大 降 水 量 の 観 測 時 刻 と

参考文献

は ，対 応 す る 地 点 に お い て 観 測 史 上 1 位 と な る[1] 近藤稔. 振動のオクターブバンド分析を用いた異常

降 水 量 を 観 測 し た 時 刻 で あ る ．特 異 点 の 検 知 時 刻検知法による車両機器の診断. 日本機械学会論文,と は ，特 異 点 検 知 手 法 に よって 検 知 し た 時 刻 で あVol. 84, No. 862, 2018.

る．最大降水量の観測時刻より，特異点検知した時刻の方が早いことがわかる．出典：気象庁ホームページ[2] 金崎拓郎, 筆保弘徳, 加瀬紘熙, 松岡大祐, 吉田龍二.

（https://www.data.jma.go.jp/gmd/risk/obsdl/index.php）
をもとに作成．
最大降水量の観測時刻

特異点の検知時刻

機械学習を用いた台風検出器の開発. 電子情報通信学会信学技報, Vol. 118, No. 197, AI2018-15, pp. 13-18,

2018.

月日

時分

月日

時分

No.1

2018/7/8

8:20

2018/7/5

7:00

No.2

2018/7/7

5:30

2018/7/4

7:00

No.3

2018/7/7

8:40

2018/7/3

13:00[4] Ito, K. and H. Hontani (2006). Analysis of cal-No.4

2018/7/7

9:50

2018/7/4

10:00

No.5

2018/7/8

8:00

2018/7/5

7:00

ibration performance of networked sensors using

No.6

2018/7/8

3:40

2018/7/5

3:00

No.7

2018/7/8

8:50

2018/7/2

10:00[3] Bishop, C. M Pattern Recognition and Machine
Learning, Springer (2006).

measurements graphs. SICE-ICASE International
Joint Conference, 1320-1325.

No.8

2018/7/8

0:10

2018/7/5

0:00

No.9

2018/7/8

6:20

2018/7/2

5:00

No.10

2018/7/8

7:20

2018/7/2

6:00

tributions, and the Bayesian restoration of images.
Readings in computer vision, 564-584.

[5] Geman, S. (1987). Stochastic relaxation, Gibbs dis-No.11

2015/9/9

17:20

2015/9/6

14:00

No.12

2015/9/11

8:50

2015/9/8

7:00

No.13

2015/9/10

4:50

2015/9/7

7:00

No.14

2015/9/10

6:20

2015/9/6

13:00

with an extended chain graph for human activity
recognition and image segmentation. IEEE Trans-

No.15

2015/9/10

6:40

2015/9/7

5:00

No.16

2014/8/9

15:50

2014/8/6

17:00

No.17

2014/8/5

10:50

2014/8/2

13:00

No.18

2014/8/17

5:50

2014/8/13

22:00

No.19

2014/8/17

5:50

2014/8/14

6:00

No.20

2014/8/4

4:50

2014/7/31

5:00[6] Zhang, L. (2011).

Probabilistic image modelingactions on Image Processing, 2401-2413.
[7] 宮川雅巳. グラフィカルモデリング. 朝倉書店（2011）[8] 井出剛, 杉山翔. 機械学習プロフェッショナルシリー
ズ 異常検知と変化検知. 講談社（2015）

気象庁から得られる降水量データは１時間ごとである
ため，特異点検知時刻も１時間ごとの単位となる．問題点として，降水量を基準とした異常検知手法であるため，
記録的豪雨などの異常気象しか捉えられない．異常気象には気温など関係するものもあり，降水量の他の環境物
理量を利用して特異点を検知する必要がある．

6.

まとめ

本稿では，物理量の時空間分布推定と時空間特異点検知を提案した．予備実験において，提案手法を用いて特異点検知手法を適用した結果，気象においても特異点は検知可能であることがわかった．また 20 地点において，最大降水量を観測する前に検知することが可能である．
ただし，降水量データのみに着目しており，温度や湿度など他の物理量を活用し，記録的豪雨といった特異点とはまた異なる特異点を捉えていく必要がある．本研究[9] Shieh, J and E. Keogh (2008) . iSAX: Indexing andMining Terabyte Sized Time Series. Proceedingsof the 14th ACM SIGKDD international conferenceon Knowledge discovery and data mining, KDD 08,
623-631.