地磁気観測所テクニカルレポート 第２巻第１号 １−17頁 平成16年3月
Technical Report of the Kakioka Magnetic Observatory Vol.2, No.1, pp.1‑17, March2004

確率差分法を用いた火山性全磁力変動の抽出手法
藤井郁子（調査課）
2004年2月6日受付，2004年2月25日改訂，2004年2月26日受理

要旨

地磁気観測値から火山性の長周期変動を抽出することを目的に，超高層や外核起源の広域的変動を確率差分法を使って除去する手法を開発した．三宅島の全磁力に対し，神津島の全磁力，柿岡の磁場3成分を参照することで，局所的な磁化や誘導を含む幅広い波長を持つ変動に対応した．
また，確率差分法の係数の最適化にはAICを取り入れ，客観的な判断を自動で行うようにした．
残差に残留する電離層潮汐，海洋潮汐による成分は，Sqと海洋潮汐の周期を持つ正弦波の重ね合
わせをロバスト最小自乗法により見積もって除去した．従来の神津島との単純差を求める手法より，火山性ではない変動を有意に減らすことができた．
新しい解析手法を観測値に容易に適用するため，必要な前処理，後処理を含むプログラム群を作成し一連の作業を自動化して，2002年より三宅島の4観測点の全磁力毎分値に対して適用を開始した．阿蘇山，草津白根山など，他の火山にも順次適用を広げている．

1．はじめに火山活動によって火山体内部で温度や圧力が変動すると，熱消磁現象や圧磁現象を通じて地表で磁場変動が観測されることがある（Johnston, 1997）．これらの変化は比較的ゆっくりとしたもので，表面現象ではわからない火山深部の状態を捉えるのに適していることから，地磁気観測から火山活動をモニターする試みが多数行われている（例えばSasai et al.,2002）．
しかし，地表で観測した地磁気には，火山活動起源ではない成分が多く含まれており，火山監視にとっては深刻なノイズになっている．日本周辺での平均的な地磁気の強さは約50000nTであるが，これは外核で流体鉄がダイナモ運動することによって生じている．主磁場と呼ばれるこの成分は全球的な波長（〜数万km）を持ち，数nT/year程度で変動する（永年変化）．数年以下の周期帯には，電離層や磁気圏な
どの超高層に流れる電流による速い変動が常時存在
し，日本周辺では最大振幅にして数百nTに及ぶ．超高層電流変化によって地中に誘導された電流による誘導磁場も存在する．また，火山岩は強磁性体を含んでいることが多いので比較的強い磁性を示し，火山では局所的な地殻磁化の影響が数千nTに達する
ところもある．最近，地表付近の磁化が気温によって変化したという報告もあった（Utada et al., 2000）．
この他，海岸近くの観測点では，海水が高電気伝導度であることから，海流が地磁気中を流れることでダイナモ発電する効果も知られている．
熱消磁現象や圧磁現象に起因する火山性変動は1〜10nT/year程度で他の成分に比べて振幅が小さい
ので，磁場で火山活動をモニターするためには，観測値から火山性磁場変動を抽出する作業が必要になる．これまで使用されてきた抽出手法は，火山活動の影響がない場所で同時に測定した参照データを利用し，火山観測点データとの差を取ることによって，局所的な火山性変動だけ取り出すものが多い．これは，最も大きなノイズ源である超高層電流と永年変化が広域的（波長=数百〜数万km）という性質に着目したものである．
単純差は最も広く利用されている方法で，参照点と火山観測点の距離rが火山性以外の磁場変動の波
長Lに対しr≪Lである場合に有効である．この場合，
火山観測点と参照点で火山性以外の変動は同一と考えられるので，2点の磁場の差が火山起源のものと推定される．適用例を見ると，r≦10km程度で良好な結果が得られる事が多い．
しかし，観測された磁場が地磁気の大きさである
全磁力成分の場合，単純差には観測点周辺の局所的

2

藤井郁子

な 地殻磁 化の 影響が 現れ ること が知 られて いる（Beahn, 1976）．地磁気ベクトル B が時間変化しな
い成分 B と時間変化する成分 ∆B で構成されている

日変化が小さい深夜2〜4時の観測値を用いて単純差を計算していたが，観測点が海岸に近いことから
海洋潮汐の地点差が残留していた．長期的傾向を見とすると， B >> ∆B であることから，全磁力Fは近るために夜間単純差を日ごとに単純平均して1日値にしていたが，それによって海洋潮汐成分がエイリ

似的にアスし約15日周期の見かけの変動を導入していた．
これらの問題点を解決するため，確率差分法を用いた広域変動除去と電離層・海洋潮汐の除去を組み

F = B + ∆B ≈ B + ∆B ⋅

B
B

(1)

と表され，Fの時間変化（式1の右辺第2項）は B の方向によっている．B に寄与するのは主に主磁場と地殻磁化であるが，比較的強い磁化のため B が局所的になり，結果として全磁力変化が局所的になる．
全磁力は手軽に精度良く計測できることから火山観測に盛んに用いられており，観測点によっては単純差を補正する必要がある．田中他（1978）は，式1のFの時間変化が ∆B と観測点固有の定数係数ベク
トルの内積であることから，∆B が2点間で同一なら補正すべき成分も同様の形式になることに注目して，DI補正法を提案した．
火山観測点と参照点の距離rが大きくなったとき，あるいは，観測点の地下に誘導された磁場が無視できないとき，∆B が2地点で同一という仮定が成立しなくなり，単純差やDI補正では誤差が大きくなってくる．しかし，r≦Lであれば同一起源の変動は誘導
成分も含めて相関があり，マグネトテルリック法や磁気探査法では磁場や電場の間に各周波数で場所に固有の応答関数の存在を想定する．これは時間領域
でインパルス応答が存在することと同義であるので，Davis et al．
（1981）は，観測値Fを参照点の磁場3成分と同じ火山で観測された全磁力F1，F2，‥にインパルス応答を作用させた確率差分方程式でモデル化した．火山では変化を見つけるため狭い範囲に複数の観測点を設けることが多いので，この手法は広域的な特徴を参照点，局所的な特徴を火山観測網でモデルすることを考えている．火山観測網の中で火山性の変化が観測点ごとに独立でない場合は適用が難しい．
気象庁地震火山部火山課と地磁気観測所では，
2001年に三宅島4点，神津島1点の全磁力連続観測点を設置し，三宅島2000年噴火後の火山活動を監視してきた（藤井・大和田，2004）．火山性変動を抽出するために，地磁気観測所では，三宅島から30km離れた神津島観測点との単純差を採用したが，いくつか問題があることがわかった．まず，昼間に電離層が太陽熱で温められることから発生する日変化（電離層潮汐）が2地点で同一でないため，全日のデータの単純差では十分除去できていなかった．そのため，合わせた手法を開発し，2002年初めから三宅島の全磁力データ処理に適用を開始した．
その後，阿蘇山，
草津白根山の全磁力データに対しても，順次新手法が適用されつつある（気象庁地磁気観測所，準備中）．
本稿では，地磁気観測所で採用している新手法について解説すると共に，採用後2年間にわかってきた
問題点とその改善策についても簡単にまとめる．
2．新手法の概説
2.1 確率差分法を用いた広域変動除去手法
除去すべき変動の中で最も大きな振幅を持つのは永年変化と超高層電流であり，どちらも広域的な特徴を持つ．しかし，超高層電流は観測点下の地下構造を反映する誘導磁場を伴うことから，除去精度をあげるなら相関する短波長成分をも視野に入れる必要が出てくる．また，全磁力の場合は，広域的変化に局所的な地殻磁化が上乗せされる．
三宅島の観測を例にとって，永年変化と超高層電流の影響を除去する方法を考える．2001年2月に三宅島には金層（MYK1），大路池北（MYK2），手島牧場（MYK3），新澪池跡（MYK4）の4点，参照点として神津島（MYK5）の計5点の全磁力観測点が設けられた（図1）．最も近い磁場3成分観測点は，約180km離れた柿岡にある．神津島・柿岡と三宅島の距離は
永年変化の波長より十分短いので，永年変化に関してはこれらの観測点では同一とする．超高層電流は，
日変化の例でもわかるとおり，同一ではない．
ここから，時刻tiにおける三宅島のある観測点の全磁力変動を∆FMYK(ti)，神津島の全磁力変動を∆FKOZ(ti)，柿岡の磁場 変動のベク トル3成分を ∆X(ti), ∆Y(ti),∆Z(ti)と表記する．データは時間間隔∆tおきに計Ν個あり，ti=t1+(i-1)∆tである．
神津島，柿岡と三宅島の磁場変動の関係を調べるため，2001年3月〜11月の毎時平均値を使って，∆FMYK(ti)と∆FKOZ(ti)，∆FMYK(ti)と∆X(ti), ∆Y(ti), ∆Z(ti)の2組についてそれぞれコヒーレンスを計算した．三宅島の観測点のうち，島の東部のMYK1と南部のMYK4の場合を図2に示す．どの観測点・磁場成分の組み合
わせの場合でも，周期10000時間以下でコヒーレン
スはほぼ0.86を越える高い相関を示す．ほとんどの

確率差分法を用いた火山性全磁力変動の抽出手法

図１

三宅島と神津島の全磁力観測点（MYK1，MYK2，MYK3，MYK4，MYK5），柿岡の3成分観測点の分布図．

周期で，神津島の全磁力のほうが柿岡の磁場3成分よりも三宅島の全磁力とよく相関しており，特に10時間以下の短周期で顕著である．三宅島の変動の大部
分は柿岡と共通するような長波長のものであるが，柿岡にはないが神津島と共通の短波長成分もあることが示唆される．従って，神津島の全磁力，柿岡の3成分を同時に参照すると，短波長から長波長の成分までをカバーできることになる．
図2より，∆FMYK(ti)と∆FKOZ(ti)，∆X(ti)，∆Y(ti)，∆Z(ti)の間には周波数ごとに安定した関係が成立している
ことがわかるので，その関係を時間領域で確率差分法を使って表現すると，
Ｋ

Ｋ∆FMYK(ti)＝ ∑ aj∆XKAK(ti＋j)＋ ∑ bj∆YKAK(ti＋j)Ｋ

j = -MＫ

j = -Mてよい．従って，一度係数を決めておけば次からはフィルタリングの演算だけを行えばよいので，効率よく広域的変動を除去できる．実際には，同じ超高層電流でも静穏時と擾乱時では電流系の空間構造等に違いがあり，フィルター係数も全く同一にはならないが，その差は計算誤差に比べて無視できる程度とみなす．
最適な係数は2つの段階を経て決められる．まず，MとKの値について探査範囲を定め，その範囲中のすべてのMとKの組み合わせについて独立に，

σ2=

1 N
2
∑ δ (t i ) が最小になるよう最小自乗法で係N i =1数を決定する．次に，AIC（Akaike, 1973）を

j = -M＋ ∑ cj∆ZKAK(ti＋j)＋ ∑ dj∆FKOZ(ti＋j)＋δ(ti)j = -M3(2)のように書ける．M，Kは過去・未来のデータを参照する個数，aj，bj，cj，djは係数である．神津島と柿岡は三宅島の火山活動の影響を受けていないと考えられるので，残差δ(ti)に火山性変動が含まれる．
式２は，神津島の全磁力を使って広域的な変動を概算し，柿岡の磁場3成分を使って全磁力では表現できない成分を補正することを狙っている．
式2は，∆FKOZ(ti)，∆X(ti), ∆Y(ti), ∆Z(ti)に長さM+K+1のフィルターを作用させて∆FMYK(ti)を予測することを表現していると考えることもできる．フィルター
係数aj，bj，cj，djは，変動のタイプや変動ソースに対する観測点の位置関係，それぞれの観測点下の電気伝導度構造や磁化分布によるので，同種の変動で地下構造が変わらないなら時刻によらない値と考え

AIC = N ′ ln(2π σ 2) + 2 E ( M + K + 1) + N ′(3)のように計算し，AIC最小となるMとKの値の組み合わせを最適とする．ただし，Eは参照点での磁場成分の数で，式2のモデルではE=4である．N’はデータの個数で，データ個数による変動をなくすために
N’=N- (Mの最大値)- (Kの最大値)に固定した．本来の定義上はN’=N-M-K個であるが，実際に計算してみるとMやKが大きいほうがAICが小さくなる傾向があり，異なるモデル間の比較がうまくいかなくな
ることが多いので使わなかった．
AICは，右辺第1項が最小自乗法でのフィッティン
グ誤差，第2項がフィルター係数の総数の2倍，第3項がデータ個数（定数）という構成になっている．
一般に，フィルター係数の数が増えると，モデルが複雑になって細かな変動にも対応できるため，フィ

4

藤井郁子

ッティング誤差も小さくなる．しかし，不必要にモデルを複雑にすると本来説明する必要のない変動（例えば，観測誤差）にも対応してしまうので，モデルの複雑さとフィッティングの良さとの間に最適なバランスを見つける必要がある．AICはそのための統計的指標で，AICを最小にするということは，
データを十分説明するできるだけ単純なモデルを見つけることを意味する．
2.2

潮汐除去手法図2を見ると，特にMYK4で，周期約24時間と約12時間で三宅島と参照点の磁場のコヒーレンスの減少が見られる．図2と同じ期間のMYK4の全磁力のパワースペクトルを見ると，およそ24時間の整数分の1の周期に鋭いピークがある（図3）．これらの周期帯には，電離層潮汐と海洋潮汐が存在し，他の周期帯と異なる振る舞いを示すことが知られている．コヒーレンスが低いと前節の確率差分法はあまり有効でないため，残差δ(ti)に電離層潮汐や海洋潮汐の影響が残る可能性が高い．
前節の残差δ(ti)に電離層潮汐や海洋潮汐がある場
合，それらを除去する手法を本節で考える．潮汐が存在しない場合，あるいは，観測点と参照点で潮汐がよく相関している場合は，この節の手法を適用する必要はない．
電離層潮汐や海洋潮汐は，なぜ参照点を使ってうまくモデルできないのだろうか？電離層潮汐で顕著なSqは短波長で，誘導磁場の空間構造も複雑になる
ので，観測点間で差が出やすい．また，海洋潮汐では月や太陽の引力で海水が全球的に移動することか
ら広域的な海洋ダイナモ現象が引き起こされるが，海岸からの距離や海底地形の影響も加わるので，観図2測される磁場は局所的な要素を持っている．
電離層潮汐と海洋潮汐による磁場変動を簡単にモデル化することを考える．本節で必要なのは残差
δ(ti)に含まれる潮汐成分であり，観測された磁場に含まれる潮汐成分とは必ずしも同じではないが，基本的な特徴は変わらないとして考察する．電離層潮
汐と海洋潮汐で最も特徴的なのは，図3にも示された
とおり，周辺周期に比べて変動振幅が大きく，周期性が非常に高い波であるということである．
Sqは24時間周期の正弦波とその高調波の重ね合
わせでモデル化されることが多い（例えば，Malinand Winch, 1996）．この場合，磁場変化をSとすると，
P Q
  nmt + θ nm  ， (4)S (t ) = ∑ ∑ Anm sin 2π  +  T 1 T 2 
n =1m = −Qと単純化できる．ただし，T1=24時間，T2=8760時
間=365日で，n，mは日変化と年変化の高調波番号，
P，Qは高調波の最大次数である．Anm，θnmが各分潮の振幅，t＝0のときの位相を表している．年変化の影響を考えるのは，Sqの振幅等に顕著な年変化があることが知られているからである．
電離層にも引力による潮汐変化があることが知られているが，熱によるものに比べてはるかに小さいので，ここでは考えないことにする．
海洋潮汐は特定の周期と起潮力を持った分潮の集合として記述できる（萩原，1978）．M2（周期＝12.42059
時間），S2（12時間），K1（23.93452時間），O1（25.81924
時間）が4大分潮と呼ばれ，M2の起潮力が最も大きい．
その他多くの分潮が知られており，カタログ化されている．単純化のため周期だけに注目すると，磁場中の海洋潮汐成分S’は各分潮の周期を持つ正弦波の重ねあわせとして表現することができる（式5）．

2001年3〜11月の毎時値から計算した，MYK1の全磁力と神津島の全磁力（実線），MYK1の全磁力と柿岡の磁場X・Y・Z成分（点線）の間のコヒーレンス（左）．MYK4について同様の計算をしたもの（右）．

確率差分法を用いた火山性全磁力変動の抽出手法

5タに対して考慮した．次に順位の高いグループは，

図3

2001年3〜11月のMYK4の全磁力毎時値のパワースペクトルL  1 S ' (t ) = ∑ Al sin 2π  t + θ l  ，  T l 

l =1

(5)

ただし，lは分潮の通し番号，Lは分潮の数，Tl，Al，
θlはl番目の分潮の周期，振幅，t＝0のときの位相である．起潮力が大きいほうから15分潮を候補とする．
式4と5を比べてみると，全く同じ周期（例えばS2）の正弦波成分が含まれていることがわかる．したがって，周期だけでは電離層起源か海洋起源かを区別することはできない．地磁気，地電位差の潮汐周期
帯 に は こ の よ う な 特 徴 が あ る た め ， BAYTAP-G（Ishiguro, 1983; Tamura, 1991）を使った海洋潮汐分
離はうまくいかないことが多い．
そこで，ここでは式4と5をまとめてδ(ti)中の潮汐成分TをL  1 T (t ) = ∑ Al sin 2π  t + θ l   T l 

l =1

(6)

と表現し，式4と5に含まれる分潮周期から重複を除いた統一カタログを作成して，最小自乗によってそれぞれ振幅・位相を決定することにした．分潮間の周期の差が微小な場合は短い区間のデータでは区別がつけられないこと，起潮力の差によって分潮間に振幅の差があることが予想されること，などを考慮してカタログ中の周期に優先順位をつけ，順位の高い分潮からδ(ti)中に有意に含まれているかどうかを調べる方式にした．
調べる分潮に共通するルールは，（１）周期が解析期間の3分の1以下，（２）優先順位の高い分潮の中に周期の近いものがあるときそれに対して解析期間
中に0.22周期以上の差がつくこと，の2つである．Sq
起源のm≠0のものは，解析期間が1460時間以上の場合に考慮した．最も優先順位の高いグループは，Sq
起源でm=0かつ1≦n≦8のもの，海洋潮汐のうちS2を除く4大分潮で，ルールを満たす限りすべてのデー

海洋潮汐のうちQ1（26.86817時間），P1（24.06587
時間），N2（12.65832時間），K2（11.96726時間）， M3（8.2804時間），Sq起源のうち1≦|m|≦3のものである．最後のグループは，海洋潮汐のM1（24.833248
時間）， J1（23.098477時間）， OO1（22.306074時間），
2N2（12.871758時間），L2（12.191620時間）とSqの4≦|m|≦8で構成される．第2，第3のグループに分類される分潮は，ルールを満たせば考慮する．
具体的には，ある分潮周期の正弦波について，データにロバスト最小二乗法を適用し正弦波の振幅と位相を求め，それらが有意かどうか，見積もりエラ
ーの大きさ・残差の振幅の減少度・他の成分に与える影響から判定する．有意とみなされれば，データ中からその成分を取り除く．これを，選ばれた分潮について順に行う．ロバスト手法は，フィッティン
グ残差のうちガウス分布に従わないものに対して重みを作用させる方式を採用し，計算値が収束するまで繰り返した．重みには，始めの3回はHuber weight
を，それ以降は自然対数weightを採用した（Chave etal., 1987）．
3．データ処理プログラム
前章で解説した新手法を中核として，三宅島の全磁力観測値を自動的に処理するために使用されているプログラムについて解説する．このプログラムは，
データ形式，参照点名などに適当な修正を施せば，他の火山の観測データや他観測点の磁場・電場デー
タなどの処理にも利用することができる．
プログラムの目的は，以下の処理をできるだけ自動的に行うことにある．
（１）観測された三宅島・神津島の全磁力毎分値と柿岡の3成分毎分値から，異常値を除いて毎時平均値を計算する．
（２）三宅島の各観測点の全磁力毎時値について，
神津島・柿岡を参照点とした確率差分法のフィルター係数を決定する．
（３）三宅島の各観測点の全磁力毎時値に（２）の確率差分法フィルターを作用させ，参照点と相関のある広域変化を除去する．
（４）（３）の残差から，電離層潮汐，海洋潮汐起源の変動を除去する．
（５）（４）の残差にローパスフィルターを作用させ，毎日値を計算する．
（６）三宅島の各観測点の全磁力毎日値（残差）を図示する．
これらのデータ処理の流れをまとめると，図4のよう
になる．

6

藤井郁子

処理を3つのシェルスクリプトに分けたのは，すべての処理を毎回行いたいとは限らないので，自由に組み合わせられるようにするためである．例えば，環境が変わらないならフィルター係数は一度計算すればよいのでFprediction.shは一度だけ使用すればよく，2回目以降はmintohr.shとfiltering.shを連続して使用すればよい．
図4のうち，前章で取り上げていない処理a，b，fについて簡単に解説する．
3.1 異常値除去
生データから異常値を検出し，欠測扱いにすることを目的とする．毎時平均値を計算する前に明らかな異常値を除くことで，毎時平均値の質を高めることができる．プログラムmintohr.fに含まれ，三宅島と神津島の全磁力観測値に適用し，標準観測所である柿岡の観測値には適用しない．
2001年3月から11月の三宅島と神津島の全磁力データを調べたところ，異常値は主にプロトン磁力計のミスカウントによるものであった．ミスカウントによる異常値の特色は，（１）継続時間の短い単発の異常，（２）正常な値に対する異常な値のずれ幅がおおよそ決まっている，の2点が挙げられる．三宅島と神津島のシステムの場合は，ずれ幅は約45nTか，それの整数倍であった．ミスカウントを1つ見逃すと，毎時平均値に約0.75nTの差が生じることがわかる．
これらを踏まえ，
以下の3段階からなる検出法を採
用した．

図4

データ処理の流れ

図4の処理を行うため，6個のFortran77プログラム，
それらをまとめる3個のcshスクリプトを作成し，入力として最初のデータの年月，期間，観測点名を与えれば，すべての処理がほぼ自動的に進むようにした．図4の処理a，bには，三宅島と神津島用・柿岡用に そ れ ぞ れ 一 つ の プ ロ グ ラ ム （ mintohr.f,
mintohr.jma.f）が割り当てられ，シェルスクリプト
mintohr.sh で 動 作 す る ． 処 理 c は プ ロ グ ラ ム
Fprediction.fとシェルスクリプトFprediction.shを使う．処理d〜fとデータ図示はシェルスクリプト
filtering.shで一括して行われ，処理d〜fに対してそれぞれプログラムFfilter.f，detidesq.f，lowpass.fを開発した．図示はフリーソフトgnuplotをシェル中で起動して実行する．

（１）生データfi（i=1,…,N）から，1分前のデータとの差の数列（ffi = fi - fi-1）を作る．ffiとしきい
値 mcscale （ 現 設 定 で は 40nT ） を 比 較 し ， ffi>mcscaleなら1，ffi <-mcscaleなら-1，それ以外
は0と記録する．
（２）次に，1分後のデータとの差の数列（fbi = fi fi+1）を作り，同様に比較結果を記録．
（３）（１）と（２）の積が0でない番号iのデータを異常値とみなし，欠測値999999.0 nTに入れ替え
る．
2種類の前後差を使うのは，ミスカウントによる異常値とsscのような自然の信号を区別するのに，単発のスパイク状のものだけ選び出すためである．sscの場合は，スパイクというよりも急に上がってだらだら下がるような変化であることを考慮している．2
本のスパイクにはさまれた場合は，真ん中の値も異常としてしまうが，こういう場合の正常の判定はそもそも難しいので，そのままにしてある．
この方法によって対処できないタイプの異常は，ステップ状の変化，振幅の小さい変化，などである．

7

確率差分法を用いた火山性全磁力変動の抽出手法

3.2 毎時値計算
ミスカウントによる異常値を除去した後，毎分値
から毎時平均値を計算する．各時の0〜59分までの60個のデータのうち，欠測でないものを足し合わせて，足し合わせた個数で割る．欠測が30個を越えた
ら毎時値を欠測とする．
毎分値ではなく毎時値でこの後の処理を行うのは，目的とする火山性変動が長周期だからである．また，World Data Centerなどで流通する地磁気の毎時値は0〜59分までの毎分値を平均したものであることを考慮して，毎時値としては標準的な平均値を採用することにした．

図5

ローパスフィルターの係数

ローパスフィルターによる毎日値計算

フィルター長が長く，その分，データの始めと終わ

広域的変化，潮汐変化を取り除いた毎時値の残差v
から，
エイリアジングのない1日値を計算することを目的とする．火山性変動を議論するには必ずしも1
日値を計算する必要はないが，ゆっくりとした変動を強調して図示したいときのために，データ量を減らして扱いやすくしようとしてこの処理を取り入れている．

りが失われることである．この場合は，4日分ずつ失われることになる．

3.3

図3でも見たように，地磁気のデータには，電離層
潮汐や海洋潮汐など周期1日付近に鋭いパワースペ
クトルのピークを持つ変動があるので，1日値の計算には細心の注意が必要である．例えば，24時間分の生データを単純平均すると，24/n時間周期の変動は打ち消されるが，24時間周期の年周変動の高調波や
O1（25.81924時間）などの海洋潮汐成分は完全に0にはならない．他の変動に比べてこれらの変動はもともとの振幅が大きいため，単純平均には残った変動がエイリアジングして現れる．
ここでは，確率差分法による広域変動の除去，潮汐除去を行った全磁力毎時値vに対して，0時を中心にカットオフ周期が48時間のローパスフィルターを作用し，周期2日以下の変動を除去して，1日ごとにリサンプリングすることにした（式7）．
L

V (t i ) = ∑ α j v(t i + j ) ,i =1, 25, 49, …

(7)

j =− L

ただし，V (t i ) はローパスフィルターによって得られた毎日値，iはデータ開始からの通算時，αはフィル
ター係数，2L+1がフィルター長である．
ローパスフィルターはnon-recursive型で周波数領域で滑らかな特性を持つタイプ（Hamming, 1977）を採用した．図5に使用したフィルター係数を示す．
カットオフ周期が48時間の場合のフィルター長は147時間である．長所は，フィルターによる位相回
転がなく，データ中の短周期パワーピークによる影
響が小さいところ，短所は，non-recursive型のため24時間平均はフィルター長とフィルター係数の値がそれぞれ24と1/24に固定されているローパスフ
ィルターの一種であり，夜間値平均は同様に5と1/5の場合であるが，周波数特性は曖昧である．平均値
フィルターはフィルター長が短く演算が単純という利点があるが，結果に高い精度が必要な場合，地磁
気や地電流の毎日値には利用しないほうが無難と思われる．
4．三宅島の全磁力データの解析例
この章では，3章で概説したプログラムを使って三宅島の4観測点の全磁力観測値をデータ処理し，火山性磁場変動を推定した結果を示す．2001年3月から11月までの毎分値を用いて計算を行った．
4.1

異常値除去

異常値のしきい値mcscale=40nTとし，MYK1から
MYK5の全磁力毎分値について異常値判定を行った．
どの観測点においても異常はランダムにおこっており，特に異常な振る舞いはなかった．
図6にMYK2における2001年12月29-30日の全磁力
毎分値とその前後差を示す．40320-43200分に，急に増加して緩やかに減少する変化1回と，スパイク状の変化が5回見られ，前後差には計6個のスパイクがある．このうち後半4本が異常と判定された．
最初の2本のスパイク状変化が異常でないのは，ス
パイクではないのでffiとfbiが同時にしきい値を越えていない，ffiとfbiがしきい値より小さい，という理由による．しきい値を小さくすれば異常値と判定されるデータ数は増加するが，正常な値を異常と判断するリス
クも上昇する．見逃す影響と，正常を異常とする影響を比較しながら，しきい値を定める必要がある．

8

藤井郁子

図6

4.2

2001年12月29‑30日のMYK2における全磁力（上）と前後差ffi（点線）とfbi（実線）（下）．

確率差分法のフィルター係数の決定

三宅島の4観測点の全磁力データについて，
それぞれAIC最小となる確率差分法のフィルター係数を計算した．用いたデータは，大きな磁気嵐がなかった
2001年5月〜9月の4か月分である．期間が比較的短いことから，得られた係数は年周などの長周期変化には対応が不完全だが，数日以下の短周期変化には十分である．MとKは，それぞれ4〜30時間の間で変化させた．MとKの最大値の設定は，短周期変動に日変化が顕著なことから24時間が目安となると予想をたて，周期2日以下を重点的に除去できればよいと考えて，余裕を見て30時間に決めた．
4つの観測点とも，M=K=30のとき，AIC最小を示した．図7にMYK1に対する柿岡3成分，神津島全磁
力のフィルター係数を示す．最も大きな係数は神津島の全磁力の同時刻（j=0）の時であることがわかる．つまり，神津島の全磁力との単純差は，最も単純で効果的な近似である．ただし，図7から，j=0の係数は1ではなく，かつ，j≠0でも有意に0でない係数が求まっていることから，観測点間の距離や地下構造の違いを反映するとみられる時間差があることが示唆される．
島の南西の海岸にあるMYK4やMYK3では，異なる成分の係数間に強い依存が見られ，j≠0で比較的大きな振幅を持っていることがわかった．図8にMYK4のフィルター係数を示す．神津島のFと柿岡のXとZが，お互い補償するように動いている．これは，3つの参照データが完全に独立でないことや，確立差分法の係数決定法が係数の変化の滑らかさに関する項目を含まないことから表れる見かけの効果であると思われる．このような係数を用いて参照データと相関する変動を計算すると見かけ上大きな振幅になってしまうが，異なる成分間で相殺するので，最終的な残差には影響が表れない．
もし，参照データに他の成分に従属な成分が含まれているとしたら，多くの参照データを使ってモデ

確率差分法を用いた火山性全磁力変動の抽出手法

9の採用を見送った．
4.3 確率差分法による広域変動の除去
前節で計算したフィルター係数を用いて，三宅島の4点についてそれぞれ柿岡の3成分，神津島の全磁力と相関する変動を計算し，観測値との残差を求めた．
MYK1の全磁力毎時平均値の変動分と，図7のフィルターを用いて予測した各参照データと相関する成分，予測値を除いた残差を図9に示す．MYK1の全磁力変動は神津島の全磁力変動とよく似ており，柿岡
図7

MYK1用の柿岡３成分と神津島全磁力の確率差分法フィルタ
ー係数の3成分が補完的な役割を果たしていることがわかる．残差振幅は数nT以下であり，確率差分法によって約400nTの広域的変動が除去された．
比較のため，神津島の全磁力との単純差と図9の確率差分法残差とを図10に示す．単純差では，主に昼間に約10nTの日変動が残留し，磁気嵐（720−744時間）では20nT以上の変動が見られる．これは神津島と三宅島との距離が約30kmとやや遠いため，広域変動の振幅差や位相差が無視できないためと考えられる．確率差分法では，過去から未来までの参照値を使って，これらの振幅差・位相差の影響を数nT以下に抑えている．
さらに，単純差と確率差分法の違いをパワースペ

図8

MYK４用の柿岡３成分と神津島全磁力の確率差分法フィル
ター係数

ルを複雑にしたのに無駄な情報だったためフィッテ
ィングに改善が見られない訳であるから，式3のAICで判定できる．例えば，柿岡の3成分と神津島の全磁力を参照したとき（E=4）と，柿岡のXY成分と神津島の全磁力を参照したとき（E=3）で，最小となる
AICを比較すれば，1成分多い効果が統計的に有意かどうか判定できる．MYK４について，参照データの組み合わせを変えてAICが最小となる点を求めたと
ころ，図8の4成分の場合が最小であった．4成分は独立でない部分もあるが，完全には従属していないということを示唆するものと思われる．
見かけの変化を技術的に取り除くには，Egbert(1992)の方法がある．確率差分法は時間領域でフィルター係数を計算するが，Egbert(1992)ではまず周波数領域で応答関数を計算してから時間領域に変換する．その際に，地磁気や地電位の応答関数が，潮汐周期帯を除いて，周波数に対して滑らかに変化す
る性質を利用して，周波数領域での滑らかさを条件としたフーリエ変換を適用することによって係数上の見かけの効果を抑える．一部のパラメーターの選
択に手動の部分が残っていたため，今回はこの手法クトルで比較すると（図11），どちらの手法も観測値の変動をパワーにして1〜2桁小さくする良い近似であるが，周期100時間以下では確率差分法のほうがさらに1桁よい近似であることがわかる．
確率差分法のパワースペクトルは，潮汐周期帯を除いて，周期に対してほぼ単調に増加しており，長周期変動を扱いやすい状態になっている．
4.4

Sq・海洋潮汐の残留成分除去

図11でもわかるように，三宅島の場合，確率差分法
で広域的変化を除去した後の残差には電離層潮汐・海洋潮汐の残留成分が含まれている．残差を調べ，電離層潮汐・海洋潮汐の周期を持つ変動を検出したら，ロ
バスト最小自乗法を使ってそれらを除去した．
潮汐の残留成分には島内の観測点で顕著な差があることがわかった．島の南西部にあるMYK4やMYK3の残差には，目で見て潮汐変動が判別できるほど大きく残っており，潮汐除去の効果が著しかった．図12に潮汐残留成分の振幅が最も大きかったMYK4の結果を示す．MYK4では，K2の振幅が最も大きく（1.02nT），M2，S1，K1と続き，計28分潮の合計最大振幅は約5nTに達した．残留潮汐成分除去の効果をパワースペクトルで見ると，図13のようになる．
周期24時間，12時間，8時間などの周辺にあった鋭いパワーピークのみが除去されていることがわかる．

10

藤井郁子

図9

上から下へ，2001年3月〜11月のMYK1の全磁力，神津島の全磁力と相関する成分，柿岡のX,Y,Z成分に相関する成分，残差．

確率差分法を用いた火山性全磁力変動の抽出手法

11

一方，MYK1には潮汐成分はあまり残っておらず（図14），計8分潮が検出され，S1の0.13nTが最大であった．このような残差の性質の違いも，図7〜8に示した確率差分法のフィルター係数の振る舞いの差の原因の一つかもしれない．

図10

MYK1の全磁力毎時値について，神津島の全磁力と柿岡の3成分を使って確率差分法で処理した残差（実線）と神津島の全磁力との単純差（破線）を比較したもの．横軸は2001
年3月1日0時（JST）からの時間．

4.5 毎日値
前節までの処理で，全磁力毎時値から広域的変化と潮汐残留成分が取り除かれ，火山性の長周期変動を含む残差を得ることが出来た．このままでも火山性変動について議論することはできるが，長周期変動をより強調したい時のため1日値を計算する．ここでは，全磁力毎時値の残差に，0時を中心としてロー
パスフィルターを作用させ，1日値にリサンプリング
した．
ローパスフィルターの作用をMYK4のパワースペ
クトルで示したのが図15である．毎時値残差とロー
パスフィルターによって得られた毎日値のパワース
ペクトルは，周期48時間以上でほぼ一致しており，
不適切な影響を与えることなく毎日値の演算ができたことを示している．

図11

2001年3月〜11月のMYK1の全磁力毎時値（一点鎖線），神津
島の全磁力との単純差（破線），神津島の全磁力と柿岡の3成分を使って確率差分法で処理した残差（実線）のパワー
スペクトル．

図12

2001年3月〜11月のMYK4の全磁力について，確率差分法での残差（上），残差中のSｑ・潮汐成分（中），残差からSq・潮汐成分を除いたもの（下）
．

12

藤井郁子

4.6 まとめ本研究で開発した手法全体の効果を見るため，従来の手法での結果と，本研究で開発した手法の結果を比較した（図16）．従来は，神津島の全磁力との単純差を夜間平均して，三宅島の火山性磁場変動を見積もっていた．従って，2つの手法の差は，確率差分法，潮汐残留成分除去，ローパスフィルタリング，の3つである．三宅島の4観測点のすべてにおいて，新手法のほうが従来の方法より短周期変動が有意に減少し，全体の振幅が小さくなっている．特に，磁気嵐の時の残差（30~50日）と，MYK4に顕著な15日程度の振動が，新手法では取り除かれており，従来の手法では不要成分が残留していたことを示している．

図14

図13

2001年3月〜11月のMYK4の全磁力について，確率差分法での残差（点線）と残差からSq・潮汐成分を除いたもの（実線）のパワースペクトル．

2001年3月〜11月のMYK1の全磁力について，確率差分法での残差（上），残差中のSｑ・潮汐成分（中），残差からSq・潮汐成分を除いたもの（下）．

確率差分法を用いた火山性全磁力変動の抽出手法

13は，標準観測所の近くにある観測点のデータ異常をモニターする，地磁気観測に対する傾斜や温度の影響を解析する，地磁気を参照して地電流観測のノイズを検出する，などの用途に転用できるものであり，大きな可能性を持つと考えられる．
しかし，手法採用から2年を経過し，当初想定とは異なる振る舞いなど，いくつか改良が必要な点も明らかになってきた．

図15

2001年3月〜11月のMYK4の全磁力について，確率差分法と潮汐除去を行った毎時値の残差（点線）と，残差からロー
パスフィルターによってリサンプリングされた毎日値（実線）のパワースペクトル．

5．議論
前章での三宅島の全磁力への適用結果から，新しく開発した手法は火山性変動を抽出するため観測値から不要な広域的変動を除去する役割を果たしていることが示唆された．加えて，草津白根山，阿蘇山の全磁力観測値に対しても新たに適用され，良好な結果を得ている．適切な修正を施せば，今回の手法

図16

問題点の1つは，欠測値に対する対処が手薄なため，1つでも欠測値があるとその影響が最終結果に大きく反映されてしまうことである．これは，一連のデータ処理の中に比較的長いフィルター長をもつフィルタリングが2回（確率差分法，ローパスフィルター）含まれていることによる．フィルター長の中に一つでも欠測が含まれていると式2と7の演算は実行できないので，フィルター長だけ欠測の影響が続くことになる．
観測環境の良い三宅島や神津島でも工事やバッテ
リーのトラブル，ミスカウントなどでデータが使えなくなる状況が時々発生しており，野外の無人観測の場合はある程度の欠測が生じることは避けられないと思ったほうがよいので，欠測値に対する対処は2001年3月〜11月のMYK1，MYK2，MYK3，MYK4の全磁力について，従来の手法（上）と新手法（下）で計算した毎日値の比較．

14

藤井郁子

考慮すべきことであろう．特に，データ2，3個分以下の短い欠測は，長周期の火山性変動の見積もりに重大な影響を与えることは少ないので，補間することを検討したほうがよい．線形補間，スプライン補間，Auto regressive法による補間など，単純な手法
で十分役割を果たすと思われる．
二つ目の問題点として，データ中に異常値があった場合の影響があげられる．今のところ，毎分値から毎時平均値を計算する際に，ミスカウントを主な対象としたスパイク状異常値の処理だけを行っている．その後の計算では潮汐残留成分除去以外はロバスト手法を取り入れていないので，仮に異常値が残留していれば計算に影響する可能性がある．確率差分法の計算時にロバスト手法を導入することも考えたが，その場合のAICの定義が自明でなかったため，
今回は見送った．三宅島は観測環境が良く大きな影
響は出ていないが，阿蘇山や草津白根山では車の通行など人工ノイズが頻繁に観測されており，今後の改良が必要であろう．
三つ目には，手法全体の構成がやや複雑であることがあげられる．本論文で取り上げた手法は，観測された全磁力変動Fobsを
F obs (t i ) = F ref (t i ) + F tide (t i ) + v(t i )(8)のように，参照点と相関のある広域的変動Fref，潮汐残留成分Ftide，長周期の火山性変動vに分け，FrefとFtide
をそれぞれ確率差分法とロバスト最小自乗法で見積
もり，残差をvとしている．2段階の方法にしたのは，当初，海から遠い火山や参照点が近い火山ではFtideを考慮する必要がないと予想し，三宅島のように必要がある火山の場合のみの対応と考えていたことによる．その後，内陸の火山に適用したり他の観測データを調べるうちに，多くの場合にFrefとFtideの両方が必要であることがわかってきた．また，2段階の方式にしたことによって，データによっては確率差分法でのAICの収束があまり良くなくて長いフィルターを最適とする傾向があった．
以上のような問題点を踏まえ，将来的な手法のバージョンアップの検討をしているところである．現在検討中の手法（Fujii and Kanda, 2003）では，
F obs (t i ) = F ref (t i ) + F tide (t i ) + v(t i ) + ε (t i )(9)のように，観測誤差εを設けvの長周期特性を強調す
る一方，カルマンフィルターなどの解法を使ってFref，Ftide，v，εを同時に見積もることで高速化，高精度化を達成しようとしている．また欠測値の補間もできる方法である．

6．まとめ
三宅島の全磁力観測値から火山性変動を抽出するために，超高層や外核起源の広域的変動を確率差分法を使って除去する手法を開発し，前処理，後処理を含むデータ処理全体をプログラム化した．確率差分法では，三宅島の全磁力に対し神津島の全磁力，
柿岡の磁場3成分を参照して，局所的な磁化や誘導を含む幅広い波長を持つ変動に対応できるようにし，係数の最適化にはAICを取り入れて判断の客観性と自動化を達成した．さらに，残差に残留する電離層
潮汐，海洋潮汐による成分についても，Sqと海洋潮汐の周期を持つ正弦波の重ね合わせをロバスト最小
自乗法により見積もって除去した．これらの手法により，従来の神津島との単純差を求める手法に比べて，火山性ではない変動を有意に減らすことができた．新しい解析手法は，阿蘇山，草津白根山など，
他の火山にも順次適用が広がっており，適用例が蓄積してきた．欠測処理，異常値処理，全体の単純化などで，今後の改良が必要と思われる．
謝辞
本手法の開発にあたっては，瀧沢技官と伊藤技官に協力していただきました．長谷川技官，山田主任研究官には，運用時の問題点の指摘や原稿のチェックなどをしていただきました．また，編集委員の皆様
からは原稿について有効な助言を受けました．ここに感謝します．
参考文献
Akaike, H., Information theory and an extension of the
maximum likelihood principle, in: B. N. Petrov
and F. Csaki, eds., 2nd Inter. Symp. on Information
Theory, Akademiai Kaido, Budapest, 267-281,
1973.
Beahn,
T.
J.,
Geomagnetic
field
gradient
measurements and noise reduction techniques inColorado, J. Geophys. Res., 81, 6276-6280, 1976.
Chave, A. D., D. J. Thomson, and M. E. Ander, On therobust estimation of power spectra, coherences,and transfer functions, J. Geophys. Res., 92,
633-648, 1987.
Davis, P. M., D. D. Jackson, C. A. Searls, and R. L.
McPherron, Detection of tectonomagnetic eventsusing multichannel predictive filtering, J. Geophys.
Res., 86, 1731-1737, 1981.
Egbert, D.E., Noncausality of the discrete-timemagnetotelluric impulse response, Geophysics, 57,
1354-1358, 1992.

確率差分法を用いた火山性全磁力変動の抽出手法

Fujii, I., and W. Kanda, New procedures to extract
geomagnetic field variations caused by volcanic
activities, Earth Plant. Sci. Lett., submitted.
藤井郁子，大和田毅，全磁力観測，気象庁技術報告「平成12年（2000年）三宅島噴火及び神津島・
新島近海地震調査報告」，気象庁，印刷中.
萩原幸男，地球重力論，共立出版，東京，242pp，
1978.
Hamming, R.W., Digital Filters, Prentice-Hall,Englewood Cliffs, 1977.
Ishiguro, M., H. Akaike, M. Ooe, and S. Nakai, A
Bayesian approach to the analysis of Earth tides,
Proc. 9th Int. Sympos. Earth Tides, New York, pp.
283-292, ed. Kuo, J. T., Schweizerbart’sche
Verlangsbuchhandlung, Stuttgart, 1983.
Johnston, M. J. S., Review of electric and magnetic
fields accompanying seismic and volcanic activity,
Surv. Geophys., 18, 441-475, 1997.
気象庁地磁気観測所，阿蘇山における地磁気全磁力変化，火山噴火予知連絡会会報，第86号，準備中.
気象庁地磁気観測所，草津白根山における地磁気全磁力変化，火山噴火予知連絡会会報，第86号，
準備中.

15

Malin, S.R.C., and D. E. Winch, Annual variation of thegeomagnetic field, Geophys. J. Int., 124, 170-174,
1996.
Sasai, Y., M. Uyeshima, J. Zlotnicki, H. Utada, T.
Kagiyama, T. Hashimoto, Y. Takahashi, Magnetic
and electric field observations during the 2000
activity of Miyake-jima volcano, Central Japan,
Earth Plant. Sci. Lett., 203, 769-777, 2002.
Tamura, Y., T. Sato, M. Ooe, and M. Ishiguro, A
procedure for tidal analysis with a Baysian
information criterion, Geophys. J. Int., 104,507-516, 1991.
田中良和，増田秀晴，河村まこと，大地洸，加藤誼
司，馬場広成，吉野登志男，阿蘇火山地域での全磁力測量―(1977)，阿蘇火山の集中総合観測（第１回1977）報告，31-39，1978.
Utada, H., M. Neki, and T. Kagiyama, A study of annual
variations in the geomagnetic total intensity withspecial attention to detecting volcanomagnetic
signals, Earth Planets Space, 52, 91-103, 2000.