地磁気観測所テクニカルレポート 第14巻 21-33頁 平成29年
Technical Report of the Kakioka Magnetic Observatory Vol.14, pp.21-33, 2017ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法
笹岡雅宏
地磁気観測所技術課
2016年10月27日受領，2017年 3 月23日改訂，2017年 3 月28日受理

要旨

地磁気脈動 Pi の諸元（発現時刻，周波数，及び振幅）を読取る手作業に代わることを目的として，ウェーブレット解析及びファジィ・ロジックを用いた自動読取り手法を開発した．2012年の女満別データを用いて Pi の諸元を自動で読取り，手動読取り記録と比較した．両者は，包括的に非常に良く一致した．一方で，手動読取り記録においては，振幅には数秒以下の地磁気短周期変動によるバイアスがあること，また，連続発生する一連の Pi の開始時刻は Pi の発現時刻と必ずしも一致しないことが確認された．

1 ．はじめに地上で観測される地磁気脈動は，ULF 波として知られ，磁化プラズマが磁力線と一緒に動いて生じる磁場の歪みが伝播する Alfvén 波や，磁化プラズ
マ密度の圧力変動である磁気音波などの磁気流体波に起因し，主に連続型の脈動現象（Pc: Continuous
Pulsation）と 不 規 則 型 の 脈 動 現 象（Pi: IrregularPulsation）に分類される．一方，Pi は，通常連続して発生し，Bay（湾型変化）もしくはサブストームに伴って出現することがよく知られており，特にオーロラ・サブストームのオンセットを示唆するとみなされている．
地磁気観測所では，1999年以来観測された地磁気脈動（Pi1及び Pi2）の諸元（発現時刻，周波数，振幅等），及び現象のクラスについてウェブ上で公開 し て い る（URL: http://www.kakioka-jma.go.jp/
obsdata/Geomagnetic_Events/Events_index.php）．
2011年以降は女満別観測点のみ地磁気脈動の諸元の読取りを継続している．図 1 に手作業で求める場合の Pi の周期及び振幅（最大変化量）の概念図を描く．
また別途，Pi の顕著さ（明瞭度）を表す Quality としてクラスをその最大変化量から判定してきた．これらの諸元は生データ（0.1秒値）から周期150sec 以
下の成分をハイパスフィルターにより抽出し，そのプロット図とスケールを用いて数値を直接読取る手作業により求められてきたが，手作業の省力化や脈©2017 Kakioka Magnetic Observatory, Japan Meteorological Agency

動現象の即時的な把握のために計算機による自動読取りの実現が望まれてきた．一方で，この煩雑な手動読取り作業がどのようにすれば計算機に移植可能となるのかについては，これまで具体的に検討されることはなかった．
本稿では，Pi1及び出現頻度の高い Pi2を対象に，

図1

Pi の周期（T）及び最大変化量（Â）の読取りの概念
図．図中にローカル・ピークの時刻と振幅を（tj, aj）
で示す．

22

笹岡雅宏

計算機による自動読取りのために開発した手法の数学的ツールと諸元の検出アルゴリズムについて述べる．そして過去データを用いて実施した自動読取り
結果と観測者による読取り記録を比較・検証する．
2 ．方法
信号のスペクトル解析を行う目的でフーリエ変換を利用すると，信号の時間の情報は失われる．スペクトルの時間変化を調べるために窓関数を導入してその区間の信号のスペクトルの時間変化を逐次的に求める短時間フーリエ変換は，窓関数にガウス窓を用いる場合は Gabor 変換と呼ばれる．窓関数の幅を大きく取ると周波数分解能は向上するが時間分解能は低下し，逆に窓関数の幅を小さく取ると周波数分解能は低下するが時間分解能は向上する．時間幅と帯域幅は同時に小さくはできない（時間周波数
解析の不確定性原理）．ウェーブレット変換による
時間周波数解析は，時間の情報を保ちつつ周波数解析を可能にする（例えば，Kaustky and Turcajova,
1994; 新井，2006）．即ち，脈動現象 Pi の発生時間とその時間の変動周期の情報を同時に把握できる．
2 ．1 連続ウェーブレット変換
地磁気脈動 Pi の諸元の自動読取りについて，求められる要件は二つある．一つは，① Pi 出現の自動検出機能であり，もう一つは②諸元の読取り機能である．Pi 及び急始型磁気嵐（SSC）の自動検出については先行研究がある（例えば，Nose et. al.,
1998; Takano et al., 1999）．これら先行研究は，離散ウェーブレット変換を用いて地磁気現象の特徴的変化を抽出することにより地磁気現象が検出できることを示した．一方，諸元の読取りにはデータに対するモデル（関数）のフィッティングが有効である．
連続ウェーブレット変換を用いた時間周波数解析により，関数フィッティングと同様な効果を得ることが可能である．本稿では連続ウェーブレット変換を基礎において，Pi の諸元を求めるアルゴリズムを構築する．
Pi の典型的な形態的特徴は，インパルス的に立ち上がりピークに達した後ダンプするが，出現の間は数十秒以上の周期の振動を伴うというものである．先ず Pi の立ち上がりからピークまでをモデル化し，そのモデル関数を用いて時間周波数解析を行う．一般的にインパルスはガウス関数で上手く近似されるため，次のガウス関数（gθ(t)，式［2］）と三角関数との積（Ψ(t)，式［1］）を連続ウェーブレット
変換に用いるマザーウェーブレット（核関数を構成するための基本参照波）とする．

[1]

[2]
ここに，ω0＝2π/T0（T0＝70sec），σ＝30sec である．
式［2］中の θ がゼロの場合，式［1］のマザーウェー
ブレットは Gabor 関数である．Gabor 関数を用いたフィルターは，一般的に線分を検出する画像処理として知られている．本稿では，パラメータ θ の導入により Gabor 関数の対称性が緩み，Pi の変動との相似性が増大することを期待している．連続
ウェーブレット変換は，時間軸に対して独立に式［1］の関数を伸縮させて時間周波数解析を行う．Gabor関数をマザーウェーブレットに用いる Gabor ウェー
ブレット変換と同様に，式［1］の連続ウェーブレッ
ト変換は時間周波数解析には適している一方，直交基底を持たないため離散ウェーブレット変換には用いられない．式［3］は式［1］がウェーブレット関数であるための条件であり，式［4］はこれを満たすよ
うに式［1］を補正したものである．
[3][4]
連続データ f に対する連続ウェーブレット変換はt→(t−b)/a と置き換えることにより式［5］により与えられる．Ψ(t) は Ψ(t) の複素共役である．パラメー
タ a はダイレーション（時間軸方向に拡大縮小するためのもの），パラメータ b はシフト（時間軸方向に平行移動するためのもの）を示す．Pi の変動と連続ウェーブレットとの相似性を最大化するように最適なパラメータを探す．即ち，式［2］中の θ に数値を選択した数だけの各ウェーブレット関数を用いた
時間周波数解析から，ウェーブレット変換 W が最大となるパラメータ（a, b, θ）を検出する．式［6］は，非定常の変動に対する相似性を評価するためのドーターウェーブレット（マザーウェーブレットを時間軸方向に拡大縮小，平行移動して生成されたもの）
である．
[5]

ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法[6]
ここで，式［5］を離散化することにより計算量を軽減する．ガウス関数の最大振幅に対する減衰比が許容誤差 ε（＝0.001）以下になる条件［7］により，積分区間を有限長にする．M は，Ψ の振幅が小さ
く積分値への影響が十分小さくなる時間幅を示す．
式［7］において，形式的に θ → θ/a と置き換える
ことにより式［8］及び［9］を導出している．Msup とMinf はそれぞれ積分区間の上限と下限を示す．

図 2 は，式［6］のドーターウェーブレットの例示であり，ダイレーション a，シフト b 及び θ の相違について示す．図 2 a は周期40sec の Gabor ウェーブ
レットの例であり，図 2 b はその対称性が緩んだ
例を示す．図 2 c は図 2 a の Gabor ウェーブレット
を半周期分シフトした例であり，図 2 d は図 2 c と逆位相の関係がある．図 2 e は周期70sec の Gabor
ウェーブレットの例であり，図 2 f は図 2 e を半周期分シフトした例である．図 2 g は図 2 e の Gabor
ウェーブレットを 1 周期分シフトした例であり，図2 h は図 2 f と逆位相の関係がある．図 2 g は図 2 e
より図 2 h と相関が強い．
2 ．2

[7]

[8]
[9]
式［10］は式［5］を離散化した近似を示す．無次元量のダイレーション a は周期に相当し，1/a は周波数に相当する．但し，時間周波数解析において a を
物理量的に周期 T（sec）として扱いたい場合には，
a→aω0/2π と置き換えればよい．

[10]

図2

23

離散ウェーブレット変換による多重解像解析
本稿の Pi 検出では，地磁気データ 1 秒値から，
Pi を含む信号成分を抽出した解析データを用いる．
Pi は定常波ではなく時間とともに変化する非定常波であるため，フーリエ変換による帯域フィルターを用いると元の Pi の形態的特長が失われることが起こりうる．このため，コンパクトサポートである
Daubechies ウェーブレットの離散ウェーブレット
変換を用いた多重解像解析によりデータを分解する（例えば，笹岡ほか，2011）．ここで，Daubechies
ウェーブレットのサポート長は16とした．また， 1
時間分のデータ毎に時間周波数解析が適用できるよ
うにデータを分解した．これら分解成分から Pi 以
外の変動成分を区別する手続きは以下の通りである．各分解成分に対して Gabor ウェーブレット変換を適用し，最大の変換係数が求められたときの周期が Pi の周期範囲150sec 以下である分解成分を選択する．また，周期数十秒の Pi の諸元の見積もり

連続ウェーブレット関数のパラメータ例．

24

笹岡雅宏

には，数秒以下の短周期変動は不要であるので，これらも除去する．こうして周期数十秒以上の Pi の形態をクリアに抽出した上で，2.1節の時間周波数解析を行う．
2 ．3 周期及び最大変化量の決定アルゴリズム
2.1節で述べた連続ウェーブレット変換 W が t＝b で最大となるパラメータ a 及び θ を用いると，式［6］の連続ウェーブレット関数における各振動ピー
クの時刻が分かる（a は，Pi の周期を含む10sec 〜
150sec の範囲で，θ は，−T/2−1sec 〜 T/2＋1secの範囲で W(b) を求める）．各振動ピークの時刻に対応する地磁気変動成分（合成波 H5＋H6）から最大変化量を求める．但し，実際の読取り作業では，地磁気変動のピークから周期及び最大変化量を算出しているので，これに合わせるために地磁気変動のピークを決定する．先ず連続ウェーブレット
関数の各振動ピークの時刻を中心とする±10sec の範囲から各ピークを検索する（暗に10sec より大きな周期の Pi を前提にこの範囲の幅を設定している．
± T/4の範囲からとした方が合理的かもしれない
が，周期数十秒以上の Pi を対象にしていること，
また経験的にウェーブレット関数と地磁気変動の各ピークの時刻の違いは5sec 程度であったので，このように設定した）．次に 3 点のピーク値とその時刻を用いて，周期及び最大変化量を求める．
2 ．4 Pi の発現時刻の決定アルゴリズム
Pi の発現時刻 t0は連続ウェーブレット関数のパラメータ（a 及び T）を用いて以下のように決定する（式［11］）．但し，整数 N は閾値 h に関する条件式［12］を満たす．図 3 に算出した発現時刻を時間軸に直交する赤線で例示する．

[11]

[12]

ここに，h＝0.0008である．
2 ．5 ファジィ・ロジックを用いた Pi 検出
観測者は，Pi の形態の概念的特徴とよく合致した信号を見分けることで Pi を検出する．即ち Pi の検
出は観測者の判断という曖昧な基準に基づくため，
本稿ではファジィ・ロジックを用いて Pi を検出す
る．先行研究に自然現象の特定においてファジィ・
ロジックを利用した例がある（例えば，Cornman etal., 1998）．ウェーブレット係数（式［5］）について a及び θ を変数として求めた最大係数 W(b) の大きさに基づくメンバーシップ関数（MA）を以下に定義す
る（式［13］）．

[13]

ここに，μA＝{0.1（X（地磁気南北成分）及び Y（東西成分）成分），0.01（Z（鉛直）成分）}，τA＝0.01である．
次に，Pi の対称性に関する θ の大きさに基づくメ
ンバーシップ関数（MB）を以下に定義する（式［14］）[14]

ここに，μB＝1sec，τB＝0.01sec である．我々が知りたい Pi の発生時間においては，W(b) の数値が大きく（地磁気変動とウェーブレット関数に強い相似性があり），且つ θ が半周期以下である（用いたウェー
ブレット関数の対称性が緩むか，もしくは逆位相を示す）．式［13］及び［14］を用いて Pi の発生時間についてのメンバーシップ関数 MA∩B を以下のようにして求める（式［15］）．
[15]
図3

Pi の発現時刻の算出例．図中に赤縦線で発現時刻
t0= －88sec（N=3）を示す．

また Pi は，W(b) のローカル・ピークが見つかる時
刻に検出されるはずなので，式［15］の MA∩B で示さ

ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法れる属性の強い時間に含まれる W(b) のローカル・
ピーク時刻 bj を検出すればよい．ここでローカル・
ピーク W(bj) は以下の条件を満たす（式［16］）．
[16]
ここに，δt＝1sec である．具体的な Pi 検出の手続
きについては，MA∩B の閾値を0.5として，それ以上の数値を示す W(b) のローカル・ピーク時刻 bj を抽
出する．

図4

25

3 ．結果
Solar Cycle 24が始まって太陽活動が比較的活発であった2012年の女満別観測点の Pi の諸元について調査した．用いた地磁気データは，生データ0.1
秒値ではなく毎秒値である．図 4 は，地磁気毎秒値
Y 成分の多重解像解析による分解例（H1〜 H8及び
L0）を示す．この各分解成分に対して Gabor ウェー
ブレット変換を適用した結果，H1から H4については各最大変換係数が H5以降と比べ一桁以上小さく，
且つ H7，H8及び L0は周期が150sec より大きいことから，H5＋H6が Pi を検出する対象成分となった．

Daubechies 離散ウェーブレット多重解像分解を用いた地磁気データ分解例．女満別観測点の地磁気毎秒値 Y（東西）成分：2012年 4 月 5 日08:00UT ～09:00UT を示す．

26

笹岡雅宏

図 5 a は Pi の諸元を検出する対象データであり，
図 5 b はそのウェーブレット係数である．図 5 b の黒線は比較のための Gabor ウェーブレットを用いた変換係数であり，赤線は式［4］を用いた結果である．図 5 b で示される比較は，Pi の非対称性を考慮した方が Pi とウェーブレット関数の相似性は高まることを示す．図6a は，W(b) のローカル・ピーク
を示し．図 6 b はメンバーシップ関数 MA∩B を示す．
図 6 c は，図 6 a のローカル・ピークに対応する
図 6 b の MA∩B のプロットである．MA∩B が0.5以上の W(b) のローカル・ピーク時刻 bj が Pi の検出時
刻として抽出されている．
自動読取りで検出された Pi の諸元リストを表 1に示す．表 1 の Pi について，図 7 に点線囲みで示
す．表 2 は観測者による Pi の手動読取り記録を示
し，図 7 に矢印で示す．Pi の自動読取り手法は，
手動読取りより多くそして網羅的に Pi の特徴が見られる信号を検出しており，ほぼ Pi の形態のイメー
ジ通りの検出結果である（検出数を少なくすることも可能であるが，あえて網羅的に Pi が検出可能と

図5なるようにメンバーシップ関数の閾値 μA はかなり下げて設定している）．図 8 a 及び 8 b は，それぞれPi 検出時刻が08:30:21UT 及び08:58:20UT におけるウェーブレット関数であり，図 7 中の対応する Pi
波形とよく似ている．手動読取りにおいては，発現
時刻は分単位までしか読んでおらず，また周期はPi2の場合 5 秒刻みでしか読んでいない．手動読取りの確度の限界を考慮しての措置かもしれないが，計算機処理の場合は秒単位まで対応可能なので，両者を比較する際には手動読取りの Pi 発現時刻の秒
値は00秒とした．現象の Quality のクラスについて
は，通常時と擾乱時の二つの基準（表 3 参照）があり，各基準ごとに求めた．用いた基準の記録を残していないので，手動読取り記録においてはどちらの基準が採用されたのか不明である．表 1 と表 2 の比較では，発現時刻は最大53秒の食い違いがあり，表 2 のクラス A 及び B の Pi については，表 1 では 1 クラ
ス小さく異なる．最大変化量については，手動読取
り記録が自動読取り結果よりも大きな数値となった．その主な理由は，ハイパスフィルター処理によ連続ウェーブレット時間周波数解析の結果．
（a）Pi 読取り対象データ，（b）連続ウェーブレッ
ト解析の結果．図 4 の時刻に対応する．図 5 b 中の黒線は参考値として Gabor ウェーブレッ
トを用いた結果，赤線は第2.1節中式［ 4 ］を用いた結果をそれぞれ示す．

ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法

るデータには数秒以下の短周期変動が含まれているためであろう．
4 ．手動読取り記録と自動読取り結果の比較・検証
これまで手動読取り以外の方法が無かったため，求められた Pi の諸元（発現時刻，周期及び振幅）の確度については，検証する手段がなかった．第 2 節の自動読取りによる Pi 検出を実行した結果と比較

図627することにより，初めて手動読取り記録の検証が可能となった．本稿では2012年の手動読取り記録のPi のうちクラス A 及び B の66事例について検証した．手動読取り記録のクラス C については，発現時刻以外の周期及び最大変化量の記載が無いことか
ら本検証では除く．最大変化量を示す各成分の割合
は，手動読取り記録では，X: 30，Y: 36であり，自動読取り結果では，X: 31，Y: 35であった．両者の食い違いは 5 例あった（手動読取り→自動読取り（内ファジィ・ロジックの結果．（a）図 5 b の時間周波数解析結果のローカル・ピーク，（b）Pi に関するメンバーシップ関数，及び（c）図 6 a のローカル・ピークの対応する図 6 b のメンバーシップ関
数値．図 4 の時刻に対応する．

28

笹岡雅宏
表1

自動読取り手法による Pi 検出リスト．各行の数値について，日付（yyyy：年，m：月，d：日），Pi 検出時刻
（h＊：時，mm＊：分，ss＊：秒），T＊：ウェーブレット関数の周期，Â＊: T＊に対応する Pi の最大変化量，Pi発現時刻（h：時，mm：分，ss：秒），T：ピーク間由来の周期，Â：T に対応する Pi 最大変化量，Q：通常時の Quality のクラス，Q＊：擾乱時の Quality のクラス．女満別観測点の地磁気毎秒値 Y（東西）成分（2012年 4 月 5 日08:00UT ～09:00UT）の事例を示す．黄色の行の検出結果は手動読取り記録に対応する．

表2

図7

手 動 読 取 り 記 録 に お け る Pi 検 出 リ ス ト． 女 満 別 観 測
点 の 地 磁 気 毎 秒 値 Y 成 分（2012年 4 月 5 日08:00UT ～09:00UT）の事例を示す．

検出された Pi の例示（赤点線）．女満別観測点の地磁気毎秒値 Y（東西）成分：
2012年 4 月 5 日07:50UT ～09:10UT を示す．矢印は手動読取り記録の Pi を示す．

図 8 図 7 中の Pi の検出に用いられたウェーブレット関数の例．

ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法

訳）：X → Y（ 2 例），Y → X（ 3 例））．
図 9 は，Pi 発現時刻の比較結果である．直線近似の関係から両者は非常によく対応しており，手動読取り記録の方がより早い発現時刻を見積もる傾向があることが示されている．95％予測区間を外れる
3 事例について図10に示す． 3 事例とも手動読取り
記録の方が発現時刻はかなり早い．これらの食い違
いは求め方の相違が関係している．もともと観測者の手動読取りでは連続発生する一連の Pi の開始と終了を一区切り（シリーズ）として，そのシリーズの開始を発現時刻として求めており，その区間における最大変化量を記録する．自動読取りでは，個々のインパルスの立ち上がりから発現時刻を求めており， 3 事例とも手動読取りの発現時刻前後でクラス
C の Pi を検出し，むしろ手動読取り記録はその Piの発現時刻に近い．即ち，手動読取り記録では Piとその発現時刻が必ずしも対応せず，クラス C のPi 検出数も少ない．また，手動読取り記録におい
表3

29て最大変化量を検出した時刻が記録されていないた
め，その Pi を毎秒値プロット上で確認するのがしばしば困難になった．
図11a は，第2.3節における地磁気変動のピーク間から求めた周期の比較結果であり，図11b は，ウェー
ブレット関数の周期をそのまま用いた比較結果である．手動読取り記録と同様に，X，Y 及び Z 成分のうち最大変化量を示す成分から求めた周期を選択している（最大変化量を示す成分が異なる 5 例については手動読取りと同じ成分）．直線近似は両者ともよい対応を示すが，手動読取りでは周期をピーク間から求めるため，図11b の方のバイアスが大きい．

地磁気観測所女満別観測点における Pi のクラスの判定基
準．Â は振幅（最大変化量）を示す．擾乱時とは当該 Piに他の現象が混在し波形を崩している場合を言う．

図9

自動読取り 結 果 と 手 動 読 取 り 記 録 の 発 現 時刻 比較．2012年 の 女 満 別 観 測 点 Pi の66事 例 の 比 較 結果．各図中に直線近似及び95％予測区間について記す．発現時刻は正時を基準に秒換算し，例えば，手動読取り23:13:00UT と自動読取り23:12:05UTの比較では23:00:00UT を基準に780sec と725secと し て， ま た， 手 動 読 取 り16:59:00UT と 自 動 読
取り17:00:50UT の比較では16:00:00UT を基準に3540sec と3650sec としてプロットした．

図10

Pi の発現時刻の比較．各図中，点線楕円で示された各 Pi に対して，青縦線は手動読取り記録，赤縦
線は自動読取りについてそれぞれ記す．紫縦線は，手動読取り記録に近い自動読取り（クラス C の Piの発現時刻）について記す．

30

笹岡雅宏

図11a の95％予測区間から大きく外れた 4 事例について，図12に示す． 4 事例とも，手動読取りの周期
は Pi 波形がダンプする区間で検出されている．また，図12c 以外の 3 例については，手動読取りの発現時刻の方が遅れているが，何故だろうか？図12cは，発現時刻以前の時間の地磁気が殆ど乱れていな
いため手動読取りでも地磁気変動の始まりがよく視認できたのだろう．それ以外の 3 事例については，
手動読取りの Pi 発現時刻以前の地磁気が乱れてい

図11

図12

るため，もっと大きな地磁気変動を視認できた時刻
を発現時刻として記録したのだろう．このため周期の見積りは発現時刻より後の変動周期から読取ることとなり Pi のダンプ区間で探す結果となったのだ
ろう．図12c については，両者とも Pi のピーク付近の周期を記録しているが，Pi の立ち上がりからピークまでの区間とダンプ区間では周期に相違が見られるため，両者が食い違う結果となった．一方，
これら手動読取り記録の遅れについては，ハイパス

自動読取り結果と手動読取り記録の周期比較．（a）ピーク間から求められた周期，（b）ウェーブレッ
ト関数で求められた周期（パラメータ）．各図中に直線近似及び95％予測区間について記す．

Pi の周期の比較．各図中，青線は手動読取り記録，赤線は自動読取り結果について，周期及び発現時刻をそれぞれ示す．

ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法

フィルターの位相特性が原因となっている可能性があることを記しておく．
図13は，X，Y 及び Z 成分についての最大変化量の比較結果である．直線近似は良い対応関係が見られることを示す．図11における手動読取り記録と自動読取りの食い違いは図12を用いて検証した通りであり，図13のばらつきの主な原因は同じ地磁気変動のピークから最大変化量が算出されていないことであろう．また，各バイアスについては，自動読取り
結果は手動読取り記録と違って数秒以下の短周期変動の影響を受けていないことが反映されているであろう．
手動読取り記録における Pi のクラスについては，
最大変化量と Pi 基準（表 3 ）との関係から，通常時の基準を用いた事例が53，擾乱時の基準を用いた事
例が 1 ，両方の基準で同じクラスになる事例が12ある（基準が不明）．基準が不明である事例に通常時の基準を用いると，66例中20例の食い違いが見つかった（手動読取り→自動読取り（内訳）：A → B（ 2 例），B → C（18例））．一方，基準が不明である
事例に擾乱時の基準を用いるとクラス C が増えて，
66例中23例の食い違いが見つかった（手動読取り→
自動読取り（内訳）：A → B（ 2 例），B → C（21例））．
概ね自動読取りの方の最大変化量が小さいことに由来する相違であろう．
5 ．まとめと考察
2012年における女満別観測点の Pi の諸元（発現時刻，周期，及び最大変化量）について，開発した
自動読取り手法を用いて求めた．自動読取り手順では，Pi の発現時刻と周期（パラメータ）をウェーブ
レット関数から得て，ピーク間の周期と最大変化量を多重解像解析の分解データから直接読取った．本手法においては，最大変化量の読取りを手作業に似せるために，振幅をウェーブレット係数から推定し

図13

31

なかった．アウトプットとしては，表 1 のフォー
マットの Pi リストのほか，時間周波数解析の全ての計算結果もファイルで出力しており，結果を詳細にチェックすることが可能である．手動読取り記録と自動読取り結果の比較については，包括的に非常に良い対応関係が見つかることが示された．周期については，ピーク間による時間幅ではなくウェーブ
レット関数で求められたパラメータそのものを採用してもよいと思うが，これまでの手動読取りの数値とバイアスが生じるため避けたほうがよいかもしれない（図11）．一方で，両者の比較により手動読取
り記録の内容を初めて検証することができた．手動
読取りは地磁気擾乱の影響を著しく受けていることがあるため，Pi の手動読取り記録の確度は低い可能性がある．自動検出では，Pi の形態的特徴を検出することにより地磁気擾乱の影響を受けずに Piの諸元を求めることが可能であり，この点は自動検出の利点であろう．一方で，手動読取り記録においては，振幅には数秒以下の地磁気短周期変動によるバイアスがあること，また，連続発生する一連の Pi の開始時刻は Pi の発現時刻と必ずしも一致しないことが分かった．Pi とその発現時刻の対応が不明瞭なことがあるため，利用には注意が必要であろう． 1 成分を 1 時間計算するのに要する CPU 時
間は約15分であり， 3 成分で46分程度であった．図
7 の事例の Pi 数であれば，自動処理のほうが手作業よりも短時間で解析が可能で効率が良い．勿論，最終的には観測者が Pi リストを取りまとめるため，自動検出した Pi が全て Pi として採用されるとは限らない．
現象全体の開始と終了の時刻決定については，Piの選択に依存する．例えば，図 7 の事例では開始
08:28UT，終了09:10UT と記録されているが，開始については08:13UT に見られる地磁気変動を Pi とみなすかどうかで変わるであろう．終了については，

自動読取り結果と手動読取り記録の X（南北），Y（東西）及び Z（鉛直）成分の比較．2012年の女満別観測点
Pi の66事例の比較結果．各図中に直線近似及び95％予測区間について記す．

32

笹岡雅宏

もっと早い時刻09:06UT とすることができるであろう．これについては観測者判断の方が優先することなので，今回の自動処理の議論の対象とはしなかった．
参考文献
新井康平，独習ウェーブレット解析，近代科学社，東京，
2006．
Cornman, L B., R K. Goodrich, C S. Morse, and W L.
Ecklund, A fuzzy logic method for improved moment
estimation from Doppler spectra. J. Atmos. Oceantic
Technol., 15, 1287-1305, 1998.
Kautsky, J. and R. T urcajova, A matrix approach todiscrete wavelets, in Wavelets: Theor y, Algorithms,

and Applications Volume 5, edited by C. K. Chui, L.
Montefusco and L. Puccio, Academic Press., Boston,117-135, 1994.
Nose, M., T. Iyemori, M. Takeda, T. Kamei, D. K. Milling, D.
Orr, H. J. Singer, E. W. Worthington, and N. Sumitomo,
Automated detection of Pi2 pulsations using wavelet
analysis: 1. Method and an application for substorm
monitoring, Earth Planet. Space, 50, 773-783, 1998.
笹岡雅宏，山崎 明，小山崇夫，火山性磁場の長期トレンド抽出のためのウェーブレット解析，Conductivity
Anomaly 研究会論文集，2011．
Takano, S., T. Minamoto, H. Arimura, K. Niijima, T. Iyemori,and T. Araki, Automatic detection of geomagnetic
sudden commencement using lifting wavelet filters,
Lecture Notes in Computer Science, 1721, 242-251,
1999.

ウェーブレット解析及びファジィ・ロジックを用いた地磁気脈動 Pi の諸元の自動読取り手法

An automatic detection method that uses wavelet analysis
and fuzzy logic to characterize elements of Pi
geomagnetic irregular pulsationsby
Masahiro SASAOKA
Kakioka Magnetic Obser vator yReceived 27 October 2016; received in revised form 23 March 2017;
accepted 28 March 2017

Abstract
To create an alternative to manual detection, an automatic detection method is developedthat uses wavelet analysis and fuzzy logic to characterize elements of Pi geomagnetic irregularpulsations (onset time, period, and amplitude). The values obtained by our new method arecompared with the records from manual detection using Memambetsu geomagnetic data for2012. Although the results show very good agreement overall, the biases due to components ofperiods less than 10 s in the high-pass filtered data with which we detect Pi by eye might cause
amplitude differences between manual and automatic Pi detection. Also, in the manual Pi
detection records, the start time of the Pi series does not necessarily agree with the onset timeof the Pi.

33