62

２．２ ノイズ除去手法
2．2．1 デジタルフィルター
2．2．1．1 はじめに地殻活動に関連する電場変動は必ずしも他の原因による変動を凌駕するとは限らないので，検出のためには他の成分を正しく評価し，時には除去することが必要となる．地殻活動に関連する電場変動を観測したという報告は数えきれないほどあるが，一方で，安易に異常な変化を地殻活動に結び付けているという批判が多いことを考えあわせても，他の成分の評価を慎重に行いどのような信号ならば観測可能かを調べることは不可欠と考える．
特に，これまでに示されたように，淡路島の電磁場データには電車によると思われる非常に振幅の大きい人工ノイズ（100mv/km，1nTのオーダー）が含まれる．外部磁場変動に伴う成分も，広い周波数帯に常に存在し，10mv/km，10-100nTのオーダーをもっている．一方，これまでに地震に伴う電場変動として報告されたものはおおよそ10mv/km，1nTのオーダーであり，観測した生データでは識別できないと予想されるので，人工ノイズと外部磁場変動に伴う成分をできるだけ正確に除去することが本研究の重要な課題となる．
徳本（1999）は，大谷を基準とした各測線の電位差について同時刻のデータを各月ごとに平均し，人工ノイズに同時刻同位相の成分があることを示した．しかし，平均で得られた成分は観測されたものに比べてはるかに振幅が小さく，これだけではノイズ除去はできないとされた．第2．1章でもわかるよ
うにノイズは複雑な周期特性を持っており，調和的な結果と考えられる．さらに徳本（1998）は，回帰分析による方法も試みている．大谷-奥畑測線の地電位差の毎秒値について，同時刻の他測線の地電位差や大谷の地磁気3成分を参照することで変動を説明しようとしたところ，地磁気変化より地電位差変化のほうが残差を小さくすること，地電位差変化のなかでも距離的に近い測線のほうが残差を小さくすることがわかった．この結果は，第2．1章からも示唆されるように，地電位差と地磁気の人工ノイズ間は位相差のない成分が卓越しているので，同時刻のデータを参照することで効果をあげることができたと解釈できる．ノイズ成分の残りや誘導成分には位相差があるので，同じ物理量で距離的に近い観測値のほうが変動が似ているのは自然である．別の時刻の観測値も取り入れたり，人工ノイズがない観測値（例えば，柿岡の磁場）を使ったりすれば，これらの性質がより明らかになるであろう．自己回帰分析
による方法の結果（徳本，1998）も示唆に富んでい

る．地電位差の毎秒値に対して，20-60個程度の回帰係数で振幅の約90％が説明できることがわかる．
図2.1.1や図2.1.10でもわかるように，ノイズの連続スペクトル部のパワーは周期100秒前後で頭打ちとなるので，自己回帰係数の個数が主に反映したのはこの部分だと考えられる．観測値に著しいラインス
ペクトル成分が存在する場合はその成分の周期までの自己回帰係数を取り入れると変動を効果的に説明できることから，淡路島の人工ノイズについては
500秒程度まで調査すればよいと予想される．徳本（1998）では最大何個の係数まで調査したのか不明
であるが，AIC（Akaike, 1973）を用いて最適個数を決めているので小さめの個数が選ばれた可能性もある．
本節では，地殻活動に関連する成分を取り出すた
め，できるだけ単純な方法で淡路島の電場変動から不要な成分を除くことを目的とし，定常時系列解析
によるディジタルフィルターの考え方を適用する．
非定常性が非常に強い場合などにはうまく対応でき
ないが，多くの場合におおまかな見積もりをするには有効な方法である．この単純な方法でどの程度まで人工のノイズやその他の不要な成分が見積もれる
か調べることで，より高度な方法をどのように適用すればよいかの指針が得られることが多い．
2．2．1．2 電磁場変動のモデル化
本節で用いる考え方を整理するため，電磁場のモデル化を行う．一般に，観測された電場Eobs，磁場
Bobsは，=Etc +Esqt +Eind +Ecn +εe
obs =Btc +Bsqt +Bind +Bcn +εbobs（2.2.1）（2.2.2）のように表せると考えられる．右辺の4項のうち，最初の3項が地球物理学的なシグナルで，残り2項がノイズを表す．シグナルとは，地殻活動に伴って生，Sqと潮汐に関連する成じる電磁場成分（Etc，Btc）
分（Esqt，Bsqt），外部磁場変動により地中に電流が誘導される現象に関連するもののうちSq・潮汐を除いた成分（B ind，E ind）で，本観測での目的はE tc（B tc）を検出することである．ノイズに2種類あるのは，通常の独立な観測ノイズεe，εbと，電場磁場に相関を持って現れるコヒーレントノイズE cn ，Bcnを区別したからである．観測ノイズεe，εbは互
いに無相関であり，理想的にはガウス分布する白色
雑音だが，現実には，電極の変動や原因不明の異常値なども含まれる．コヒーレントノイズEcn，Bcnを生み出す代表的な原因としては直流電車による漏え63

い電流があげられる．
E tcを検出するには，他の成分と正しく区別しな
くてはならない．区別するためには，それぞれが際立った違いを持っているか，あるいは独立な参照データが必要である．それぞれの成分について順番に性質を考えてみる．
E tcとB tcは生成メカニズムにより様々なパターンを示す．本観測では周期3秒以上のシグナルしか計測していないので，E tcとB tcとしては比較的ゆっく
りした変化を狙っていることになる．また，E tcとB tc 以外の成分が非常に大きい（100mv/km，10100nTのオーダー）ため，ある程度の振幅を持つ変化でないと検出できないと予想される．というのは，例えば，なんらかの解析により不要な成分の99％が取り除けたとしたら数値的にはノイズ除去成功と言えるが，それでも不要な成分の振幅が100mV/kmの場合は1mV/kmが取り残されることになり，十分な
S/Nを持ってE tcを検出しようとすればE tcが残った
ノイズレベルより有意に大きい必要があるからである．これらの要請から，本観測でE tcとB tcを検出できるとしたら，周期にして秒から年，振幅にして10mV/km以上か5-10nT以上のものと考えてよいであろう．Johnston（1997）のレビュー論文からは，（1）岩石中の水の移動に伴う流動電位現象，（2）圧磁効果，（3）電荷を生む物理過程（圧電効果など）
をこれらの要請を満たすメカニズムとして挙げることができる．（1）-（3）について簡単に検討する．
流動電位現象では数10mV/km，2-3nT程度のゆっくりした電磁場変動を期待でき，かつ，実際に水の動きは地殻活動に深く関わっているのではないかと考えられている．磁場変動の振幅は検出限界程度で観測できるか疑問だが，電場変動のほうは観測できる可能性がある．
（2）の地殻活動に伴う圧磁効果は，理論，観測共に数nT以下の磁場変動を報告している．本観測では電場変動を捕えるような設計がなされている点や，淡路島北部が磁化と磁気異常が小さ
く圧磁効果が出にくい地下構造をしている点（Electromagnetic Research Group for the 1995
Hyogo-ken Nanbu Earthquake, 1997）を考えると，
この効果を検出できるかどうかはわからない．（3）の電荷を生む物理過程では理想化された系で実験やモデル計算がされており，実際の地球にどのようにあてはめられるのか判断が難しい．多くのモデルでは，応力変動などによって断層面等に生じた分極が地中を伝播して電場変化につながる．地表に到達する前に減衰する，振幅が小さい，などの欠点がある
が，減衰を抑えるような地下構造や補償回路を考えればよいとの反論もある．

上の議論から，本観測として検出できる可能性があるのは（1）あるいは（3）で，S/N比を考慮すると，第一次近似として電場変動Etcのみを考え，
Btc ≈ 0（2.2.3）を想定できることが示唆される．本節では，簡単のため，式2.2.3が成立すると仮定する．
Sqと潮汐に関連する成分E sqtとB sqtは，およそ24
時間，12時間，8時間，6時間，…の特定の周期にのみ存在する成分で，すべての周期は既知である．
Sqや潮汐のうちいくつかの周期は全く同じであり，
電磁場データの中から2つのソースによる成分を分離するのは現実的には簡単でない．例えばBAYTAP-G（石黒他，1984）などは海洋潮汐を分離するプログラムだが，電磁場データの場合，Sqと海洋潮汐が同時に存在する24時間，12時間の潮汐の値がうまく求められないことがあることが知られている（第3．1章参照）．このような事情から，本研
究では，なんらかのモデルをたててSqと潮汐成分を求めるのではなく，電磁場データ中の特定の周期の成分を取り除くという方針でE sqtとB sqtの除去を行うことにする．
統計的な観点から平均的な外部磁場変動を考えると，Sq・脈動などの特定の周波数や磁気嵐の時を除いて，仮想的に磁気圏の定常環電流を想定すればソースとしてよい近似であることが知られている．
都市ノイズがない地域で電磁場変動を観測した場
合，周期2年以下で卓越するのはB indとE indである．
平均的なソースの場合，地表でのソース波長は数千kmであり，地域的な探査では空間的に一様な水平磁場と近似でき，BindとEindの水平成分ベクトルBhindとEhindは周波数領域で線形関係
Ehind (ω)=Z(ω)Bhind (ω)（2.2.4）

ただし，
Z(ω)=

Zxx

Zxy

Zy x

Zy y

で結ばれている．Z（ω）はインピーダンステンソ
ルと呼ばれ，外部磁場変動の波数と地下の電気伝導度構造から一意に決まる．ソースが定常的な磁気圏
環電流で，地下構造も時間的に変わらないなら，Z(ω)も時間的に不変である．そのような定常的なZ(ω)を10-5〜102Hzで求めて地下の電気伝導度構造を推定するのがMT法である．Bindの3成分の間には
経験的に，

64

とにする．即ち，周波数領域で定常的に
ind
ind
B ind
z (ω)= A(ω)B z (ω)+B(ω)B y (ω) （2.2.5）

Ehcn(ω)=Zcn(ω)Bhcn(ω)の関係が成立することが知られており，地磁気変換
関数A，Bから地下の電気伝導度の不均質を推定する手法をGDS法という（第3．3節参照）．
式2.2.4を時間領域に変換すると，例えばX成分に
ついて，
N

M

i= – Ni= – M

ind
ind
E ind
x (t)=∑ Zxx(i )B x (t + i∆t )+ ∑ Zxy(i )B y (t + i∆t )（2.2.6）と表される．zxx(i)，zxy(i)はインパルスレスポンス
で，∆tはサンプリング間隔，N，Mは打ち切り数である．zxx(i)，zxy(i)は，Bindに対するEindのディジタ
ルフィルターと見ることもできる．Z(ω)が定常的なら，z xx (i)，z xy (i)も時間変化しないので，一旦zxx(i)，zxy(i)を求められればどの時間tのデータにも適用でき，磁場データから電場の外部磁場変動による成分を見積もることができる．
都市部周辺で電磁場観測を行うと，主に短周期帯で異常なインピーダンステンソルが得られることがある．これは都市活動に伴うノイズが測定値に混入し，電場と磁場の水平ベクトルの線形関係がEhind (ω)+Ehcn(ω)=Zdb(ω)(Bhind (ω)+Bhcn(ω))（2.2.7）のように表される場合である．EhcnとBhcnが混入したため，見かけのインピーダンステンソルZdbはZとは異なってしまう．E hcnとB hcnが互いに無関係ならZ＝
Z dbとなるはずであるが，電車ノイズなど多くの人
工ノイズではEhcnとBhcnの間になんらかの相関関係が成立している．コヒーレントノイズという呼び名はそこから来ている．サンフランシスコなど直流電車が運行されている地域でこの現象が観測され，イン
ピーダンスが特におかしくなる0.1〜1Hzはデッドバ
ンドと呼ばれた（Fraser-Smith and Coates, 1978;
Mackie and Madden, 1992; Egbert, 1997）．電車が影響を及ぼす空間スケールは数10km程度と考えることが多いが，最近の研究で300km離れた地点で無視できないノイズをもたらしている例も報告されている（Egbert et al., 2000）．
第2．1節で示されたように，淡路島の電磁場の人
工ノイズの性質は空間的，時間的に比較的安定しており，著しい周波数特性を持っている．そこで，本節ではEhcnとBhcnの間に最も単純な関係を仮定するこ（2.2.8）が成り立つことを前提とする．このZcnは図2.1.4-6で示されたインピーダンスにほぼ相当し，それらの計算でエラーバーが小さくコヒーレンシーが非常に高かったことから，妥当な仮定であると考えられる．
ただし，昼と夜でノイズの性質が大きく違うことはわかっているので，昼と夜を区別してZ cnを求める
ことにする．
ところで，式2.2.4，8は人工ノイズと外部磁場変動に伴う電磁場成分の水平成分に関係があることを述べているのであって，Z成分が関係を持たないことを表現するものではない．磁場のZ成分はMaxwell方程式を通じて電場と関係すると同時に磁場のX，Y成分と関係する（例えば，式2.2.5）．そのため，単純に式2.2.4，8で磁場を3成分に拡張して電場水平2成分との線形関係を計算すると入力間相関が高くて計算が不安定になり，かつ，全体のコヒー
レンスは取り立てて良くならない．そこで，磁場のZ成分を利用することを考える場合は，他の成分と入れ替えるなどの方法をとることになる．
観測値に含まれるある成分を見積もりたい時，参照点での磁場B refとの関係からデータを見ると便利なことが多い．B refは，ターゲットとする成分と独立な観測ノイズしか含まない観測量で，B refと比べることによって観測値に含まれるターゲット成分を見積もることができる．
例えば，Bobs中の外部磁場誘導成分を見積もりたいとすると，外部誘導成分B'ind，独立な観測ノイズ
ε'bを用いて，Bref =B′ind +ε′b（2.2.9）と表される参照磁場B refを探す．外部磁場変動が地球規模あるいは探査地域より十分大規模な現象と考えられる周期帯では，B indとB'indの間に周波数領域
で定常的な線形関係，Bind (ω)=Z bB′ind (ω)（2.2.10）が存在する．Z bは2観測点の地下構造と外部磁場変動の波数の関数である．B'indは式2.2.10と2.2.4を通じてEindとも線形関係
Eind (ω)=Z′b(ω)B′ind (ω)（2.2.11）

65にある．
もしBrefが存在するなら，式2.2.11からBrefを使ってE obs （B obs ）中のE ind （B ind ）を見積もり，E cn（Bcn）と区別して評価することができる．そして残差のうち互いに相関のあるものをEcn（Bcn）とすればよい．Larsen（1996）やEgbert（1997）のコヒー
レントノイズ除去手法は基本的にこの考えに基づいている．一方，B refが存在しない場合は，比較する物がないのでEind（Bind）とEcn（Bcn）の統計的性質の違いに基づいてそれぞれを区別し評価することになり，一般に非常に困難である．おそらく主成分法などが使用できる手法だろう．
経験的に，E obs（B obs）の観測点に近いところで観測されたB refのほうが，Z bが安定に求まる傾向にある．特にデッドバンドのような短周期成分では参照点が遠いと非常に難しい．しかし，淡路島の場合，周辺部は都市化が進んでおり，近隣部の磁場観測点
は人工ノイズから逃げられない．そのため，B refとして数百km離れた柿岡観測点や鹿屋観測点のデー
タを使わざるを得ず，図2.1.3に示されたように，周期40秒より短周期側では淡路島のデータと有意な相関がない．そこで，本節では，柿岡や鹿屋のデータと淡路島のデータとの間に式2.2.10の関係が成立する長周期側に注目して解析を行う．
もしEcn（Bcn）に対する参照データが存在するな
ら，観測値からEcn（Bcn）を直接取り出すことも可能なはずである．実際，Ecn（Bcn）の原因が特定の路線の電車だとして，変電所の電圧のデータを用いた例や，電車の漏洩電流を矩型波モデルで表現した例がある．しかし，淡路島のデータの場合，路線も本数も膨大で，かつ，個々の電車に対応する変化が識別できないため矩形波モデルもたてにくい．ノイ
ズ源と思われる阪神地区から海を隔てているため信号が鈍る可能性があること，阪神地区には数多くの路線があり電車の本数も多いので個々の電車というより平均的にノイズが存在しているように見える可能性があることが指摘されている．
以上をまとめると，電場の観測値Eobsに含まれるEtc以外の成分のうち，Esqtは特定周期成分を除去することによって分離し，Eindは柿岡の磁場Brefと相関のある成分を除去することによって分離する．E cnを見積もるためには，先に大谷の磁場観測値Bobsから柿岡の磁場B refと相関のある成分B indを取り除き
B cnを取り出しておく．このとき，B obsにも地殻活動による成分Btcが含まれているとBcnと分離ができないが，電磁場変動発生メカニズムと観測状況の考察よりBtc=0が成立すると想定している．EcnはBcnと相関のある成分として見積もることができる．EindとE cn はルーチン的に計算できるようにするため，
ディジタルフィルターを開発する．
2．2．1．3 フィルター係数の計算手法
時系列x(t)と参照時系列y(t)にFFTを適用して求めたX(ω)とY(ω)が，変換関数P(ω)を用いて線形関係
X(ω)=P(ω)Y(ω)（2.2.12）と表される時，式2.2.12をインバースFFT変換すると，
Nx(t)=∑ p(i )y(t+i∆t)
i=–N（2.2.13）と表される．このようなフィルター係数pを計算することにより，yからxを予測することは，よく行われてきた．例えば，BAYTAP-Gの参照データと相関する成分の計算や，確立差分法（Mori, 1987），
徳本（1998）が用いた回帰分析などは直接xとyからpを最小自乗的に求めている．しかし，これらの方法では有限データ長の離散化データを使うことか
らギブス現象がおこる可能性がある．
周波数領域のPから時間領域のpを求めることも可能である．Egber t（1992）は，周波数領域から時間領域の変換に単純なインバースFFTを使うと，変換関数の計算値が存在しない周波数で振動するような解が得られることを示し，周波数領域での平滑化を提案した．また，電磁誘導における電場と磁場のように因果関係が成立している物理現象で観測値が参照値の過去の値だけに依存している場合でも，データの有限性や離散性から因果関係が破たんしてしまうこと，そのため最適なフィルター係数としては現在・過去・未来に依存するとしたほうがよいことも報告している．そこで，Egber t（1992）の方法を用いて，平滑化のペナルティ付きで，観測されたK個の変換関数Pを最もよく説明するように，フィルター係数pを計算することにする．具体的には，pの周波数領域での計算値 P̂とすると，
Kˆ k ))2/ε(ωk )2χ2 =∑(P(ωk )–P(ω
k=1（2.2.14）のようにχ2を計算し，χ2が自由度2Kになるpを求める．εはPの誤差である．平滑化のためのペナル
ティ関数fは，P̂が周波数領域で直線的に変化するな
ら一次微分の和，微分が滑らかに変化するなら2次微分の和というように変えることができる．この場合は，図2.2.1からもわかるようにlog周波数領域で

66

図2.2.1 柿岡の磁場X，Y成分に対する大谷の磁場（a）X成分と（b）Y成分の変換関数

滑らかであることが望ましいので，2ˆ

f = f ( d P2 )dw

（2.2.15）あるいはˆ

f = f ( dP )dw

（2.2.16）

ただし
w=log10ωと定義する．ペナルティの寄与の大きさを決めるた
めには，fに対するラグランジェパラメータαを導入する．実際の計算では，フィルターの長さ2Nを固定して，αを0.1-0.0001まで動かしながχ2が2Kより数％程度大きくなるαを求め，そのαでの計算値
pを長さ2Nでの最適値とした．理想的なχ2＝2Kとしなかったのはオーバーフィットを防ぐためである．次に，Nを変えて，上記の計算を繰り返す．数
種類のNについて最適フィルターが求まったら，周
波数領域，時間領域でそれぞれを目で比較し，周波数領域でのフィットがよく時間領域での振動が少な
いものの中で，最もNが小さいものを最終的な最適フィルターとした．Nの選択はAIC，αの選択も含
めたNの選択はABIC（Akaike, 1980）などの尺度もあるが，ここでは使わなかった．これは，それらの尺度がややNを小さめにとる傾向があり，時間領域での振動を選択の基準に反映しづらいためである．
2．2．1．4 大谷の磁場10秒値のモデル
短周期（40-1000秒）の人工ノイズを精度よく除去することを目的として，10秒サンプリングデータを使ってノイズ除去を試みる．長周期の変動に対しては，1分値を使ったモデルをたてる予定である．
ここでは，手法を適用する例として柿岡と大谷の
磁場3成分，野島-大谷測線の地電位差を選び，1998
年1月1日から3月31日までのデータを使った．観測された1秒値に対し，位相回転がないタイプのロー
パスフィルターを作用させて周期20秒以下の成分を取り除き，10秒値にリサンプリングした．目的の周期帯にはSqと潮汐が含まれないので，除去すべき
成分は外部磁場変動に伴う成分と人工ノイズ成分である．
まず，柿岡の磁場X，Y成分に対する大谷の磁場3成分の変換関数を計算する．計算には，第2．1．1．
2節で用いたのと同じロバスト手法を取り入れた変換関数演算プログラムを使用した．リモートレファ
レンスなしでπ prolate windowを適用しており，エラーはジャックナイフ法により算出した．
図2.2.1に，大谷の磁場水平2成分の変換関数を示す．Z成分は，短周期帯で変換関数が安定しなかったので図2.2.1には示さなかった．5.7x10-4 - 4.5x10-2
Hzの計算は，日中の強力な人工ノイズを避けるた
め深夜2:00 - 4:00の2時間のデータのみを用いて行われた．6.7x10-5 - 4.7x10-4Hzは，人工ノイズがかなり弱まっている周波数帯なので，すべてのデータを用いて計算した．深夜のデータから求めた変換関数と全日データから求めた変換関数は，滑らかに接続している．高周波帯では，ADコンバーターの機器応答の影響もあり，変換関数が大きく変化している．
しかし，低周波帯では安定し，周波数にかかわらず同じ成分同士ではほぼ大きさ1，違う成分同士では0.2前後となっている．

67

次に，図2.2.1の変換関数を時間領域に変換し，
N1N2i=–N 1i=–N 2N3N4i=–N 3i=–N 4

kak
kak
B otn
x (t)=∑ px(i )B x (t+i∆t)+∑ py(i )B y (t+i∆t)

kak
kak
B otn
y (t)=∑ qx(i )B x (t+i∆t)+∑ qy(i )B y (t+i∆t)（2.2.17）

図2.2.2 大谷の磁場X成分の柿岡の磁場X成分に対するイン
パルスレスポンスの長さ300，400，500，600のものを示す．（a）全体図，（b）時刻0付近の拡大図，（c）フィルター端の拡大図のフィルター係数p x，p y，q x，q yを計算する．N 1，
N2，N3，N4はそれぞれのフィルターの長さである．
図2.2.2（a）-（c）に，大谷の磁場X成分と柿岡の磁場X成分の間のフィルターpxを周波数領域と時間領域で示した．2N1として，300，400，500，600を選び，それぞれのN1での最適なpxが示されている．
2K=74であるのに対し，χ 2は75.01〜75.12であり，周波数領域でのフィッティングの良さはどれでも同じである．pxは時刻0周辺で大きな値を持つが，その他の時刻では急速に0に近づく形をしている．どのフィルター長でも時刻0周辺の値はほぼ同じだが，
時刻N 1周辺を拡大してみると200付近で振幅が0になり，200以上ではほぼ0が続いていること，がわかる．そこで，この場合は2N1=400を最適フィルターとした．
図2.2.3にp y ，q x ，q y を示した．いずれも，長さ400が最適と評価された．Y成分どうしでは最大値が1に近く時刻0付近に極在する形になり，図2.2.2の
X成分どうしの場合とよく似ている．違う成分どう
しでは，最大値が0.1-0.2程度であった．pyとqxの計算では，最長の2，3周期の変換関数を入れるとフィットがうまくいかなかったため，計算から取り除いている．他の周期に比べると不確かさの大きい計算値であったため，ややばらついていて，フィッティ
ングが不安定になったと思われる．
px，py，qx，qyを用いて1998年1月16日の柿岡の磁場X，Y成分から大谷の磁場X，Y成分を予測し，観測されたものと比較した（図2.2.4，5）．両成分とも約30nTの変動が観測されていたが，午前（23:00 4:00UT）と深夜（17:00 - 19:00UT）を除いて，ほとんどの長周期変化が外部磁場変動に伴う成分で説明

図2.2.3 大谷の磁場の柿岡の磁場に対するインパルスレスポ
ンス（全長400）．大谷のX成分と柿岡のY成分（太い点線），大谷のY成分と柿岡のX成分（太い実線），
大谷のY成分と柿岡のY成分（点線）の組み合わせのレスポンスの時刻0付近を示す．

68

できることがわかる．午前の長周期変化はピーク時
間や形は一定ではないもののほとんどの日に見られることから，柿岡と大谷で日変化のタイプが違うことからくる見積もり誤差と思われる．主に数時間以下の周期成分に注目した解析なので日変化に対応できないことはあらかじめ想定できた．しかし，柿岡
より鹿屋のほうが大谷の日変化に近い変動をするので，同じ解析でも鹿屋の磁場を参照した方がより良い結果を得られる可能性は残されている．
大谷の磁場Z成分では柿岡のX，Y成分との変換関数が短周期で定まらず，フィルター係数が計算できていなかった．これは，外部磁場変動に伴う成分よりも人工ノイズのほうが卓越したためと考えられる．平均的な外部磁場変動を考えると，X，Y成分
はソースの成分と誘導された成分の和でありどんな地下構造の場合でも存在するが，Z成分は地下構造の不均質によって誘導された成分のみを含む．そのため，第一次近似的に，X，Y成分で外部磁場変動に伴う成分が人工ノイズより強度があったとしても，Z成分では人工ノイズのほうが卓越することは不思議ではない．そこで，参照データに，柿岡の磁場X，Y成分だけでなく外部磁場変動に伴う成分を取り除いた大谷のX成分を加え，人工ノイズにも対応できるようにして，変換関数を計算した（図2.2.6）．
X，Y成分の変換関数（図2.2.1）と同じ期間のデー
タを用い，7.4x10-3 - 4.0x10-2Hzは深夜の値を使って計算し，9.5x10-5 - 5.7x10-3Hzは昼間のデータを用いた．大谷の補正済みX成分を加えたことでコヒーレンシーが有意に上がり，高周波側の計算が安定した．
図2.2.1に比べると，係数値が小さいこと，短周期でばたつきが大きいこと，低周波でも定数に漸近せず変化が続いていることが特徴である．
図2.2.6の柿岡X，Y成分に対する変換関数を使って図2.2.2-3と同様の計算を行い，得られた最適なフィルター係数を図2.2.7に示す．X，Y成分に対するインパルスレスポンス共に，長さ200の場合を最適と判断した．係数の最大値は，両成分共に同程度であり，X成分どうし，Y成分どうしの場合より小さく，異成分どうしの場合の1.5-3倍となった．また，
時刻0以外の寄与率が比較的大きい．

図2.2.4 1998年1月16日の大谷の磁場X成分について，（a）観測値（灰色の点線）と柿岡のX，Y成分から予測された値（実線），（b）残差69

図2.2.5 1998年1月16日の大谷の磁場Y成分について，（a）観測値（灰色の点線）と柿岡のX，Y成分から予測された値（実線），（b）残差図2.2.6 柿岡の磁場X，Y成分に対する大谷の磁場Z成分の変換関数

図2.2.7 大谷の磁場Z成分の柿岡の磁場X成分（実線），Y成
分（点線）に対するインパルスレスポンス（全長
200）．

70

図2.2.8 1998年1月16日の大谷の磁場Z成分について，（a）観測値（灰色の点線）と柿岡のX，Y成分から予測された値（実線），（b）残差図2.2.9 1998年1月から3月までの昼間 （9:00-24:00LT）における大谷の磁場3成分から柿岡の磁場と相関のある部分を除去した後のパワースペクトル71

図2.2.7のインパルスレスポンスを使って，柿岡のX，Y成分から予測した1998年1月16日の大谷磁場Z成分と観測値，残差を図2.2.8に示す．17:00UTごろの変化などはうまく説明できているが，全般に誘導成分の寄与が低いようである．これは，柿岡の磁場
X，Y成分だけを参照した解析では変換関数が安定しなかった理由として考えた仮説と調和的である．
午前の日変化と思われる変動は，やはりうまく説明できていない．
柿岡の磁場に相関のある成分を除去した残差の性質を調べるために，大谷の残差3成分のパワースペ
クトルを計算した（図2.2.9）．用いたデータは図2.1.10と同じである．外部磁場変動に伴う成分除去の影響は，X，Y成分の周期300秒以上に顕著だが，Z成分にはほとんど見られていない．X，Y成分の長周期のパワーが減ったことで，生データではわからなかったおよそ900秒と600秒のスペクトルピークが確認できるようになった．また，残差のX，Y成分のパワースペクトル構造はZ成分のもの，野島-大谷
測線の電位差のものと良く似ている．これらは，本節で適用された手法によって，主に短周期帯に見られる大谷の磁場中の人工ノイズ成分を取り出せたことを示唆していると考えられる．
2．2．1．5 野島−大谷測線の地電位差10秒値のモ
デル
ここでは，野島−大谷測線の地電位差から，まず人工ノイズ成分を取り除くことを試みる．外部磁場変動に伴う成分は対象とする周期帯では相対的に小さいので，補助的に扱うことにした．柿岡の磁場X，
Y成分に加え，外部磁場変動に伴う成分を取り除いた大谷の補正済み磁場3成分を参照データとして，
人工ノイズ用のフィルターを計算することを考える．具体的には，最適な参照データの組み合わせを

図2.2.10

見つけること，昼間と夜間のそれぞれの人工ノイズに対する応答を求めること，がこの節の目標である．
外部磁場変動に伴う成分と人工ノイズを表現する
ためには，柿岡の磁場X，Y成分と大谷の補正済み
磁場X（Z），Y成分の計4成分が最低限必要であると予想される．しかし，応答関数を計算するために用いたプログラム上の制限から参照データとして3成分しか使うことができなかったため，まずは最も振幅の大きいと考えられる人工ノイズに対する応答が精度良く求まる成分の組み合わせを探すことにした．比較したのは，以下に示すような6つの組み合わせである．
Case 1: 大谷の補正済み磁場のX，Y成分
Case 2: 大谷の補正済み磁場のX，Y成分と柿岡の
磁場X成分
Case 3: 大谷の補正済み磁場のX，Y成分と柿岡の
磁場Y成分
Case 4: 柿岡の磁場X，Y成分と大谷の補正済み磁場のX成分
Case 5: 柿岡の磁場X，Y成分と大谷の補正済み磁場のY成分
Case 6: 柿岡の磁場X，Y成分と大谷の補正済み磁場のZ成分
1998年1月1日から3月31日の10秒値を使って計算したCase 1-5の応答関数のうち大谷の補正済み磁場のX，Y成分に対するものを図2.2.10，11に示す．
Case 6はCase 1とすべての応答関数についてよく似ているため割愛した．図2.2.10は昼間の応答関数であり，顕著なラインスペクトルが見られる周波数帯では応答関数が乱されるので，その周波数帯の応答関数は取り除いてある．深夜の時間帯にも図2.1.13に見られるようにスペクトルの緩やかなピークが存

LT9:00-24:00における大谷の補正済み磁場（a）X，（b）Y成分に対する野島-大谷測線の地電位差の
応答関数

72

在し応答関数が乱されたため，やはり図2.2.11から取り除いた．
図2.2.10-11を見ると，Case 1-3の応答関数は互いに良く似ており差はほぼエラーバーの範囲内におさまっているのに対し，Case 4あるいは5の応答関数は有意な違いを示している．このことから，地電位差変動には人工ノイズのほうが外部磁場変動より大きな影響を与えることがわかる．また，地電位差の測線がほぼ南北方向であるにもかかわらず，大谷の補正済み磁場のX成分に対する応答関数のほうがY成分にたいするものより安定して求まっていて，
Case 4とCase 1-3の違いよりCase 5とCase 1-3の違いのほうが大きいことから，人工ノイズのうちでも磁場の北向き成分に現れる変動のほうが東向き成分に現れるものより寄与が大きいことが示唆される．これらは第2．1．2節で明らかにされた人工ノイズの偏向性と調和的である．ラインスペクトルによって100-1000秒の応答関数が乱されているが，X成分よ
りY成分のほうが大きく影響を受けているようである．
以上の結果から，Case 4，5は大谷補正磁場の応答関数の対象外とする．Case 1-3の応答関数を比べると，Case 1よりCase 2-3のほうが長周期帯でより安定しており，Case 2と3では信号の強さを反映してCase 2が広い周期帯にわたって安定しているので，Case 2の大谷磁場の応答関数を使用することにする．昼間と夜間のCase 2の大谷補正済み磁場X，
Y成分の応答関数を図2.2.12に示す．
図2.2.12から，人工ノイズの応答関数は約1000秒
以下の周期帯で，実部が周波数によらず一定，虚部が0というDC的な振る舞いをすることがわかる．昼
間も夜間もDC的であることは変わらないが，X成
分では昼間のほうが抵抗が大きく，Y成分では夜のほうがやや強い周波数依存性を示している．また，昼間の応答関数は15時間の連続データのスタッキン
グであり周期10000秒程度まで安定して計算できるが，深夜の応答関数は2時間の連続データのスタッ
キングであるため周期1000秒以上の計算が非常に不安定であった．そのため，昼間の応答関数ではハイパスフィルター型の周波数分布が顕著に見られるが，深夜の応答関数では実部が0に収束していく長周期側の部分が除かれた分布になっている．
図2.2.12の応答関数を説明するインパルスレスポ
ンスを，昼間と深夜についてそれぞれ計算した．レ
スポンス計算では，平滑化は式2.2.16のlogの一階微分を採用し，フィルター長さを150から450まで50ず
つ増やしながら計算を行った．X成分では長さ300の場合が，Y成分では長さ400の場合が，振動の少なさから最適とみなされた（図2.2.13a）．深夜のレスポンスは，長さ150-400について計算が行われ，X，
Y成分でそれぞれ長さ250，350の場合を最適レスポンスとした（図2.2.13b）．どのレスポンスも，DC成分を表現する時刻0で図2.2.12の短周期帯の実部と近い値を示し，その他の時刻では振幅が100分の1程度になって時刻が0から離れるにつれ値が0に近付いている．
図2.2.13（a）に示された昼間のインパルスレスポ
ンスをフィルターとして用いて，大谷の補正済み磁場X，Y成分と相関のある電位差変動の予測値を求め，実測値，残差と共に図2.2.14に示した．用いられた年月日は前節と同じ1998年1月16日で，UT0:0016:30，19:30-24:00について予測値を計算している．
大谷の補正済み磁場は，簡単のためラインスペクト
ルを含んだ値（図2.2.4−5の残差）を用いた．本来
ならば，ラインスペクトルに対する応答は図2.2.12−13とは違う値であるので，ラインスペクト
ル分の寄与は別に評価し，図2.2.4−5の残差からラインスペクトルを取り除いた値に図2.2.13（a）のフィルターを作用させる必要がある．しかし，第2．1．
1節で用いたような見積もり法だとLT依存性をうまく取り入れられないことや，見積もりに用いた時刻とそれ以外の時刻との間でラインスペクトル成分を取り除いたことによるギャップを生じることなどか
ら，人工ノイズ成分を計算する最初の試みとしては適用しなかった．そのため，図2.2.14の残差には本来と違うフィルターを作用させたことによるライン
スペクトル成分のずれが含まれている．
図2.2.14を見ると，大谷の補正済み磁場X成分から予測された電位差は振幅400-500mVで短周期成分が卓越している．Y成分から予測された電位差は，短周期成分の振幅は100mV未満でX成分に比べると小さいが，LTの午前には日変化に対応すると思われる緩やかな変化が見られる．どちらの予測値も，深夜の時間帯近辺では変動の振幅が小さくなっていて実測値の傾向と対応しているが，深夜の時間帯に地電位差のベースが著しく0に近づく傾向は反映されていない．残差をみると，短周期成分の振幅は100mV程度まで小さくなっており，矩形的な変化やベースの変化など長周期の変動が実測値よりも目立
つ．
大谷の補正済み磁場に相関のある成分を取り除いた操作について影響を調べるため，30日間の残差のパワースペクトルを計算し観測値のものと比較した（図2.2.15）．予測値は2000秒以下の周期を持つこと，
残差の連続スペクトルのパワーが実測の10分の1程度であることから振幅のおよそ90%が予測されたこ

73

図2.2.11

LT2:00-4:00における大谷の補正済み磁場（a）X，（b）Y成分に対する野島-大谷測線の地電位差の応答関数

図2.2.12

昼間と深夜の大谷の補正済み磁場（a）X，（b）Y成分に対する野島-大谷測線の地電位差の応答関数

図2.2.13

大谷の補正済み磁場X成分（実線），Y成分（点線）のための（a）昼間と（b）深夜のインパルスレス
ポンス

74

図2.2.14 1998年1月16日の野島-大谷測線の地電位差の（a）実測値，昼間のインパルスレス
ポンスを用いて予測された大谷の補正済み磁場（b）X成分，（c）Y成分と相関する成分，（d）残差

75

図2.2.15 1998年1-3月の野島-大谷測線の地電位差の実測値（点線）と大谷の補正済み磁場の影響を除いた残差（実線）のパワースペクトルとがわかる．しかし，残差の連続スペクトルは依然
周期200秒付近に自然ソースにないピークを示して
おり，人工ノイズ成分が取り残されていることが示唆される．ラインスペクトルに関しては，当初の予想通り残差中に顕著に存在するが，実測に比べてパワーが減ったもの，あまり変わらないもの，新しく現れたもの，というようにスペクトルによって異なる振る舞いをすることがわかった．
残差中に含まれるラインスペクトルの影響をおおまかに見積もるため，第2．1．1節で用いたDFT法によって平均的なラインスペクトル成分を試算した．1998年1月16日のUT0:00-15:00の結果を図2.2.16（a）に示す．ラインスペクトル成分は振幅約100mV
で，ビートが見られUT3:00と10:00に振幅が最大になる．図2.2.16（b）は，残差からラインスペクト
ル成分を取り除いたもののUT11:00-14:00を拡大して示したものである．矩形変化の他に2−3分周期のピークが見られ，0，15，30，45分などの時刻にピークが現れる傾向があるので，電車の発着などの人
工的な原因による変動と考えられる．外部磁場変動による地電位差変動を概算すると図2.2.16（b）より1桁程度振幅が小さくなり，図2.2.16（b）の変動の主要な原因とはなり得ないようである．
図2.2.14−16の結果から，ラインスペクトル成分を除いた最終残差に，まだ人工ノイズが残っていることが推測されるが，その原因としては次の2つが挙げられる．第1に図2.2.12−13のレスポンスに見積もり誤差がある場合，第2に予測された人工ノイズ
はsin波の重ね合わせで表現した平均的なものなので平均からのずれが残差中に残る場合である．第1の見積もり誤差に関しては，図2.2.12−16の一連の手続きを最終残差が小さくなるように繰り返すイタ
レーション化によって，精度を上げることができると予想できる．大谷補正済み磁場のX成分の応答関数は安定しているので大きく変わらないであろうが，Y成分は改善される可能性がある．イタレーシ
ョン化する場合は，観測網全体にデータ処理を広げるためにも，一連の手続きの自動処理化が課題になる．第2の人工ノイズの平均からのずれの見積もりについては，いくつか方法が考えられる．ずれが有限個の周期の波の重ね合わせで表現可能ならば，最も簡単な方法は，図2.2.12−13のレスポンスに時間依存性を持たせる方法である（例えば，Widrow etal.（1975）のアダプティブフィルター法）．レスポ
ンスを使わずに時間依存する成分を見積もる主成分法も候補であろう．一方，平均からのずれが，電車の発着のように本質的にパルス的あるいは矩形的な変化を有限個の波の重ね合わせで表現したことから生じたのなら，矩形波を基底関数に用いたウォルシ
ュ解析を適用することも考えられる．その延長には，
時間依存する基底関数を用いるウェーブレット解析も候補に挙がってくるだろう．
電車などの人工ノイズ以外の理由によるノイズ，
例えば電極の変化などが最終残差に含まれている可76

図2.2.16 （a）1998年1月16日昼間の野島-大谷測線の地電位差から大谷の補正済み磁場の影響を除いたもの（点線），ラインスペクトル成分（一点鎖線），残差（実線），（b） UT11:00-14:00の残差能性もある．それらについては，大部分の測線で複数の電極を用いた観測が行われているので，別の電極ペアの測定値との比較をすることで区別することができると考えられる．いずれにせよ，観測網の他のデータと交互に参照しあうようなアルゴリズムを組めば，ノイズの見積はこの節に示したものより改善されるはずである．
次に，図2.2.13（b）に示された深夜のインパル
スレスポンスをフィルターとして用いて大谷の補正済み磁場X，Y成分と相関のある電位差変動を求めた．図2.2.17に1998年1月16日のUT15:00-20:00の予測値と実測値，残差を示す．短周期の変動について
は予測と実測の傾向が一致して残差の振幅が小さくなっており，昼間のインパルスレスポンスを使って深夜の予測をした場合よりも実測に近い予測値が得られた．しかし，予測値には実測値にない大きなト
レンドがあり，それが残差のトレンドとなっている．
このトレンドのため，昼間のインパルスレスポンス
を用いて計算した残差とUT15:00付近や20:00付近でうまくつながらない．これは，深夜の応答関数（図2.2.12）に長周期成分がないためにそれを変換した
深夜のインパルスレスポンスでは長周期成分のコントロールが効かなくなっており，大谷の補正済み磁場の長周期のトレンドが予測値にそのまま反映されてしまうことを示している．
応答関数の計算に使えるノイズの小さい深夜の時
間帯は，LT2:00-4:00（UT17:00-19:00）より著しく延びることはなく，観測値から長周期成分を計算できる可能性は低い．しかし，図2.2.12を見ると，周期1000秒以上で昼間と深夜の応答関数はよく一致している．そこで，試みとして，深夜の応答関数が計算できなかった長周期の部分に昼間の応答関数を補間し，合成した深夜の応答関数を使ってインパルス
レスポンスを計算してみた（図2.2.18）．合成した応答関数から求められたインパルスレスポンスを用いて，大谷の補正済み磁場X，Y成分と相関のある電位差変動を求め実測値と比較したものを図2.2.19に示す．予測値は短周期変動が卓越しており，残差には実測値のトレンドがそのまま残る形になっている．長周期部分に急速に0になる昼間の値を使った

77

図2.2.17

1998年1月16日の野島-大谷測線の地電位差の（a）実測値（点線），深夜のインパルスレスポンスを用いて予測された大谷の補正済み磁場X成分，Y成分と相関する成分（実線，薄い実線）， （b）残差ことで，本来ならば予測されるべき長周期成分をも排除してしまった可能性はあるが，少なくとも不要
な長周期変動を取り込んでしまう可能性はなくなったと考えられる．
図2.2.19の残差を見ると，地電位差は深夜の時間
帯UT15:00-19:30に特徴的な変化をしている．平均
的なベース値がそれ以外の時間帯の約-350mVから約0mVへ変わっており，16:30で山，17:30で谷，
18:30で山というように2つの正のピークがある．
1998年の1-3月について調べてみたところ，この変動はほとんど同じ形状で日によらず存在する．図2.2.20に1月11-20日の10日間についてUT14:00-21:00の残差を示す．合成した深夜のレスポンスを用いて大谷補正済み磁場と相関する成分の除去を行い，特徴的な変化をわかりやすくするためカットオフ周期
1200秒のローパスフィルターを作用させている．深夜のベース値の変化と，2つの山とその間の谷が共通していることがわかる．この変動の原因はよくわからないが，ベースの変化はこの時間帯になんらかの都市活動（電車など）が止まっていることを示し

図2.2.18 大谷の補正済み磁場X成分（実線），Y成分（点線）のための深夜のインパルスレスポンス．深夜のイ
ンピーダンスに昼間のインピーダンスの長周期部分を組み合わせた合成インピーダンスに対して，最適なインパルスレスポンスを計算した．

78

図2.2.19 1998年1月16日の野島-大谷測線の地電位差の（a）実測値（点線），合成した深夜のインパルスレス
ポンスを用いて予測された大谷の補正済み磁場X成分，Y成分と相関する成分（実線，薄い実線），（b）残差

図2.2.20 1998年1月11-20日の野島-大谷測線の地電位差から深夜の予測値を除いた残差の長周期成分．ローパ
スフィルターのカットオフ周波数は1200秒．

ていると思われる．2つの山とその間の谷について
も，自然の変化としては日ごとの違いが小さすぎる
ので，やはり人工的なものと思われる．これらの谷の部分が徳本（1999）が示した平均値に対応するのであろう．ここでは短周期のノイズ成分を磁場から予測して取り除いているため，徳本（1999）より長周期的な傾向が鮮明に現れた．この特徴的な変化をノイズ除去に用いるためには，それぞれの日のベース値のずれなどを補正して，真の平均値を求める手法を考える必要がある．

79

2．2．2 主成分分析法を用いたノイズ除去
2．1．2でわかったように，淡路島で我々が観測している地域では地電位差変動，地磁気変動の卓越
ノイズの空間分布のスケールは大きい．それに対して淡路島内の断層運動に伴う場合のように，地殻活動に伴う電磁気変動源が観測エリア内の浅部にあるような場合には，観測されるシグナルの空間分布のスケールはもっと小さい（エリア内でのシグナル分
布のコントラストが強い）ことが予想される．従ってノイズを除去する手段として予想されるシグナルとの空間スケールの違いを利用することが考えられる．ここでは，主成分分析法を利用したノイズ除去について検討してみる．
2．1．2で行った主成分分析によって地電位差変動の広域的な変化を抽出することができた．しかもその変動成分の振幅は全体に占める割合が大きいため，主成分分析で得られた上位の主成分を取り除く
ことによって，小さな空間スケールの地電位分布の
時間変動のみを取り出すことができる．ここでは地電位差データに擬似的なシグナルを混入させて，それがこの方法でピックアップできるかどうか試験を行った結果を報告する．まず方法について述べる．
主成分固有ベクトルbkが求まると，ある地電位差時
系列データẽ(tn)の中に含まれる第k主成分の振幅akを式2.1.7に従って求めることができる．従って，
Kẽ′(tn)=ẽ(tn )–∑ak(t n )bk
k =1（2.2.18）で得られる時系列データは，元のデータから第1〜
第K主成分までを取り除いたものである．
図2.2.21に試験に用いた1999年2月18日02:00〜
02:59UTの地電位差毎秒値データを示す．ここで示したデータは2．1．2で主成分分析を実行する際に行ったように，全21点の平均からの偏差を取り出し，
周期100秒のハイパスフィルターをかけている．この図で明らかなように，ノイズの振幅は大きい地点
では200mVを越えていることがわかる．また02:46
頃から始まる激しい変動はこの時に発生したSSCに伴う変化である．ここで示した地点4（野島）のデータには測定したデータに対して地殻活動によるシグナルに見立てて，擬似的なパルスデータを加算している．その波形は図2.2.22に示すような三角波で，振幅が20mV，継続時間が10秒である．02:10と02:30の2ヶ所にこの三角波を与えているが，当然のことながら図2.2.21ではノイズに紛れて識別することができない．
2．1．2では1998年1月15日の1日間の地電流毎秒値を使って，主成分分析を行った．その結果によると少なくとも第6主成分までは，上位の主成分ほどよ り 広 域 的 な 変 動 を 表 し て い る よ う で あ る（ 図2.1.17）．従ってそれ以降の主成分についても上位の主成分ほど広域的な変動を表しているものと仮定
し，図2.2.21のデータの中の上位9主成分の変動成分を計算し，元データから取り除いた結果が図2.2.23
である．この図は図2.2.21に比べて縦軸の表示スケールを10倍に拡大している．これを見ると地点毎に効果にばらつきはあるものの，変動の振幅は非常に小さくなっていることが明らかであり，ノイズの大部分を取り除くことができたと考えられる．ところ
でSSCに伴う変動についても振幅が10分の1以下となっているが，これはSSCという地球スケールの地磁気変動によって誘導される地電流変動の大部分は空間スケールが大きいため，広域的なノイズと同様に除去されたものと思われる．なおSSCに対応する
時刻で除去しきれていない量は地電流変化の局地的な変動成分と考えられる．ところで02:10と02:30に
地点4，5，8にパルス状の変動が現れている．これが前に述べた擬似的なシグナルに対応している．この図でわかるように第9主成分までを取り除くこと
によりこのシグナルを十分に識別できる程度にノイズが取り除かれている．このうち02:10分を中心に2分間の変動を拡大した図2.2.24を見ると，地点4（及び5）の三角波が再現できていることがわかる．な
おシグナルを与えたのは地点4だけであったのに対し，地点5，8のデータにもパルス状変動が現れているが，これはそもそも主成分分析を適用するにあた
り式2.1.1によって全地点の平均からの偏差を解析に使用していることから，地点4と逆向きの変化が他の地点にも現れたものである．なぜ地点5，8だけにこのような変動が現れているのかを明確に説明するのは難しいが，地点4，5，8は互いに近接した測点であり，それらが逆方向に変化することは空間スケールの小さな高次の主成分に属する現象として残差に含まれたと解釈できるのかも知れない．
ここでもう一点特筆すべきことは，ここで試験対象としたデータと主成分の基底を算出するのに使用したデータとが全く異なる時期（1年以上）であるということである．シグナルの有無を定常業務として監視する場合には簡便な方法が望ましいが，今の方法ではいったん基底を計算しておけば，それを使って長期間のデータに適用できるということは実際の業務に取り入れやすいということが言える．
ところでこの手法にはまだ改良すべき点がある．
図2.2.23で示されている通り，ノイズの除去の程度
は観測点毎に異なっており，第9主成分まで取り除いても明らかにノイズとわかる変動が残っている場

80

図2.2.21 1999年2月18日02:00〜02:59UTの地電位差毎秒値．縦軸は地点番号，横軸は分．時定数100秒のフィルターで長周期成分は除いている．また1秒毎に全地点のデータの平均からの偏差を表示し
ている．

81

図2.2.22 テストに用いた疑似シグナルの波形．

図2.2.23 図2.2.21から第1〜10主成分を取り除いたもの．

82

図2.2.24

図2.2.23のデータのうち02:09〜02:10だけ表示したもの．

所もある．そのような点もさらに高次の成分まで取り除けばノイズは0に近づいていくが，そうするとその他の点のシグナルまで除いてしまう可能性がある．従って，第何主成分までを取り除くかということについて地点毎に最適な次数をきめることができ
れば，全地点でS/N比を向上することができると思われる．また別の問題として，ここではテストとした20mVよりもさらに小さなシグナルを検出する必要がある場合，あるいはシグナルがパルス状ではな
く時間スケールがもっと長いような場合には検出が難しくなる．現時点では断層活動の前兆としての地電流のシグナルがどの程度の振幅，時間スケール，波形であるか情報が十分ではないが，今後それが明らかとなってくれば，それに応じて主成分分析によるノイズ除去手法も改良する余地はでてくると思わ
れる．

83

参考文献（第2章）
Akaike, H., Information theory and an extension of the
maximum likelihood principle, 2nd Inter. Symp. on
Information Theory, edited by B. N. Petrov and F.
Csaki, 267-281, 1973.
Akaike, H., Likelihood and the Bayes procedure, in
Bayesian Statistics, ed. by J. M. Bernardo et al.,
University Press, Valencia, 143-166, 1980.
Chave, A.D. and D.J. Thomson, Some comments onmagnetotelluric response function estimation, J.
Geophys. Res., 94, 14215-14225, 1989.
Egbert, G. D., Noncausality of the discrete-timemagnetotelluric impulse response, Geophysics, 57,
1354-1358, 1992.
Egbert, G. D., Robust multiple-station magnetotelluric
data processing, Geophys. J. Int., 130, 475-496,
1997.
Egbert, G. D., M. Eisel, O. S. Boyd, and H. F. Morrison,
DC trains and Pc3s: Source effects in mid-latitude
geomagnetic transfer functions, Geophys. Res.
Lett., 27, 25-28, 2000.
Electromagnetic Research Group for the 1995 Hyogoken Nanbu Earthquake, Tectonoelectric signal
related with the occurrence of the 1995 Hyogo-ken
Nanbu Earthquake （M7.2） and preliminary
results of electromagnetic observation around the
focal area, J. Phys. Earth, 45, 91-104, 1997.
Fraser-Smith, A. C. and D. B. Coates, Large amplitude
ULF electromagnetic fields from BART, Radio
Science, 13, 661-668, 1978.
Johnston, M. J. S., Tectonomagnetism andtectonoelectricity, Reviews of Geophysics, 25, 983988, 1987.
Johnston, M. J. S., Review of electric and magnetic fields
accompanying seismic and volcanic activity,
Surveys in Geophysics, 18, 441-475, 1997.
石黒真木夫、佐藤忠弘、田村良明、大江昌嗣、地球
潮汐データ解析，統計数理研究所彙報，32, 7185，1984Larsen, J. C., Robust smooth magnetotelluric transfer
functions, Geophys. J. Int., 124, 801-819, 1996.
Linington, R. E., The magnetic dicturbances caused by
DC electric railways, Prospezioni Archeologiche, 9,9-20, 1974.
Mackie, R. L. and T. R. Madden, A magnetotelluric
survey around the Loma Prieta fault zone, EOS,
Trans. Am. geophys. Un., 73, 99, 1992.
三品正明, ネットワークＭＴデータに含まれる電車

漏洩電流ノイズの振幅変化について，Conductivity Anomaly研究会1995年論文集, 5964, 1995.
Mori, T., Variations in the geoelectric field with relation
to crustal conditions of the Earth, Geophys.
Magazine, 42, 41-104, 1987.
Mulargia, F. and P. Gasperini, Evaluating the statistical
validity beyond chance of VAN earthquake
precursors, Geophys. J. Int., 111, 32-44, 1992.
村上貴久, Network-MT法を用いた南部淡路島下の電気伝導度構造の推定, 修士論文, 京都大学，108
pp, 2000.
小嶋美都子, 伊豆大島で観測される地電位差異常変動の原因についての一考察, 地震, 44, 177-184,
1991．
小嶋美都子,地震予知―地電流による可能性, 地震ジャーナル, 14, 20-28, 1992．
小嶋美都子，徳本哲男，山田雄二，庄司哲也，田中智巳，地磁気・地電流データのノイズ除去手法の開発−重価差・フィルター等を用いた方法−，平成8年度地磁気観測所技術検討会，調査研究中間報告，1996．
Park S. K., M. J. S. Johnston, T. R. Madden, F. D.
Morgan and H. F. Morrison, Electromagnetic
precursors to earthquakes in the VLF band: A
review of observations and mechanisms, Reviews ofGeophysics, 31, 117-132, 1993.
Szarka, L., Geophysical aspects of man-madeelectromagnetic noise in the Earth - A review,
Surveys in Geophysics, 9, 287-318, 1988.
Takayama, H., comment on "Evaluating the statistical
validity beyond chance of VAN earthquake
precursors" by Mulargia, F. and P. Gasperini,
Geophys. J. Int., 115, 1197-1198, 1993.
高山寛美，山田雄二，田中智巳，角村悟，地磁気・地電流による地殻活動の検出手法についての調査，−BAYTAP-G を用いた方法−，地磁気観測所技術報告，第109号，39，2000．
徳本哲男，高山寛美，山田雄二，田中智巳，小嶋美都子，地磁気・地電流データのノイズ除去手法の開発，−重価差・フィルター等を用いた方法−，地磁気観測所技術報告，第106号，82，
1998．
徳本哲男，高山寛美，山田雄二，田中智巳，地磁
気・地電流データのノイズ除去手法の開発（その2），−重価差・フィルター等を用いた方法−，
地磁気観測所技術報告，第107号，37，1999．
Thomson, D.J., Spectrum Estimation and Harmonic

84

Analysis, Proc. IEEE, 70, 1055-1096, 1982.
上田誠也, ギリシャの地震予知, 科学, 55, 180-184,
1985．
上田誠也, 地震予知に成功した国―ギリシャ, 地震,44, 391-405, 1991．
Varotsos, P. and K. Alexopoulos, Physical properties ofthe variations of the electric field of the earth
preceding earthquakes, I, Tectonophysics, 110, 7398, 1984a.
Varotsos, P. and K. Alexopoulos, Physical properties ofthe variations of the electric field of the earth
preceding earthquakes, II. Determination of
epicenter and magnitude, Tectonophysics, 110, 99125, 1984b.
Widrow, B., J. R. Glover, Jr., J. M. McCool, J. Kaunitz,
C. S. Williams, R. H. Hearn, J. R. Zeidler, E.

Dong, Jr. and R. C. Goodlin, Adaptive noise
canceling: Principles and applications, Proc. IEEE,63, 1692-1716, 1975.
山田雄二，徳本哲男，田中智巳，高山寛美，淡路島
における地電流データの短周期ノイズ特性について，地磁気観測所技術報告，第109号，96，
1998．
山田雄二，田中智巳，横山恵美，高山寛美，藤井郁子，角村悟，淡路島における地電流・地磁気ノ
イズの空間分布について，地磁気観測所技術報
告，第110号，26，2000．
Zonge, K. L. and L. J. Hughes, Controlled source audiofrequency magnetotellurics, in Electromagnetic
methods in applied geophysics, vol. 2, Ed. by M. N.
Nabighian, Society of Exploration Geophysicists,
Tulsa, Oklahoma, 713-810, 1991.