第3章

3．1

自動検測手法の研究＊

はじめに

近年、地震活動の早期把握の必要性から、震源決定の即時性および精度の向上が急務となっている。これらの要求に応えるべく、観測網の拡大・充実、観測方式の近代化等が実施され、テレメー

タ方式による集中観測が一般的に行われるようになってきた。このため、多量の地震波信号を自動的に、しかも高速・高精度で処理するための方式の開発が、地震計測の重要な課題となって，きた。
地震波信号の処理は、大別すると、

i）

地動に含まれる地震波信号の検出

ii）

P波、s波の初動時刻の決定の2段階に分かれる。i）の処理は、信号系列の中に地震波が含まれているか否かを判別するもの
で、ii）の処理は、i）の処理の後、震源計算のための観測情報を得るものである。震源位置の決定精度は、観測網の幾何学的形状と関係し複雑であるが、初動時刻の決定精度に依存することは言
うまでもない。このため、ii）の初動時刻の決定は、十分な配慮の下で行わなければならない。

従来、これらの処理は人間が行っており、永年の経験と慎重な配慮の下で作業を実施するという意味から、『験測』（例えば、P波験測、S波験測）と呼ばれてきたようである。しかしながら、験測者の経験や主観の相違によってその結果が多少異なることは避けられない。最近では、このような主観的方法によらず、統計的モデルに従って処理することが可能となってきた。この新しい方法では、統計的合理性に基づき、客観的な検査・検定をして測定を行うので、『検測』と称し、従来の『験
測』と区別して使うこととする。

地震波信号の自動処理の一つに、1960年代に地下核爆発の探知のために考察されたもので、規則的に配置した群列地震計による手法がある［Carpenter（1965）、Greenet．a1．（1965）、Carpon（1969）］。この処理方式は、多数の観測点からのコヒーレント（Coherent）な信号の重ね合わせによる雑音除去を基礎として組立てられている。群列地震計は、地下核爆発の探知だけでなく、遠地地
震の観測にも広く利用されている。

1970年代に入ると、地震予知研究に関連して、多くの微小地震観測網が展開されるようになった。

しかし、遠地地震の波形信号と異なり、観測網の近傍に発生する地震の波形信号は、同一地震でも＊横田

崇：地震火山研究部

一56

し気象研究所技術報告

第16号

1985

観測点ごとの地震波信号に類似性が見られないのが普通である。そのため、通常我々が対象とする
地震観測網における自動検測処理では、群列式地震計のそれとは事情を異にし、インコヒーレント（incoherent）な地震波信号を取り扱わねばならない。したがって、コヒーレントな地震波信号を前提とした群列式地震計における自動検測処理の方式を、通常の観測網へ適用することはできない。

また、通常の観測網では、観測点の配置が幾何学的に不規則であり、群列式地震計とはその性格を異にする。このような事情から、それぞれの観測点からの信号を一つ一つ個別に検測処理する方式
6を採用する必要がある。さらに、地動（常時微動）そのものが変化するうえに、人工的振動も加わる場合もある。すなわち、地震波信号以外を通常、地動雑音と呼んでいるが、この地動雑音そのものが多様であり、自動検測方式の開発を困難にしている要因ともなっている。

このような自動処理方式の困難さにもかかわらず、従来より様々な方式が開発されてきた［例えば、渡辺ら（1975）、渡辺・黒磯（1977）、Stewart（1977）、Allen（1978）、Anderson（1978）、市川（1980）］。これらの方式の代表的なものは、ある種のフィルターあるいは特性関数からの出力の長期的平均値（10ng−termaverage）と短期的平均値（short−termaverage）との比較から地震波信号の検出および初動時刻を決定する方式である。しかし、これらの手法は経験的要素を多分に含み、必ずしも十分な合理性を持つ手法とはいえない。

地震波信号の検出手法とし七は、Freiberger（1963）の信号検出手法を基礎とし、Walsh関数を適用した手法がある［GoforthandHerrin（1981）］。これは、近似された尤度比の判別基準の設定

等に人間の介在を要するが、地震波信号の検出に関しては、最も検測的な手法と言えよう。従来の手法が以下で述べる検測的な手法からかけ離れていた要因としては、使用された計算機の能力からの制約も大きかった。

一方、地震波信号の自動処理において、最も重要な初動時刻の決定に、統計学的概念を導入した手法の開発が行われた。白井・徳弘（1979）は、ARモデルで地震波信号を表現し、ベイズの手法に基づき処理する方式を提案した。彼らの方法を適用して、浜口・鈴木（1979）、森田・浜口（1981）

は、P波、S波の初動時刻の決定を試みた。しかし、彼らの手法は、1時点毎にベイズの手法を適用

して信号か否かの判別を行う方式のため、不安定で尤度比の系列を平滑化して再評価する必要に迫られた。

横田ら（1981）は、OzakiandTong（1975）、Kitagawa

andAkaike（1978）らによる局所定常という概念を利用し、ある非定常区間を定常な2つの区間に分割する手法により、到着時刻の決定を行った。これより、数学的客観性を保ち、合理的に初動時刻の決定が実用化できるようになった。

以上、地震波信号の処理手法の変遷について述べたが、同様な問題は他の地球物理関係のデータについても言える。例えぽ、地殼変動関係データの処理では、地震発生前の異常の検出が地震予知のための重要課題となっている。Ishii（1976）は、チェビシェフ多項式を、石川・宮武（1978）はARモデルを、大内・高橋（1981）はIARモデルを、橋爪・三雲（1983）は多成分ARIMAモデル

一57一

気象研究所技術報告

第16号

1985

を用いて、それぞれ異常の検出を試みている。地殼変動関係のデータは、トレンド成分や季節変動成分等の多種類が混合したものであり、これら成分の分離も含めたより一般的な手法が開発されている［例えば、Ishiguro

et．a1．（1981）］。岡田・高橋（第6章）はIshiguroet．a1．（1981）の方法を用いて、検潮データから地殼変動の解析を行っている。また、検潮データに関しては、津波到着時
刻の決定問題などもある。

地震波信号、地殼変動関係の信号、検潮信号等の処理は、それぞれに含まれる信号を検出する点において、本質的に同じである。このようなことから本報告では、異常の検出という観点から処理の内容に一般性を持たせ、

i）信号検出ii）

信号出現時刻の決定のための自動検測手法について述べる。

本章は6節からなり、以下の各節の内容は次の通りである。
3．2節『信号検出』では、信号検出問題の基本的手法について述べる。

3．3節『雑音の中の信号検出』では、実際的な時系列データの信号検出手法を扱う。
3．4節『信号出現時刻の決定』では、横田らの手法を一般化した形で記述する。

3．5節『地震波信号における相の同定』では、特に地震波信号の際に問題となるP相、S相の同定法について記述する。

3．6節『おわりに』では、本報告における結果と将来の展望について述べる。

参考文献
浜口博之・鈴木次郎、1979：微小地震P波初動の自動検測とその信頼性について、自然災害特別研究成果、
A−54−2、62−83。

橋爪道郎・三雲健、1983：地殼変動データの多変量解析

ARIMAモデルによる予測の試み、地

震2、36、515−530。
市川政治、1980：地震記録自動処理装置、験震時報、44、55−74。

石川有三・宮武
31、隆、1978ニウィーナーフィルタの適用による地殻活動・地震活動の予測め試み、地震2、

73−86。

森田裕一・浜口博之、1981：2次元自己回帰過程によるS波初動の自動検測、地震2、、34、223−240。

大内2、

徹・高橋亨、1981：地震学で扱う非定常系列データの解析
34、予測と異常検出の試み、地震455−464。

白井克彦・徳弘一路、1979：地震波初動の自動検出、地震2、32、141−148。
横田

崇・周

勝奎・溝上

恵・中村

功、1981二地震波データの自動検測方式とオンライソ処理システム

における稼動実験、地震研究所彙報、55、449−484。

渡辺一郎・菅原正己・福井隆文・勝山ヨシ子、1975：微小地震の自動検出方法について（第3報）

一58

自動

気象研究所技術報告

検測の試み
渡辺

第16号

1985

、国立防災科学技術セγター研究報告、15、33−47。

晃・黒礎章夫、1977二阿武山地震観測所における微小地震の自動処理について、地震2、30、115−118。

Anderson，KR．，1978：Automatic

Seismic

Analysis

andanalysisofmicroearthquake

networkdata．Computer−Aided

Discrimination，ed．C．H．Chen，Elsevuier，Amsterdam，159−175．

Allen，R．V．，1978：Automatic

earthquake

recognition

andtimingfromsignaltrace．BulL

Seism．Soc．

Am．，68，1521−1532．
Carpenter，E

W．，1965：Ahistorical

Carpon，工，1979：Investigationofreviewlong

ofseismic

Periodnoise

arraydevelopment．Proc．IEEE．，53，1816−1821．

at

the

large

aperture

seismic

array．J．Geophys．

Res．，74，3182−3194．
Crampin，S．and

C．J．Fyfe，1974：Automatic

analysisoftaperecordingsfromseismic

networks．

Geophys．」．Roy．Astron．Soc．，39，155−168．
Freiberger，W．F．，1963：An
Goforth，T．and

approximate

E．Herrin，1981：An

methodautomatic

insignalseismic

detection。Q．App1．Math．，20，373−378．

signaldetection

algorithm

based

on

the

Walsh

tramsform．Bull．Seism．Am．，71，1351−1360．
Green，Jr．P．E．，R．A．Frosch

seismic

andC．F。Romey，1965二Principles

ofanexperimentalofcmstallarge

aperturearray，Proc．IEE．，53，1821−1833．

Ishii，H．，1976：Applicationofprediction

methodfor

analysis

movement．J．Geod。Soc．

Japan，22，299−301．
Ishiguro，M．，H．Akaike，M．Ooe

Proc．9th
Kitagawa，G．andandlnt．Sympo．on

S．Nakai，1981：A

Earth

H．Akaike，1978：A

Tide，New

Procedure

for

fittingofBayesianapproachtotheanalysis

ofearthtldes．

York．
the

modelingofnon−stationary

timeseries．Ann。

Inst．Statist．Math．，30−B．
Ozaki，T．and

H．Tong，1975：0n

analysis．Proc．8−th
Stewart，S．W．，1977：Real

Bu1L

the

Hawaii111t．Conf。ontime

detectionandnon−stationary

SystemlocationofautoregressivemodelsintimeseriesScience，224−226．
local

seismic

eventsincentralCalifomia．

Sesm．Soc．Am．，67，433−452．

Stewart，S．W．，W．H．K．Lee

quakes

along

the

San

and

J。P．Eaton，1971：Localion

Andreasfault

systeminandcenralreal−time

detectionofmicroearth−Califomia．Roy。Soc．New

Zealand，9，205−209．

3．2信号検出問題
信号検出の問題は、通信やレーダーなどに関連して古くから研究されており、そのモデルは種々考案されている。ここでは続計的手法による検出モデルを考えるが、これは一種のクラス分類（識別問題）で、比較・判断の基準をいかに合理的に設定するかという問題［Wald（1950）、Rao（1977）］

一59一

気象研究所技術報告

第16号

1985

に帰する。

信号検出問題は，信号か否かを判別する2者強制選択問題であり、その判別・決定結果は、
i）

正報

（Hit）

二信号を信号と正しく識別ii）誤報（Miss）二信号をノイズとして見落とすiii）虚報（False）ニノイズを信号として誤って判断

iv）実報（Quiet）ニノイズをノイズと正しく識別の4種に分類できる。これらを真のクラス（母集団）と、決定したクラスとの対応で示すと次のようになる。

．．夏のクラスー一＼、、
決定したクラス＼

S（信号）N（ノイズ）

S（Signal）

正報（Hit）

虚報（False）

N（Noise）

誤報（Miss）

実報（Quiet）

統計的検定問題として、これらHit，False，Miss，Quietの確率を考えると、帰無仮説として信号

S、対立仮説としてノイズNを仮定した場合、Hitは帰無仮説を採用する信頼区間に、Missは帰無仮説を棄却する第一種の誤りに、Falseは第二種の誤りに、そしてQuietは検出力に対応している。
帰無仮説と対立仮説を入れ替えた場合は、これらの対応は逆になる。

ここでいう信号SとノイズNは、情報源空間として2つの排反する事象を含むものを考え、一方をクラスS、他方をクラスNと言葉上分類したものであって、クラスNは通常の雑音源だけを想定しているわけではない。これら2つの情報源のうち、特に知りたいものを信号Sとし、他方をノイズNとしており、クラスSは刺激とも呼ばれる。このため、信号検出問題においては、識別結

果をHit，False，Miss，Quietと分類している。

3．2．1ベイズ（Bayes）の識別手法
ある情報をもとに、排反な有限個のクラス（群、class）のうち、・それが属すると思われる1つのクラスを決定する問題は、識別問題（discrimination

problem）と呼ばれている。

情報（観測量あるいは特徴パラメータ）ベクトルをx、識別されるべきクラス数をk個とし、各クラスをA、，A2，…，Ak（Ai∩Aj二φ，i≠j）で表す。また、Aiの出現確率をP（Ai）、Aiからxが出現する確率をP（xlAi），xであるときAiと判断される確率をP（AiIx）とする。P（Ai）、P（Aiはそれぞれ事前確率（priorprobability）、事後確率（posteriorprobability）と呼ばれ、P（x

l

l

x）

Ai）

は尤度（1ikelihood）と呼ばれることもある。

識別問題において、一般的に損失という概念が導入される。これは、誤って別のクラスのもので一60一

気象研究所技術報告

第16号

1985

あると判断したときに被る損害である。いま、A、からの情報をA、のものであると識別したときの損失をCijで表すことにする。Ciiは、正しく識別された場合の損失である。識別することにより利益（得点）が得られる場合は、負の損失と考えることもできる。損失は、経済的観念から、費用と呼ばれることもある。損失Cijを要素とする行列C＝［C、j］は、損失行列と呼ばれている。

一般に正しく判断された場合の損失は最も小さく

Cii＜Ci」orC」i，（i≠j）（3．2。1）である。
弩いま、P二（P（A、lx），P（A、lx），…，P（Ak

1＝（4，あ，…

と定義する。ここで記号

l

x））

，lk）

とおき、損失ベクトル1を、＝CP（3．2．2）は転置を表す。損失ベクトルの成分るは、Aiと識別したときの平均的損失を表わしている。

合理的な一つの識別方法は、判別することによって生じる損失が最小となるように行うことである。識別のための決定ベクトルをd＝（d、，d2，…

，dk）とする。diは、情報xがAlに属していることを意味し、

0≦di≦1、Σdi二1である。このベクトルdを用いると、全体の平均的損失関数Lは、（3．2．3）L＝d71と書ける。したがって、情報xがA、，A2，・一，Akのいずれのクラスに属しているかを識別する問題は、Lが最小となる決定ベクトルdを決める問題となる。
、

いま、決定ベクトルdをd＝（廟、，虜2，…

とすると、吻n，廟k），虜j二〇（i≠」），虜j＝1（i＝j）

L＝桑となり、損失ベクトルの要素のうち最小となるものを求めることと、情報xの判別が同等になる。これらのことを整理すると、

4，あ，…

Jkのうち最小のものが4→x∈Ai（3．2．4）となる。ただし、x∈AiはクラスAiと識別したことを意味する。
（3．2．4）の識別手法は、事後確率P（AiIx）により構成されている。すなわち、事後確率が既知

り場合は、合理的に識別基準を設定することができることを意味する。事後確率を知ることは、一般には容易ではないが、事前確率P（Ai）が既知の場合はベイズの定理により推測することができる。
事後確率P（Ailx）をP（Ai）P（xlAi）で置き換えた識別基準は、ベイズ解と呼ばれている。

一61一

気象研究所技術報告

第16号

1985

識別基準を幾何学的立場から見ると、n次元のxの空間をk個の部分空間に分割することと同等になる。部分空間の境界面は識別境界面と呼ばれている。一般的には、事前確率は末知な場合が多い。

そのため、この識別境界面をいかに設定するか（決定方式）という問題が生ずる。識別境界面さえ合理的に設定できれぽ、識別は機械的に行える。しかし決定方式のうち、優劣が単純比較不可能なもの（許容的決定方式）が存在する。これら許容的決定方式の中から先験確率が既知という条件のもとで選択された最適決定方式がベイズ解である。
（3．2．4）の識別手法を信号検出の二者択一問題について具体的に適用する。母集団（情報源）を

排反事象の信号SとノイズNとし、損失行列Cを

CC
CC
S
S

N N
S N

〔 〕
一一

C

NS

とする。Css，CNNは正しく判断された場合の損失、CsNとCNsはそれぞれFalseとMissの場合の損失である。損失ベクトル1＝（IS，JN）の成分は、

Js＝Cssp（SIx）十CsNp（NI

x）

lN二CNsp（SIx）、十CNNp（NI

x）（3．2．5）となる。したがって、識別は、

Js＜JN→x∈S（3．2．6）ls〉IN→x∈Nで行われる。ここで、Js−JN二｛（CsN−CNN）P（N）P（xlN）一（CNs−Css）P（S）P（xlS）｝／P（x），

よって（3．2．1）式より、（3．2．6）式の識別手法は次のように書ける。

A＞β→x∈S（3．2．7）A＜β→x∈N
ただし、

A＝P（x

I

S）／P（xI

N），

β＝（CsN−CNN）P（N）／（CNs−Css）P（S）。

Aは尤度比（likelihoodratio）とよばれ、（3．2．7）は統計学では尤度比検定として知られている。

また、Aの対数は情報量の測度となる。

一62

気象研究所技術報告

第16号

1985

識別基準の設定としての立場から、ベイズ解を見てみよう。簡単のために、信号検出の二者択一問題を考えることとし、Sの識別領域をRs、Nの識別領域をRNとして考える。識別結果の期待値；

Hit：P（H）一P（S）みSp（xlS）dx

Miss：P（M）一P（S）みNp（xlS）dx
False：P（F）一P（N）録sp（xlN）dx

Quiet：P（Q）一P（N）遜Np（xlN）dx
なので、損失関数L。は
Le＝CssP（H）＋CNsP（M）＋CsNP（F）＋CNNP（Q）（3．2．8）となる。したがって、L，を最小とするような領域Rs，RNを決める問題となる。いま、識別境界面をx＝θとすると、Leを最小とする条件は、（CNs−Css）P（S）P（θI

S）＝（CsN−CNN）P（N）P（θl（3．2．9）N）となる。これがベイズ解である。したがって、ベイズ解を用いた識別手法は、

A＞β→x∈S（3．2．10）A＜β→x∈N
ただし
A＝P（x

l

S）／P（xiN），β＝（CsN−CNN）P（N）／｛（CNs−Css）P（S）／、となり、（3．2．7）に相当している。xがスP（xlS〉

P（xiN）

カラーの場合の概念図を図3．2．1に示す。

ベイズ解に基づく識別手法は、統計学における最良尤度比検定に相当する。これは、二者択一問題において、事前確率が既知でXe（㌔←一一一一RN一一一→1←一Rs

あれば、Neyman−Peasonの定理に基づく
図3．2．1 情報xが1次元連続で、μS＞μNの最良棄却域を設定できることを示している。

場合の識別境界面の概念図。

次に事前確率が未知の場合の識別基準のP（x設定を考えてみる。このときの制約条件は、P（S）十P（N）＝1

lN）、P（xlS）は情報源NおよびSの密度関数。RN、Rs、はそれぞれのクラスの識別領域を示す。島は、（CNs−Css）＝（CsN−CNN）で、P（N）；（3．2．11）

P（S）の場合のベイズ解を示す。

63

気象研究所技術報告

第16号

1985

である。この制約条件のもとに（3．2．8）のL，の最小化は、婦、p（xlS）dx＋蟻P（xIS）dx（3．2．12）

一C頭sp（xIN）dx＋CNψ，p（xIN）dxのとき実現される。識別境界面をx二θとすると、この場合の識別基準βは、
β二P（θlS）／P（θlN）（3．2．13）となる。ただし、θは（3．2．12）式を満たす。これは、窺初一窺砿解と呼ばれている。

識別手法として、ベイズ解を用いた手法はベイズの識別手法、窺1n一窺砿解を用いた手法はアn吻

琢一吻砿識別手法と呼ばれている。しかし、窺初一㎜解はあるベイズ解と一致するので、本章では、これら両者合わせてベイズの識別手法と呼ぶことにする。

3．2．2

識別不能状態信号検出問題を一般化し、ここでは信号検出問題の決定出力に、識別不能状態Uを設けた場合について扱う。すなわち、情報源はクラスSとクラスNの2つとし、最終判断として、S、N、Uの3つの状態のいずれかを選択する問題について述べる。
簡単のために、正しく判断する場合の損失をCss＝CNN二〇、誤って判断する場合の損失をCNS二CsN〒Ce＞0、判断不能とする場合の損失をCus二CuN＝Cuとする。また、決定ベクトルをd＝（ds，dN，du）フとすると、損失関数Lは（3．2．2）、（3．2．3）式より、

L＝Ce（1−L）

P（3．2．14）

ここに〜

L＝dsP（Slx）十dN（NIx）十du（C，十Cu）／Ce（3．2．15）となる。Lの最小化は、Lの最大化と同値である。du＝0の下での最大化、すなわちS、Nのいずれと識別するかについてはすでに述べたので、ここでは識別不能状態Uと決定される場合を問題とする。

識別は損失C，とCuの大小により、次のように分類される。

（1）Cu＞C。
常に、
P（S

l

x）〉（C，一Cu）／C，〈P（N

したがって、識別不能の場合は生じない。
（II）C。＞Cu＞0

一64

I

x）

（

気象研究所技術報告
（i）

P（S

I

x）＜（C．一Cu）／C，＞P（N

I

第16号

1985

x）

識別不能。
（ii）

（i）

以外の場合

識別不能ではない。実際P（SIx）＝P（NIx）二〇の場合は生じないであろうから、C，二Cuの場合も識別不能ではない。
（III）Cu≦00≦P（SIx），P（NI

x）≦1

したがって、常に識別不能。

これらのことから、識別不能状態の損失が、誤りを犯した場合の損失以上であれば、識別不能状
態は生じない。また、識別不能の状態を損失としなければ、常に識別不能とすればよく、何もしないことと同じである。

実際問題として、信号検出問題はクラスSか否かを知ることが目的であり、最終的識別クラスとしての識別不能状態（U）は、後の変化を生せず好ましくない。したがってクラスUと判断された
場合は、人問が介在するにしろしないにしろ、情報の抽出法あるいは情報源の状態を変更する等、次に何らかの処理を経由して、クラスSか否かの判断を下すことになる。
クラスUと判断された場合の次の処理の一つに、判断の保留という処理が考えられる。これは、次の新たな情報が得られるまで判断を保留し、その情報も含めて再び識別を行う手法である。すな

わち、識別結果の確度が低い場合は、次の情報も含めて標本数を多くすることにより識別の確度を高める手法である。このような識別手法は、逐次識別手法と呼ばれている。この手法を適用するためには、同じ情報源から連続して情報が発生されることが大前提となる。

情報源に識別不能状能のクラスを設定した場合は、3つのクラスからの情報を3つのクラスのうちいずれかに決定することであり、もはや識別不能の問題ではない。

参考
Rao，C．R．，1977：Linear

statistical

Wald，A．，1950：Statistical

3．3

inference

decision

and

文
its

献
applicalions（2ndfunctions．Wiley，Newed．）．Wiley，New

York．

York．

雑音の中の信号検出

本節では、得られた情報が単に情報源Nからの出力が加わった出力

N＋S

N

か、これに情報源Sからの出力

S

なのかを識別する問題を考える。

時刻tにおける情報源Nの出力をnt、情報源Sの出力をSt、我々の測定系への入力を簸とする。
時刻パラメータtは整数の全体Zを動くとする。識別に用いる時刻を／1，2，．．．解／とし、識一65

気象研究所技術報告

第16号

1985

別区間と呼び、記号丁で表す。

情報源N、Sの母集団が共に知られている時、識別クラスをクラスN二Xt…nt，t∈Z（3．3．1）

クラスS：Xt≡nt十St，t∈Zとして、検定すべき仮説を帰無仮説

H。：Xt≡nt，t∈Z（クラスN）

対立仮説

H1：xt≡nt＋st，t∈Z（クラスS）（3．3．2）と設定する。

H。（クラスN）かH、（クラスS）かの識別はベイズの決定手法に従い次のようになる。

A＞β→H、（クラスS）（3．3．3）

A＜β→Ho（クラスN）
ここで

A一＿P（Xt，t∈TIH1）
P（Xt，t∈T『H。）（3．3．4）であり、βは識別基準を与える定数である。

もし、情報源Sの母集団が不明ならば、対立仮説H、（クラスS）は、

H1二Xt羊nt，t∈Z（3．3∫5）3．3．1正規定常情報源に対する信号検出手順
情報源N、Sが互いに独立な正規定常過程に従う場合に識別手順を具体化しよう。
最も簡単な場合を考える。すなわち

条件A二叛（渉∈Z）は互いに独立で平均0分散σ祷の正規分布ψ（0、品）、St（孟∈Z）も互いに独立で平均μs分散σ§の正規分布ψ（μ，、σ書）にそれぞれ従う。

まず各時点毎に識別を行うことを考える。この場合識別区間丁は1時点｛1／よりなる。仮設H。
（クラスN）、H、（クラスS）を（3．3．2）のように設定すると、クラスNに対する尤度は、P（xllH∂一嘉exp｛一2導貧｝で与えられる。クラスSに対する尤度も簡単に求めることができ、

一66一（3．3．6）

気象研究所技術報告

第16号

1985

P（為陶』。．d鉦烹呵一（Xl嘉）2／毒￠功／一（s蓑）2／一2π（議＋σ粛）吻／一捲漏／『となる。対数尤度比ln（3・3・7）Aは

lnA−1％詣2（σ繹σ蕎）＋σ尋σ畜｛＋2讐斎x多／（3・3・8）ゆとなり、識別手方は、
愚

L＞α→H1（クラスS：信号有り）
L＜α→H。（クラスN

l信号無し）（3．3．9）で与えられる。ここでσ甕 2
L＝x1μs十2x1
2σ N（3．3．10）α一（σ§＋σ貧）／lnβ一ln詣／＋ヂである。この場合識別は真の信号の平均値と観測値との相関および観測値のエネルギーにより行われる。したがって、平均信号の絶対値が大きく分散が小さければ識別能力は高い。そうでない場合、各時点毎の観測値を用いる識別手法の確度は低い。

今、信号をm個連続して受けるとしよう。この場合仮説（3．3．2）に対する尤度はP（一一痴隅）一邑烹呵一2弩昌／（3．3．11）

P（一・商IHD一且蒲ゆ｛一2潔詣／となる。したがって対数尤度比ln

Aは1％A−m／l％識一2（σ搾σ為）／＋（σ§｝σ畜）／藤＋2繋xそ｝（3・3・12）となり、（3．3．3）式にしたがって識別が行われる。これは、前節の逐次識別手法に相当する。
二つの判定方式の優劣を具体例で見ておこう。μ，＝0、P（N）二P（S）、Css＝CNN＝0，CsN＝CNsかつ
σ蕎＝σζとする。各時点毎の識別方式では、P（Hit）≒0．49，P（Miss）≒0．51，P（False）≒0．24，P（Quiet）≒0．76となる。いいかえれば信号無しを信号有りと誤るのが4回に1回、信号有りを信号無しと誤るのが2回に1回でその識別結果は殆ど信用できない。m時点判別方式では同じ前提条件下

一67一

気象研究所技術報告

第16号

1985

でm二100とすると、P（Hit）≒0．99，P（Miss）≒0．01，P（False）≒0．01，P（Quiet）≒0．99となる。

mは検出したい信号とノィズのS／N比と、期待する検出の信頼度により決まる。

次に条件Aを要求しない、すなわち各時点毎の独立性を仮定しない場合にm時点識別方式を以下に記しておこう。正規定常過程という仮定からn二（n、，n2，…，nm）
ΣN＝E［m］のm次元正規分布、s二（s、，s2，…，sm）

列Σs二E［（s一μ）（s一μ）

は平均ベクトル0，共分散行列

は平均ベクトルμ＝（μs，μs，…，μs）、共分散行］のm次元正規分布に従う。ここでE［・］は平均を表わす。nとsの独立性を用いるとx二n＋sは平均ベクトルμ共分散行列ΣN＋Σsのm次元正規分布に従う。仮説（3．3．2）の下で、クラスN及びクラスSの確率はP（一・ヲ服）一（2π）叢≠ふ）券呵一麦どΣ虻x｝（3．3．13）P（一・㌔IHD一（2π）糀！複＋ふ）｝参吻｛一壱（x一μ）7（糸＋ふ）一1（x一μ）｝
で表わされる。したがって対数尤度比は1％A一麦1㌦渉窓転）＋去｛緩1x一（x一μ）・（糸＋ふ）一・（x一μ）／（3・3・・4）となり、識別方式は

L＞α→H、（クラスS）

L＜α→H。（クラスN）（3．3．15）となる。ここにL＝xΣユx一（x一μ）

N（ΣS十Σw）一1（x一μ）

伽ΣN

α＝21箆β一1％46渉（Σ、＋ΣN）（3・3・16）（3．3．15）の識別手法がスペクトルを反映していることは、Wiener−Khintchineの公式より明らかであろう。ΣN二σ鳥1，Σs二σ§1（1：単位行列）ならば条件Aを満たし既に述べた場合に帰着する。

3．3．2情報源のモデル
前節において情報源の確率構造は既知とした。しかし実際上は観測値からそれを推定しなければならない。一つの手がかりは観測時系列のスペクトルであるが、周波数領域からのアプ・一チは非定常あるいは非線形現象に対して心ずしも有効ではない。我々は時間領域からアプ・一チする。

時間領域での情報源モデルとして取り扱いが簡単なものにマルコフ的情報源がある。情報源Xの時刻tの出力をx此する。

IXt／が多重マルコフ過程、すなわちある自然数Lが存在して時刻t＋1

における出力Xt＋・の確率分布がXt，Xt一・，…Xt−L＋・のみに依存してXt−L，Xt−H，…にはよらない時、Xをマルコフ的情報源とよぶ。いま、t＝（xt，xt一、，…，xt−L＋、）

一68と考えると、xt＋1の分布はxtのみに依

気象研究所技術報告

第16号

1985

存し、単純マルコフ過程になることを注意しておこう。
さて条件付き確率P（xtI

xt一、，xt−2，…，xt−L）はxt，xt一、，…，xt−Lのみからきまるのだから適当なL変数関数fと正規確率変数εtを用いて
xt二f（xt−1，xt−2，…，xt−L）＋εt（3．3．17）と表されるであろう。fが線形関数のとき線形時系列モデル、非線形関数のとき非線形時系列モデルという。線形時系列モデルで表される情報源は正規定常過程に従い、（3．3．1）項で検討した範疇にはいる。
乃

線形時系列モデルとしてよく知られているものに、自己回復（AR）モデル、移動平均（MA）モデル、これらの混合型である自己回帰移動平均（ARMA）モデル、積分型自己回帰（IAR）モデル、

積分型移動平均（IMA）モデルおよびこれらの混合型である自己回帰積分型移動平均（ARIMA）モ
デル等がある。

AR，MAおよびARMAモデルは共分散が発散しない、つまり平均値のまわりで振動する時系列
データを表現するのに適している。もちろんこの3者の中ではARMAモデルが最も優れているが、MAモデルも含めて推定は困難である。これらに対してARモデルは最も取り扱い易いモデルである。IAR，IMA，ARIMAモデルはトレンド成分をもつような時系列データを表現するのに適しており、ARIMAモデルが最も優れている。IAR，IMA，ARIMAモデルのうち、最も取り扱い易いのは
IARモデルである。これらのうち、どの時系列モデルを採用するかは観測データにより決まる。例

えぽ地震波信号に対してはARモデルが、地殼変動データにはIARモデルが適している。
線型マルコフ的情報源は、次の状態空間表現により一般化して表される［例えぽ、Akaike（1974）、
Kalman（1968）］。

Yt＋1二FtYt十GtVt（3．3．18）

Xt＝HtYt十Wt
ただしYt

二状態ベクトル

Xt

二観測値ベクトル

Vt

二遷移確率ベクトル

Wt：観測誤差ベクトル

Ft，Gt，Ht二状態ダイナミックス
状態空間表現は、ダイナミックシステムとして見ることができる。Ft，Gt，Htが時間によらず一定であるものは時不変システム、時間と伴に変化するものは時変システムと呼ぽれている。（3．3．19）の表現により、上述の時系列モデルを含むより一般的な時系列データを表現することができる。本章

一69一

気象研究所技術報告

第16号

1985

では、通常用いられる定常・非定常の概念とは異なり、状態ダイナミックスが時不変として表現され、かつ遷移確率ベクトルの従う分布が不変である場合を定常、それらのいずれか一方でも時変の場合を非定常と呼んでいる。

3。3．3

エントロピー最大化原理と赤池の情報量基準AIC

前項で想定したマルコフ的情報源の内どれを選ぶべきかという問題を考える。これは、想定したモデルのうちどのモデルを選択するかという問題である。これに対する有力な基準の一つに赤池による情報量基準AIC［Akaike（1973，1977）］がある。これはエント・ピーを最大にすることによって最良のモデルを選択しようとするための基準量である。以下、AICについて解説する。

真の分布と近似分布の差を測るため、kullback情報量

1（f：9）一乃z伽姜1劣（3・3・・9）を導入する［kullback（1959）］。ここでfは真の密度関数、gはあてはめたモデルの密度関数である。
この情報量は性質；
1（f：9）＞O（3．3．20）I（f：9）二〇⑬f＝9を持つため、近似モデルgは1（f；g）→minとなるように選べばよいことがわかる。一1を一般化されたエソト・ピーと考えれぽこの手続きはJaynes（1957）によって提唱されたエント・ピー最大化原理にほかならない。

今近似モデルgはパラメータθのみに依存するとしよう。すなわちg二g（Z，θ）とする。
一1（f：9）一一∫旨zf（z）1盟f（z）一∫ヒxf（z）1％9（z）（3。3．2・）において右辺第1項は真の母集団のみに関する項でgにはよらないので、第2項のみを最大化すれぽよい。

真の分布密度fを持つ確率変数をN回独立に観測して値Xi（i二1，。．N）を得たとする。大数の法測より

撫N一・ゑ1宛9（xbθ）一∫dzf（z）1％9（ろθ）（3・3・22）となる。ところで平均対数尤度を1（θ）とすると
1（θ）＝N｝1Σ1％g（Xi，θ）

i二1

70

（3．3．23）

気象研究所技術報告

第16号

1985

したがって∫dzf（z）1箆9（励を最大にすることは＋分大きいN耐して畑撮大化することを意味し、最尤法は近似的にエント・ピー最大化原理と同等である。

さて・撃∫dzf（Z）1％9（励を実現するθを島とし真のパラメータと呼ぶことをこする・こ砥g（z，島）＝f（z）とは限らないことを注意する。また対数尤度1（θ）を最大にするθをθとする。θは最尤定量とよばれる。θは島に十分近いとして∫dzf（z）lng（か∫dzf（z）」κ9（鴫）＋声zf（z）〔蓄1（賜）〕（秒一島）＋吉（み一島）∫dzf（z）〔∂語び1（3・3・24）（賜）〕（3一島）と展開する。ここでθ。の定義より
z声f（z）一1多z∂g（z，θ。）＝0∂θ

（3．3．25）

したがって（3．3．25）は∫dzf（z）1π9（吻一∫dzf（z）1（脇）一去（碗」（命）（島一島）（3・326）となる。ただし」（島）一一∫dzf（z）∂藷〆η9（筋）（3・3・27）

ここで声zf（z）1η9（z，島）は∫dzf（z）1％9（筋）一西zf（z）1％9（賜）一N一・菖1％9（x謝十N−1Σhz
i＝1g（x茎，θ。）一N−1Σln g（xi，のi＝1

ハ（3．3，28）ハ十N−1Σ1銘9（Xi，のi＝1と書ける。N−1呈1％g（x、，島）を3のまわりで展開すると

iニ1
N一・菖1銘9（xb島）雲N一・急♂％9（xi・初＋N−1鵠1％9（畑（曲（3．3．29）＋⊥（仇一初・N−1至∂21銘9（x、，あ（島一初2i一・∂θ∂θ
最尤推定量の定義から

N一・呈一塗一1％9（Xi，あ一〇
i−1∂θ

一71一

（3．3．30）

気象研究所技術報告

第16号

1985

またNが十分大きいとし、θ≧島より
N−1Σ N∂21η9（xi，の2−J（θ。）＾（3．3．31）i−1∂θ∂θ

一方中心極限定理より
N一麦呈』L lπg（x、，島）雲平均0共分散1（島）の正規確率変数i−1∂θ

（3．3．32）

ここに1（θ。）はFisherの情報行列で1（島）一∫dzf（z）券1箆9（謡）島1％9（協）（3・3・33ン

で定義される。したがって、（3．3．29），（3。3．30），（3．3．31），（3．3．32）より〉下「（み一θ。）は平均0

共分散J−1（θ。）1（θ。）の正規分布に従うことが分かる。故に、E［（θ一θ。）

J（θ。）（θ一θo）］ンN−11zα66（J『1（θ。）1（θ。））（3．3．34）

ここにE［・］は平均を表わす。またN一・菖lng（X離∫ヒzf（z）1η9（幅）よって（3．3．24）式はNが十分大きい時に声zf（z）1κ9（吻＝諭一N一∫勉66（J−1（島）1（島））（3・3・35）となる。
もし9（z，θ。）＝f（z）ならぽIJ（θ。）＝1（θ。）となるので

加06（J（島）一11（θ。））＝K

ここでKはパラメータの次元である。これを（3．3．25）式に代入すれば
z

＾
一＾
声
f（z）1錫9（z，θ）21（の一一KN（3．3．36）を得る。これを一2N倍したものが赤池の情報量基準AICである。
AIC二一21（の十2K
AICが小さい程dz

f（z）ln（3．3．37）g（z，θ）は大きくなり9は真の分布に近いとことがいえる。

赤池の情報量基準AICによりパラメータの次元の異なるモデル間の比較も可能になりパラメトリッ
クなモデル構成がパラメータの個数を含めて行える。AIC導出の過程からわかるようにデータ数が少ない場合や想定したモデルが真の分布と掛け離れている場合はAICをそのまま適用することはで72

・

気象研究所技術報告

第16号

1985

きない。

最後に統計的検定問題をエント・ピー最大化原理の立場から見ておこう。仮説H。、H、に対する確率分布をP（・lH。）、P（・lH1）、真の分布をP（・）とする。Nが十分大きい時

一1（P（・）：P（・iHo））2N『1Σ1解P（xilHo）一N−1Σ1％P（xi）
i＝I

一1（P（・）二P（・li＝1H1））＝とN−1Σ1％P（xii＝1lH1）一N一1Σ1％P（xi）
i＝1

したがって
N｛一1（P（・）二P（・lnP（x五iH、）H1））十1（P（・）：P（・lH。））／＝伽i看1

11P（Xi旧。）
i＝1
これは尤度比検定と同等であり、先に述べたベイズの識別手法とも等しい。

3．3．4

マルコフ的情報源モデルによる信号検出xt（t∈Z）において、x、，x2，…，xnはクラスNであるという条件のもとに、x。＋、，xn＋2，…，xn＋mがクラスSか否かを識別する問題を考える。時刻パラメータtに関して、T＝｛1，2，・一，n＋m／、T・＝｛1，2，…，n／、T2＝／n＋1，n＋2，…，n＋m｝とおくと、識別クラスは、

クラスN：Xt二nt，t∈T（3．3．38）

クラスS：Xt＝nt十αtSt，t∈T
ただし儂謝となる。したがって検討すべき仮説は、、帰無仮説H。二Xt・＝nt，t∈丁対立仮説H、：xtニnt＋αtst，t∈T｛1：1詔
このとき、それぞれの丁上での結合確率は、
P（xt，t∈TIH。）＝P（x1，x2，…，xn，…，xn＋mlN）（3．3．39）P（xt，t∈TIH1）＝P（x1，x2，…，xnlN）P（xn＋1，xn＋2，…，xn＋mlS①N）となる。ここにP（・lN）、P（・lS①N）は、それぞれ情報源Nおよび合成された情報源S①N

一73一

気象研究所技術報告

第16号

1985

からの出力と考えたときの結合確率である。

情報源モデルを通して情報を測る立場から上の識別問題を考えてみる。まず、情報源モデルの推定も含め同時に上の識別問題を解決する手法について述べる。記法を簡単化するため、H・に対する

情報源モデルの状態ダィナミクスおよび遷移確率の共分散列をそれぞれFT，GT，CTで表わす。また

H・に対しては、FT1，GT1，CT1、およびFT2・GT2・CT2とする・各情報源モデルのパラメータ次元は異なるため、AICを通じてこれらパラメータ推定の実現を図る。実現されたパラメータに記号をつけて表わすと、H。，H、に対するAIC（H。）、AIC（H、）は、
AIC（H。）＝（n十m）1η4αCT十2κT十const．
（3．3．40）
もAIC（H1）＝nln46渉CT、十m吻4切CT、十2（κT、十κT2）十const・となる。ここに、κT，κT、，κT、は、推定情報源モデルのパラメータ数である6AIC（H・）・AIC（H・）の大小によりクラスN、クラスSの識別が行える［OzakiandTong（1975），Kitagawa

andAkaike（1978）］。

3．3。3項で述べたよ鬼呵一がIC掩尤度の不偏推定量を目指したものである・そ砥尤度をAICにより表わすと、ベィズの識別手法は

A＞β→H1（クラスS）（3．3．41）

A＜β→Ho（クラスN）
ただし、

A−exp｛一青AIC（H・）＋去AIC（Hの／となる。

次に、情報源N（マルコフ用情報源）は既知という仮定のもとに｛Xt，t∈T／T＝／1，2，…，m／がクラスSか否かを識別する方式について考えてみる。この場合の識別クラスは、（3．3．38）においてTをT2に限定したものとなる。情報源Nによる予測誤差をWtとすると、新たな識別クラスは、
クラスN二｛Wt，t∈T2／⊂x：ψ（0，CN）（3．3．42）

クラスS：IWt，t∈T2／d〔ψ（0，CN）

ここにψ（u，C）は平均u共分散行列Cの白色正規密度関数を表わし、｛・／俣ψは／・／が、ψに従うことを意味する。クラスSに対してWtは新たなマルコフ的情報源｛Fw，Gw，Cw／で表現されるとすると、クラスN、クラスSに対するAIC（H。），AIC（H1）は、AIC（Ho）＝ml％4αCN十ΣwlC西lwt十const．
＾
t＝

AIC（H1）＝ml箆4αCw十mr十2κw十const．

一74一（3．3．43）

気象研究所技術報告

第16号

1985

となる。したがって、この場合の識別手法は、AICを用いて尤度を表わすことにより（3．3．42）と同様になる。

さて、（3．3．43）の識別をクラスSの情報源モデルの推定により行ったが、E［w，wl＋．］＝0、（τキ0）と仮定し、単に共分散行列の検定のみを行っても十分実用に耐える。この場合wl

C西l

w，が

π2分布に従うことを利用して検定を行なってもよい。また、共分散行列Cが未知であるとして、F分布を用いることもできる。

ここで述べた手法は、次節の信号出現時刻の決定に応用されている。

参考文献
横田

崇二周勝

奎・溝上

恵・中村

功、1981：地震波データの自動検測方式とオソライソ処理システム

における稼動実験、地震研究所彙報、55、449−484。
Ak＆ike，H。，1973：Information

Intemational

Symposium

theory

on

and

an

extension

lnformation

of

the

maximum−1ikelihood

Theory，eds．B．N．Petrou

andp口nciple．2nd

F．Csaki，Akademia

Kiado，Budapest，267−281．
Akaike，H．，1974：Markovianreprensentationofstochasticprocessesandofautoregressive

Akaike，H．，1977：0N

movingaverageentropyits

applicationtothe

analysisprosses．Am．Inst．Statist．Math．，26，363−387．

maximization

principle．Proc．SymposiumonAmsterdam，Norm−Holland，27−41．
Jaynes，E．T．，195711nformation

Kalman，R．E．，1968：New

theoryandstatisticaldevelopmentinsystemsmechanics．Phys．Rev．，106，620−630。

theoryBiology，ed．M．D．Mesarovic，Springer−Verlag，New
Kitagawa，G．and

H。Akaike，1978二A

procedure

for

the

relevanttobiology．Systems

TheoryandYork，222−232。

modelingofnon−stationaly

timeseries．Am。

Inst．Statist．Math．，30−B
Kullback，S．，1959：Infomation

theoryOzaki，T．andthe

H．Tong，1975：0n

analysis．Proc．8−th

3．4

Hawai

and

fitting

statistics．Wiley，Newof

lnt，Conf．on

non−stationary

Sytem

York，

autoregressivemodelsintimeseriesScience，224−226。

信号出現時刻の決定

信号出現時刻は、地震波信号にあっては、初動時刻あるいは別種の波の混入という意味で位相混

入時刻［横田ら＜1981）1とよばれている。それを客観的に精度よく知ることは震源を定める上で不可決であり、決定手法について数多くの研究がなされている。

非統計的手法としては、渡辺ら（1975）、Stewart（1977）、Allen（1978）、Anderson（1978）等がある。そこではある種のフィルター出力による特徴関数を内外挿して初動時刻が求めてられてい一75

気象研究所技術報告

第16号

1985

る。統計的手法による試みとして、白井・徳弘（1979）によるベイズの手法に基づくものがある。

しかしこれは、1時点毎に信号か否かの判別を行う手法のため、すでに前節で述べたように、その確度は十分とはいえない。横田ら（1981）は、初動時刻＝入力波の定常性が破れる時刻であることに注目し、非定常時系列の一般論を応用した決定手法を開発した。

信号出現時刻の決定は、信号検出問題においてクラスSと識別された後に実行される。基本的には前節までの議論の応用により解決される。白井・徳弘（1979）の手法を1時点識別手法、横田ら（1981）の手法を多時点識別法として区別する。

3．4．1信号出現時刻決定の統計モデル

亀横田ら（1981）に従って述べる。横田ら（1981）では、特に地震波信号に主眼をおいているため、情報源として最も取り扱いやすいARモデルを採用しているが、ここではマルコフ的情報源という一般的制約のみをおく。

信号出現時刻の決定問題は、定常・非定常状態を判別するという立場から次のように考えることができる。図3．4．1に示すように任意に与えられた区間ACにおいてB時点から信号が混入した場

合を考えよう。すなわち定常な情報源からの出力にB地点からやはり定常な情報源Sからの出力が混入したとする。区間AC全体は非定常であるが、区間AB及びBCのみでは、定常な時系列である。言いかえると、信号出現時刻の決定は任意に与えられた非定常な区間を定常な2つの区間にわけることに他ならない。

3．3節と同様に、時刻tにおける情報源

A

B

l

ド．stati。nary

・

l

＿→く。＿

c
stati。nary

Nの出力をnt・情報源Sの出力をSt我々のl一一．．→i
じ

1

測定系への入力をXtとする。時刻パラメ

一タtは整数の全体Zを動く。判別に用いる時間を／1，2，……，m／としてこれをTで表わす。k時点から信号が出現したとす
しこ

1←一一一一一…

！

nonstati・narγ

1

ると

一…一…一一ゆ1：（

xt二nt十αtSt、t∈Z（3．4．1）

図3．4．1 定常な情報源Nからの信号に定常な情︾﹂報源Sからの信号が混入した場合の概念図。情報源Sからの信号の混入時点
をBとすると、区間ABは定常な情報
源Nからの信号で、区間BCは合成された定常な情報源S①Nからの信号となる。区間AC全体は、1つの定常な情報源では表現できず、非定常となこる。

一切持っていない。したがってk時点付近

一76

で

一｛慧一1（3・4・2）

さてSk．、Sk．2……についての情報を我々は

気象研究所技術報告

第16号

1985

でSを通常のマルコフ的情報源とみなすことは妥当性を欠く。むしろ時間の向きを逆転してStの分
布はSt＋、、St＋2、……で決まると考えるべきである。実際、信号出現時刻の決定は信号有りの判断後に行われるもので、情報源Sを表現するに十分長い区間｛k，k＋1，……，m／にわたる情報が既に我々の手にある。このように時間の向きを逆転させて作った情報源を後ろ向き情報源と呼び、記号←をつけて表すことにする。

Xtはt〈kでは情報源N、t＞kでは合成された後ろ向き情報源S㊦NFからの出力である。それぞれの多重マルコフ性の次数をLN、LSNとすると、Xtは条件付き確率の系
P（xtI

xt＿1，xt＿2，…，xt＿L）

t＜kN（3．4．3）P（xtlxt＋・，xt＋2，…，xt＋LNs）t≧kで表わされる。信号出現時刻の決定範囲を丁二｛1，2，…，m／に保つため観測範囲を前後に拡張してTo二｛1−MN，…，1，2，・…・・，m，・一，m＋MsN｝としておく。

｛x｝のTo上での結合確率はP（xt，t∈Tolk）二P（x1＿MN，…，x1，x2，…，xk＿1iN）

×P（xk，Xk＋1，…，xm，…，…xm＋MSMとなる。ここにP（・IくN）、P（・I

←一

I

S㊦N）

（3．4．4）S㊦N）はそれぞれ情報源N、情報源S㊦Nからの出力と考えた時の結合確率である。記法を簡単化するため｛1−MN，…，1，2，…，k−l／＝T9、
m，…，m十MsN／＝T8N、／l−MN十LN，…，1，2，…，k−1／＝丁歯、／k，k＋1，…，／k，k十1，…，m，…，m十MsN−

LsN／＝丁杢Nとし、丁歯∩T§N二丁1とする。また、情報源を3，3，2項の状態空間表現により／F，G，C／

で表わし、情報源N，S㊦Nに対してそれぞれ｛FN，GN，CN／、｛FsN，GsN，CsN／とする。
（3．4．4）に対して、t∈T1上での条件付き確率をP（xt，t∈Tllk，icx）＝P（xt，t∈Tolk）／｛P（x1−MN，…，xLN−MNIN）（3．4．5）×P（xm＋MsN−LsN＋b…，xm＋MsNIS㊦N）／とする。
，

信号出現時刻kの決定は（3．4．4）あるいは（3．4．5）を最大1こするkを求めることに帰着する。その具体的手法を以下に3つ述べる（図3．4．2参照）。

手法1二情報源N、N＋Sを推定するのに判別区間丁の全データを用いる。
まずt∈丁畠に対してP（xt，t∈丁畠IFN，GN，CN）

t∈Ts箔に対して

一77一（3．4．6）

気象研究所技術報告

（1）

第16号

1985

．N

k

t

k−1

S㊥N
（五）

WN

一一一一一→k−1

（皿）

N

皿N
一聡

kぐ

｛

t

︺↑酬

N

WN

一一一一一＞k−1t

k一一一一一一→〉DWN
図3．4．2 信号出現時刻kの決定のための情報源NおよびS㊦Nの推定区間の概念図。

（1）は手法1、（II）は手法2、（III）は手法3に対している。波線矢印はNおよびS㊤Nからの予測誤差Wn（t）、Wsn（t）を計算する方向を示す。（HDの手法3においては、S㊦Nのみの情報源を推測した場合はこの図とは逆向きになる。

P（xt，t∈TgNIFsN，GsN，CsN）

を計算し、最大値を実現するFN、GN、CN、FsN、GsN、CsNを記号（3．4．7）をつけて表わす。このとき、（3．4、4）に対する赤池の情報量基準AICはAIC（k）＝（k−1十MN）1％4αCN十（m−k十1十MsN）1％惚CsN
十2（細十κSN）十const．

（3．4．8）で与えられる。ここで撤、κSNはそれぞれ情報源N、S㊦Nのパラメータの個数である。このAICを最小にするkをもって信号出現時刻とする。

この手法は、判定する時点k毎に2つの情報源を作成しなければならず、計算時間がかかる。

手法2二判定区間の最前部および最後部のデータから情報源N、S㊦Nを推定する。
この手法は次の2つに大別される。
1）

予測誤差の共分散行列をCN、CSNに固定する。

2）

予測誤差の共分散行列AN、ASNを計算される予測誤差より推定する。

78一

気象研究所技術報告

第16号

1985

t時点のNおよびS①Nからの予測誤差をそれぞれwN（t）、wsN（t）とすると（3．4．5）に対するAICは、
1）の場合は、

AIC（k）＝（k十MN）1％惚CN十（MsN−k十1）1％惚CsN（3．4．9）
k−1

＾

MSN

＾

十Σ．w長（t）C『歯wN（t）十Σws脳（t）C一§NwsN（t）十const．

t＝一MN

t＝k

となる。ここでMN＝MN−LN−1
MSN二m十MSN−LSN
とおく。

また2）に対するAICはxの成分数をrとしてAIC（k）＝（k十MN）1％惚AN十（MsN−k十1）ln襯AsN＿＿十（MN十MsN十1）r十2（r十r2）十const．

（3．4．10）となる。これらのAICが最小となる時点kを信号出現時刻とする。
手法1と比較して、計算時間は飛躍的に短縮され、精度もさほど見劣りしない。手法2−1）において、近似情報源の共分散行列によっては結果が多少不安定になることがある。

手法3：判別区間の最前部あるいは最後部のみの情報源を推定する。

簡単のため、最前部のみの情報源Nを推定する場合を考える。t時点の情報源Nからの予測誤差をwN（t）とするとt∈丁爵とt∈Ts翁ではwN（t）の共分散が異なることに注目する。t∈Ts㌔におけるwN（t）の推定共分散をAwとし、（3．4．5）を最大化する。ただし、t∈丁吾NにおいてもE［wN（t）w長（t＋τ）］二〇（7キ0）と仮定する。手法2と同様、2つに大別される。

1）t∈丁爵の共分散行列を推定情報源の共分散行列に固定する。このときAICはAIC（k）＝（k十MN）1η惚CN十（MsN−k十1）Aw
k−1
十 Σ

＾
＿
wl（t）C一凸wN（t）十（MsN−k十1）r（3．4．11）

t二一MN

十r2十r十const．

2）t∈丁轟の共分散行列を推定する。このときのAICはAIC（k）＝（k十MN）劾44AN十（MsN−k十1）肋44Aw＿ ＿十（MN十MsN十1）r十2（r2十r）十const．

これらのAICが最小になるkをもって信号出現時刻とする。

一一79一（3．4．12）

気象研究所技術報告

第16号

1985

手法3はただ1つの情報源のみを推定するだけなので、SまたはNいずれかの情報源が推定困難な場合にむいており、計算時問は3つの手法の中で最も短い。以上、大別して3種の手法を述べた。
ここで述べた手法は、バラメータ数をゼ・と仮定することによりOrigina1なデータをそのまま用いることもできる。手法3は、手法1〜3との組み合せとして用いることもできる。つまり、片方のモデルによる逆フィルタを通すことと同等である。このホワイト化フィルタの出力に対して再び手法1〜3を適用するのである。これは、手法2についても同様である。このようにあらかじめホワ
イト化フィルタを使用する手法を、仮に手法0と呼ぶことにしよう。
実際的問題として観測量が多成分の時、計算量は膨大になるが、次のように簡略化できる。まず第1成分の事後確率を第2成分の事前確率とし、ベイズ流に第2成分の事後確率を得る。これをr成分まで順次繰り返して最終的な事後確率を得る。このときのAICを最小とする時点kを信号出現時刻とする。この方法は、本質的に共分散行列の対角成分を0と仮定することに相当しているが、十分実用に耐える。実際に適用するときは、第1成分で十分な精度が得られないとき第2成分を…と満足すべき精度が得られるまで順次成分数をふやしていく。
さてこれまで述べてきた信号出現時刻は、等時間間隔でサンプリソグされた観測量をk

用いて、最初に信号が出現した時刻として、定義される。しかし、これは従来のアナ・

1

グ可視記録によって人間が信号出現時刻と判断していた時点とは異なる。例えば、地

一

響

k−1

震波信号の初動時刻としては、正に信号が始まらんとする時点を初動時刻としてアナ↑・グ可視記録から人間が読み取っている。

このように、時間を離散的に見た場合と連続量的視点から見た場合とでは信号出現時

図3．4．3

刻のとらえ方が多少異なる。この両者の視覚的な差異は、出現した信号が大きければ離散時間のモデルにおいて信号が出現した
時点kと連続時間的視点からみた信号出現
時刻との違い。図中矢印の時点（k−1）
を視覚的な初動時刻と呼ぶ・大きい程顕著である。そこで、ここで定義
した信号出現時点を1サンプル時点前にもどした（k−1）時点を、視覚的な初動時刻として定義する（図3．4．3）。

3．4．2

ベイズの識別手法として見た信号出現時刻の決定

前項で述べた信号出現時刻の決定手法を、ベイズの識別手法の立場から見直してみる。
信号出現時点kすなわち（3．4．4）あるいは（3．4．5）を最大化するkは確率過程｛Xt｝により決まる

一80

気象研究所技術報告

第16号

！985

確率変数である。しかし、実際には各kにより逆に確率過程／Xt｝が決まる。なぜなら、情報現Nからの出力／nt／および情報現Sからの出力｛St｝の確率過程がそれぞれ独立に存在し、これら2つの出力端子を接続させた時刻がkであり、その結果として｛Xt／が決まるからである。そこで、これらの出力端子を接続させる事象の全体をAとし、時刻iに端子を接続させる排反事象をAiで表わ
すと、信号出現時刻の決定問題は、情報／Xt／が事象Aのうちどの事象からの情報であるかを識別する問題となる。

いま、A＝IA1，A2，…，Am／、（Ai∩Aj＝φ，iキj）とし、／xt，t∈To／の全体をX、損失行
列をC＝［Cij］、P二（P（A、IX）…P（Am

I

X））とおくと、損失ベクトル1＝（1、」2…塩）は（3．2．2）と同様に表わされる；
1＝CP（3。4．13）

したがって、ベイズの識別方法は

郷初li二lk→X∈Ak（3．4．14）iとなる。ここにX∈AkはXがAkからの情報であると識別された、すなわち信号出現時刻がkであると判断されたことを表わす。

（3．4．14）のベイズの識別手法において、前項の多時点識別手法を見てみよう。前項で述べた種々の手法は、本質的には尤度関数の違いと考えられる。そこで、これらの手法におけるi時点に信号が出現する尤度をP（XIAl）で代表させることにする。前項のAlc（i）を用いるとP（X陶一呵一麦AIC（i）／（3・4・・5）

により与えられる。したがって、3．4．1項の手法は、
吻砿P（X

I

Ai）＝P（XI

Ak）→X∈Ak（3．4．16）孟により信号出現時刻の決定を行ったことになる。これは、事前確率および損失をP（A1）＝P（A2）＝一一二P（Am）

C11−CQ『

Cij＝Cji＝Ce（i＝1，2，…，m）（i＝1，2，…，m，j＝1，2，…，m，iキj）

Co＜Ceとしたときの（3．4．14）のベイズの識別手法と同等になる。ここに、事後確率はP（A、IX）〜P（Ai）P（X

I

Ai）／P（X）

一81（3．4．17）

気象研究所技術報告

第16号

1985

を満たすと仮定した。

次に、白井・徳弘（1979）による信号出現時刻の決定手法をベイズの識別手法の立場から見てみる。あるパラメータ（いまの場合信号出現時刻）を求めるために通常なされるように、（3．4．4）あるいは（3．4．5）式を時間について微分（正確にはいまの場合は離散的なので差分）してゼ・とおき
ズぐ
極値を求めることが考えられる。情報源NおよびS①Nを固定した前項の多時点識別法の手法2−1）

において差分〜AIC（k＋1）＝AIC（k）＝0→信号出現時刻はk時点とすれば、森田・浜口（1984）が述べているように白井・徳弘（1979）と同等である。しかしながら、これをもって白井・徳弘（1979）と横田ら（1981）の多時点識別法が同等であると主張（森田・浜口（1984））することは正しくない。実際、極値は多数存在し、平滑化などの手続きが必要となる。

白井・徳弘らの手法は、1時点しか情報がなく情報源N、S㊦Nが既知の場合のベイズの識別手法
による3・3節の信号検出問題である。
一般化された信号出現時刻の決定は、（3．4．14）のベイズの識別手法により行われることを示した。

すなわち、各信号出現時刻k（クラスAk）の尤度は3．4．1項の手法で求め、その識別は（3．4．14）のベイズの識別手法により行われる。

さて、信号出現時刻の決定を行うために重大な問題として、出現時刻kを決定する区間の設定がある。しかし、これもベイズの識別手法により解決される。いま、識別クラスをA＝｛…，A。＋・，
A・＋2，…，A・＋m，…／と拡張する。そしてこれら各事象の事前確率がP（Ai）＝0（i≦7，i≧7十m十1）となるτおよびmを設定することにより判定区間の設定が行われる。

3．4．3

信号出現時刻の決定誤差

ベイズの識別手法による信号出現時刻の決定手法について3．4．3項で述べた。ここでは、決定された出現時刻に対する誤差の評価について考えてみる。
通常パラメータθの推定誤差を評価する方法は、標本数Nが十分大ぎいとき〉午「（θ一6も）が斬近
的に平均0共分散」一1（洗）1（島）の正規分布に従う（3．3．3項参照）ことを利用して行われる。しかし、いまの場合、信号出現時刻kは（3。4．4）の尤度関数に陽として含まれていない。さらに言えば、信号出現時刻kについての標本は1個しか得られていないのである。
ベイスの識別手法にもどって考えてみよう。識別は損失が最小となるクラスを選ぶことにより行われている。そして、これは各クラス間に棄却領域を設定して検定を行うことと同等であった。このとき、事前確率に基づき棄却領域が設定されたものが、ベイズ解である。したがって、信号出現時刻の決定誤差を、棄却領域の設定法の違いにより識別されるクラスの全体とすることができる。

一82

気象研究所技術報告

第16号

1985

しかし・これは現実的な方式ではない。そこで識別されたAkに対して独自に棄却領域を設定し、この棄却領域でAkと他のクラスとの検定を行い、棄却されないクラスを誤差として採用することにする。すなわち、

レい

AA
＝
A＝A
HH

説説
仮仮
無立
帰対

σ＝

L2

m

キ

k

の検定により、帰無仮説が採用されたクラスの全体を信号出現時刻の決定誤差とするのである。したがって、問題はAkに対して棄却領域をいかに設定するかということになる。
3．4．1項で述べた各手法に対してみてみることにする。まず手法1について考えようこれは、情報源モデルの推定を含めて同時に識別を行う方式であった。そこで、モデルの推定誤差すなわち
π「（a一疏）が平均0、共分数J『1（島）1（島）に従うことを利用することにする。つまり、パラメー

タの推定誤差の共分散行列がX2分布に従うことを利用して、AIC（k）であるいは尤度の推定誤差の棄却領域を設定するのである。しかし、この場合J『1（島）1（島）は末知であることに注意しなければならない。

次に手法2について考えよう。まず手法2−1は、情報源は既知としてその予測誤差の系列により識別を行っている。森田・浜口（1984）は、この手法2−1）と同等の方法で信号出現時刻を決定し、このときの誤差を、既知とした情報源のパラメータの推定誤差を利用することにより評価することを提案している。ここでは別の視点から考えてみる。手法2−1）は、本質的には予測誤差の共分散が既知とした情報源の共分散と等しいか否かを検討している。そこでΣw

lC−1Wtがκ2分布に従うことを利用して検定を行うことにする。このとき、Ak自身についても検定を行うことができる。
次に手法2−2）ついて考えてみよう。これはシステムダイナミックスは既知として予測誤差の共分

散行列を推測することにより識別を行っている。そこで、AkとAiのそれぞれのΣwlC『1Wtの比がF分布に従うことを利用して検定を行うことにする。

手法3−1）については簡単ではない。手法3−2）は、予測誤差に対する仮定は手法2−2）と基本的には同等である。そこで手法2−2）と同等にF分布を利用して検牢を行うことにする。そして手法3−1）についても手法3−2）と同様にF分布を利用して行う方式をとる。

次に、Ak自身についても同様の目安を考えよう。Akは、信号が存在するという条件下のもとで識別されたクラスである。そこで信号がないとした場合のAICとAkのAIC（k）との差を、エント・ビーの意味での信号出現の確度として用いることにする。

誤差に対するより実際的な評価の目安を考えてみよう。Akが識別された後、TN、Tsにっいての予
測誤差の共分散行列をそれぞれ計算し（情報源モデル自身も再計算してよい）、それぞれの情報源モデルからの予測誤差をS①Nについてはt＝k−1，k−2，…、Nについてはt二k，k＋1，…と両方向に計算しながらF分布を利用して検定を行い、信頼区間を超えた時点のt1，t2をもって、信号出現時刻

一83

気象研究所技術報告

第16号

1985

に対する区間推定［t1，t2］とする。これは手法1）〜3）について共通に行うことができ、より実際的な方法である。

参考文献
森田裕一・法口博之、1984二自己回帰過程による地震波初動の自動検測とその信頼区間、地震2、37、
281−293。

白井克彦・徳弘一路、1979：地震波初動の自動検出、地震2、32、141−148。

横田崇・周勝奎・溝上恵・中村功、1981二地震波データの自動検測方式とオソライン処理システム
における稼動実験、地震研究所彙報、55、449−984。
Allen，R．V．，1978：Automatic

earthquake

recognition

andtimingfromsignaltrace．Bull．Seism．Soc．

Am．，68，1521−1532．
Anderson，K．R．，1978：Automatic

Analysis

and

along

of

microearthquake

network

data，Computer−Aided

SesmicDiscrirination．ed．C．H．Chen，Elsevier，Amsterdma，159−175．

Stewart，S．W．，W．H．K．Lee

quakesanalysis

the

San

and

J．P．Eaton，1971：Location

Andreasfault

systeminandcentralrea1−time

detectionofmicroearth−Califomia．Roy．Soc．New

Zealand，9，205−209．

3．5

地震波信号における相の同定

時系列データの信号検出手法については、2者強制選択問題として3．3節で述べた。検出したい信号が複数個ある場合についても、基本的には3．2節で述べたBayesの決定手法により解決される。
このような複数のクラスの信号検出問題の一つとして、地震波信号がある。

地震波信号は、P波、S波等の混合波からなり、特に震源を決定する上で重要なことは、これらP波、s波の識別である。2者強制選択問題としての地震波信号検出は、通常次の3段階の手順により
行われる。

i）

地動雑音からP波信号の検出

ii）

地震波信号の終了検出iii）

地震波信号の中からs波信号の検出

地震波信号は、発震機構や波の伝播経路の地殼構造により非常に多様な様相を呈する。例えば、震央からある方位の観測点では、P波信号が非常に微弱となり、S波信号しか識別できないようなことがある。このような場合、単に最初に検出された信号をP波信号とする上述の手順では、S波をP波と誤認し、震源決定に重大な悪影響を及ぽすことになる。

しかし・P波・S波の情報モデルを設定することは容易ではない。なぜなら、これらの波は、発震機構とその伝播経路により決まるもので・情報源モデルはこれら両者をあわせた伝達関数により表

一84

気象研究所技術報告

第16号

1985

現されるからである。そこで、これらの波を、振動方向による特徴パラメータにより識別することを考える。

地震波の振幅は、P波の場合は上下動成分に、S波の場合は水平動成分にそれぞれ卓越することが多い。またP波は震央方向に振動する波である。このような特徴を利用し、P波、S波を識別することおよび震央方位角を推定することは古くから行われている。最近では、竹山ら（1977）が、計算機に接続された観測網で、また、彼らと類似の手法により中村ら（1984）が、1観測点のみでP波、o

S波の識別および震央方位角を推定することを試みている。一方、Flim（1965）は、3成分波形を
particle

motionによりP波的、SH波的、SV波的な波に分離するためのpolarization

案し、Montalbetti

filterを提andKanasewich（1970）は、この手法により遠地地震波形の中に含まれる構造的な波形の抽出を行った。

ここでは、振動方向の特徴パラメーターの抽出のために、Flinnらによるpolarizationfilterの概念を利用した。以下に、地表面に入射したP波、S波の振動方向（3．5．1項）、主成分分析の立場か
ら再編成したpolarizationfilter（3．5．2項）、震央方位角の推定とその推定誤差（3．5．3項）、そして最後に、P波、S波の特徴パラメータによる識別基準について記述する。

3。5．1地表面でのP波、S波の振動方向
地震波信号は、通常、地表面での振動を観測することにより得られている。P波、．S波の識別を振動方向による特徴パラメータで行う準備として、これらの波が地表面に入射したときに観測される見掛けの振動方向について簡単に記述しておく［詳しくは、例えばEwing

etaL（1957），Aki

andRichards（1980）］。

P波は伝播方向に振動する縦波であり、S波は伝播方向に対して直交する方向に振動する横波で、鉛直面内の振動成分（SV波）と水平面内の振動成分（SH波）に分けられる。いま、水平成層構造
で異方性を持たない場合には、それぞれの直達波が地表面に入射する直前の振動方向は、図3．5．1のようになる。したがって、これらの波のparticle

motionは、入射角（ここでは、水平面から下向きに測った角度とする）と、震央方向（北から時計まわりの方向に測った角度）に対する情報を持っている。

地表面（自由表面）への入射問題として、最も簡単な半無限弾性体の自由表面に平面波が入射した場合を考えよう。SH波が入射した場合は振幅が2倍に増加されるだけで何ら特別なことは生じない。しかし、P波またはSV波が入射した場合には、反射P波、反射SV波が生成されSH波の入射の場合とは事情を異にする。

波が入射する角度をf、反射P波、反射SV波のそれぞれの角度をfp、fsとする。ただし、これらの角度は、水平面から下向きに測るものとする（図3．5．2）。またP波、S波の速度をVp、Vsとする。このときSnel1の公式より

一85

気象研究所技術報告

P

u

Ne

第16号

1985

／R

R

SV

u
f・iH

u

N

fa

fs

E

R
SH

u

fp

A

P
P・rSv

N

図3．5．2Sv

半無限弾性体の自由表面入射の概念

図。fは入射角（P波入射のときfp、
SV波入射のときf、）、f。は見かけの入図3．5．1 地表面へ入射直前のP波、S射角を示す。

V波、

SH波のparticle motionの概念図。
Rは震央方向を示し、θは方位角で
ある。

60sf

V

60sfp
Vp

60sfs
Vs

（3．5．1）

ただし

v一｛ll：認翻となる。

P波が入射した場合の見かけの振動方向f、は、上下動振幅Uと水平動振幅Hとの比として簡単に求めることができ
伽fa＝一60渉2fs
ただし

U／H＝如％faとなる（図3．5．2参照）。fsは（3．5。1）より聡一・・ガ1（蜘6・s鼻〉

一86＿（3．5．2）

気象研究所技術報告で与えられる・このとぎ、Vp＞Vsより、常に、60s

第16号

1985

fs＜1である。（3．5．2）を変形すると、より直接的な関係式

fa＝2fs一一2で表わされる。

次にSV波が入射した場合を見てみる。このとき（3．5．1）よりVp
60sfp＝ 60sf≧1
Vsとなる場合がある。等号が成立するときの入射角f。＝60s−1（Vs／V，）はcriticalangleと呼ばれている。f＜f・の場合には、いわゆる全反射となりphaseshiftを生じる。このため水平方向の変位は、f＝45。でゼ・となり、f＜45。では見かけの振動方向は、入射SV波の振動方向とは逆方向の振動となる。

すなわち、振動方向はP波の振動方向と同一の象限となる。
さて、f＞f，の場合を考えてみよう。この場合の見かけの振動方向は簡単に求めることができ

伽fa＝2鰭fp／（1一伽2f）（3．5．3）となる。ここで、振動方向として、U／Hをとっているので、（3。5．3）からすぐ解かるようにf、＜0である。すなわち、P波入射の場合とは逆の象限で振動している。これに対し上述したf＜45。の場合は、fa＞0である。

3．5．2主成分分析の立場から見たpolarization

filterの構成

時刻tにおける水平2成分（通常は南北、東西）、上下動成分の振幅値をN（t）、E（t）、U（t）とし、x（t）＝（N（t）E（t）U（t））

描く軌跡は、particle

とおく。t＝t。、t。＋1、…、t。＋nに対するx（t）が3次元空間内で

motionと呼ばれている。

particlemotionの主たる振動方向を求めることを考えてみよう。これは、z・（t）＝xてt）e・という変換を考え、ΣIZ・（t）／2→maxとなるベクトルe・を求めることにより実施される。但し、Il

e・II2二1とする。同様にして、第2番目、第3番目の振動方向を求めることが考えられる。ただし、各振動方向は直交するようにとる。このような手順で各eiおよびZi（t）を求めることは主成分分析と呼ばれており、これはx（t）の共分散行列の固有値問題として解決される。
いま、簡単のためにT＝／1，2，…，n｝とし、x（t）、（t∈T）は同じ波群からの出力と考えよう・

ここで、（x）tの標本共分散行列をR、Rの期待値をCとする
C＝E［R］＝E［x（t）x

ここにx＝（NE（t）］二xx（3．5・4）U）・とする。以下、共分散行列Cについて考える。このことによって一般性は損

一87一

気象研究所技術報告

第16号

1985

なわれない。

Cは正値対称行列なので、次のように分解される。
（3．5．5）

C二丁ΣT

ただし

Σ一〔詫〕
λ1〉λ2＞λ3＞0

丁二［ele2e3］、TT＝1e、、e2、e3は固有値λ、、λ2、λ3にそれぞれ対応する固有ベクトルである。λ、、λ2、λ3に対応する主成分をZ＝（Z、Z2Z3）

また、主成分得点（因子スコア）をZ（t）二（Z、（t）Z2（t）Z3（t））

z＝T

とすると、

x
（3．5．6）z（t）＝Tx（t）

E［z（t）z（t）］二zz

で表わされる。また、（3．5．7）＝Σ

となる。すなわち、各主成分の分散はそれぞれ対応する固有値と等しい。各主成分の分散の全体に占める割合は、寄与率と呼ばれ
γi二λf／trace（3．5．8）Cで表わされる。またγ（」）＝Σ質は累積寄与率と呼ばれている。

i冨1

いま、主成分の標準偏差に境界面を設置すると、x（t）のparticlex

Cx＝1

motionは、（3．5．9）の楕円体で表現されることになる。したがって主たる振動方向は、第1主成分に対応する固有ベク
トルで表わされ、全体の波形を代表する程度は、第1主成分の寄与率で見ることができる。また、第2、第3主成分の影響度は、寄与率範、％を用いて評価できる。
さて、各主成分を構成するもとの変量の割合は、効果ベクトルと呼ばれており、基底ベクトルを［eNeEeu］とすると
xcj＝竹［NeN

EeE

Ueu］

一88

ej

（3．5．10）

気象研究所技術報告

第16号

1985

で定義される。また、主成分得点に対応するものは、
xc」（t）＝必［N（t）eN

E（t）eEU（t）eu］窩

（3．5．11）で定義される。この効果ベクトルを用いて、

ゴ
又き＝ΣXci一・ji（3．5．12）

又き（t）二ΣXc（t）
i＝1iを定義する。ただしj二1、2、3である。

これは、累積寄与率に対応する効果ベクトルである。j＝1の場合が、Flimらのpolarizationfilterに相当する。我々は（3．5．12）で定義される量をpolarizationfilterと呼ぶことにする。
実際の波形へこのfilterを適用する際は、共分散行列が時間とともに変化していることを考慮して、

例えば共分散行列Cを
R（t）二一Σx

1 M（t十i）x（t十i）
2Mi一一M（3．5．13）で置き換えればよい。

3．5．3

震央方位角の推定

震央方位角の推定に、サンプル時点毎の3成分の振幅値を用いる手法が採られることが多い。しかし、各時点毎の値は非常にバラつくため、上下動方向の振幅を重みに採用して平滑化することにより計算されている［例えば竹山ら・（1977）、中村ら（1984）］。この方法では、震央方位の推定誤差を得ることは困難であり、さらに最小2乗法で求めた場合は、真の振動方向とはならない。また、中村ら（1984）は、S波群についても、P波群と同様の計算により震央方位を推定しているが、これは3．5．1項で述べたように常に正しいとは限らない。すなわち、SH波の場合は±90。、SV波の場合は全反射の領域では正しいが、そうでない領域では一180。ずれることになる。

ここでは、固有ベクトルを用いることにより、P波、S波群からの震央方位の推定法およびその推定誤差について記述する。

（1〉P波群からの震央方位の推定
P波群からの震央方位は、パーティクルモーショγが入射方向に振動することから、第1主成分の固有ベクトルがP波群の入射方向を表わしている。また第1主成分が全体の波群をどれだけ代表しているかは寄与率で測られる。第1主成分の固有ベクトルをe、二（e、N

e、E

e、u）

とすると、北から時

計回り方向で表わした方位角命は命二伽一1（e・E／e・N）一δqUπ

一89（3．5．14）

気象研究所技術報告

第16号

1985

で与えられる。ここで、

鑑一〔1：盤とする。e、u二〇のとき方位角は不定である。

（2）S波群からの震央方位の推定
S波群からの震央方位の推定はP波群からすでに述べたように単純ではない。すなわち、SV波が全反射領域でphaseshiftを生ずるからである。さらに、SH波だけしか存在しない場合は、振動方向に直交する直線上に震央が存在することは言えるが、そのいずれの方向であるかを規定することはできない。したがってここで述べる推定法は、あくまでも目安にすぎない。

さて、S波は、SH波、SV波に分解して考えることはできるが、実際問題としてはその両方を含んだ形でS波は観測されている。そこでSH波、およびSV波を同程度含んだ場合について考えよう。

まず入射角がCriticalangleよりも大きい場合について考える。この場合は、第1、2主成分の固有ベクドルによって張られる平面は、震央方位に対して垂直で、かつ震央の位置とは逆方向に傾いている。したがって、震央方位角θsは、第3主成分の固有ベクトルe3＝（e3N

e3E

e3u）

によりP波群と同様の式：
醜＝伽一1（e3E／e3N）一δ馬Uπ（3．5．15）で推定される。

入射角がcriticalangleのときは、上下動成分はゼ・となり、S波は、一見真下から入射したと同様になる。criticalangleから450の間も上式で求まる。

さて、入射角がcriticalangleより小さいときは、phaseshiftにより、S波の振動面は、震央の方向に傾いた面となる。したがって震央方位角は、＆二勧一1（e3E／e3N）一（1一δ免U）π（3．5．16）となる。

（3．5．15）と（3．5．16）は、第3主成分の固有ベクトルの上下動方向の成分の正負により震央の位置が異なる点である。これはSV波の地表面における効果によるものである。したがって、もしもSV
波しか存在しない場合は第1主成分の固有ベクトルにより、入射角が45。以上のときは（3．5．16）で、45。以下のときは（3．5．15）でと逆の式で求めれぽよい。

いずれの式で求めるかは、第1主成分、第2主成分の寄与率で判別すればよい。また震源に対して何の情報も持たない場合は、SV波に対応する固有ベクトルが寝ていれば（3．5．15）で、立ってい90

気象研究所技術報告

第16号

1985

れば（3．5．16）で求めればよい。しかし、これらはあくまでも目安にしかすぎない点を注意しておく。

（3）方位角の推定誤差推定された方位角の誤差について考えてみる。震央方位を示す固有ベクトルを含む切断平面を水平面に投影した2次元平面z＝（z、z2）で考える。この平面において、z、を推定された震央方向にとることにする。このとき水平面に投影された固有値をそれぞれ6、6とすると、

唄鵬〕（3．5．17）となる。

方位角が正しく推定されているか否かは、上式が正規性の仮定のもとに成立しているか否かである。すなわちZ、とZ2の独立性の検定問題になる。

いま相関係数をρ、標本数をnとすると、相関係数ゼ・の仮定のもので、〉冠ρ／〉マは自由度n−2のt分布に従う。t分布より求まる100α％の危険度での相関係数をραとすると、

zガー〔のρα（巧σゑ

ρα（巧（，2（3．5．18）2
σ2となる。この（3．5．18）から求まる角度を震央方位の推定誤差と定義すると、推定誤差θ。は2

屍一⊥伽一・2ρ襯（3．5．19）61一品で与えられる。6二6のときは、震央方位角は不定である。

3．5。4P波、S波の特徴表現

観測されるP波、S波は、すでに述べたように地殻構造と震源パラメータにより非常に多様な様相を呈する。3．5．1項で述べた最も単純な半無限弾性体の自由表面への入射の場合ですら、振動方向は入射角度により大きく異なる。すなわち、地震波信号の振動の様相は、震源パラメータが決まった後の比較判断により正しく評価されるべきものである。しかし、震源位置を正しく決定するためには、震源位置が未知の状態でP波、S波の識別を行わなければならない。したがって、正しい震源位置は、震源位置とP波、S波の識別とそれらの初動時刻のフィードバックにより評価され決定される。

震源パラメーターが末知な状態での第1義的なP波、S波の識別として、従来から行われている
ように、上下動成分が卓越すれぽP波、水平動成分が卓越すればS波という方式を利用しよう。これは、入射角が立っている場合は正しい。そして一般に地表に近づく程入射角は立ってくるからで

一91一

気象研究所技術報告

第16号

1985

ある。

さて、P波のparticlemotionは、第1主成分により代表され、S波はSH波、SV波を同程度に含むとすれば第1、第2主成分により代表される。したがって、P波的な振動を上下動で見ることにすると、P波度を表現するものは効果ベクトルの上下動成分、S波度を表現するものは効果ベクトルの水平動成分とするこ，とができる。そこで、第i主成分の固有ベクトルの上下動方向の方向余弦をeiUとし、P波度、S波度を、
Dゴ（P）＝Σ箕eぞu（3．5．20）t＝1

まDゴ（S）＝Σ％（1−eぞu）（3．5．21）t＝1で定義する。ただし、j二1または2。ここで、上下動方向および水平動方向の合計の効果を表現する
ためj二1または2とした。Dノ（P）、Dゴ（S）はparticlemotionの形状のみに依存している。次に固有値の大きさに依存する量として
D蓬（P）＝Σγ抗∫eぞu（3．5．22）i＝1ゴ
D蓬（S）＝Σ％λ，（1−eぞu）（3．5．23）∫＝1を定義し、これらをそれぞれ固有値依存なP波度、固有値依存なS波度と呼ぶことにする。D差（P）、
D差（S）は、それぞれ上下動方向、水平動方向の振幅のエネルギー的な効果を表現し、これらの平方根
は振幅的な効果を表現したものである。Dゴ（P）、D，AI二D」（S）およびD蓬（P）、D蓬（S）の比をそれぞれ（P）／D」。（S）（3．5．24）

AE二D蓬（P）／D蓬（S）（3．5．25）

Dpls二如盟一1A（3．5．26）とおき、を定義し、DplsをP／S波度と呼ぶことにする。ただし、A二AiまたはAE。j二1のとき、A【＝AEとなり、Dplsは第1主成分の見かけの入射角を表わす。

A（二AIまたはAE）あるいはDplsを用いた第1義的なP波、S波の識別手法は

A＞β（またはDpls＞如π一1β）→P波
A＜β（またはDpls〈勧一1β）→S波

一92一（3．5．27）

気象研究所技術報告

第16号

1985

となる。

震源が求まった場合について考えよう。この場合、3成分の振幅を直接計算することによる比較判断も行えるが、ここでは方位角を利用した評価量を導入する。

いま求まった震源位置に対する方位角ベクトルをac、P波群およびS波群から3．5．3項の手法で推定した方位角ベクトルをap、asとする。地殼構造に異方性がなければ、これら方位角ベクトルは一致しなければならない。そこで、P波の確度r（P）、S波の確度r（S）をG

r（P）二a6ap（3．5．28）r（S）二a6as（3．5．29）と定義する。これらの絶対値が1．0に近いときacと、apあるいはasは一致し、最も確度は高い。また、r（P）、r（S）の誤差は、（3．5．19）式で与えられる方位角の推定誤差を持つ方位角ベクトルとの内積で与えられる。

ところで、S波群からの方位角の推定は、3．5．3項で述べたようにSV波に依存しており、その推定式は入射角により異なる。そこでSH波に相当する固有ベクトル（上下動方向の方向余弦がほとんどゼ・）の水平面への投影方向を持つ大きさ1のベクトルをaSHとし、これとacとの関係を表わす
量を考えよう。aSH、acは、直交することが期待されるので、r（P）、r（S）に相当する量として、SH波の確度r（SH）を、r（SHl）＝1−

I

a6asH

I

（3．5．30）で定義する。r（SH）一・のとき最も確度が高し・。ただし、一・≦｛r（P）・r（S）／≦・に対し・・≦
r（SH）＜1である。

次に、震央は不明であるがP波であることは正しく識別された場合のS波の確度を考えてみよう。

P波の第1主成分の固有ベクトルをep、SH波に対する固有ベクトル（第1主成分と第2主成分の固有ベクトルのうち、水平に近い方）をeSHとすると、これらは直交することが期待される。そこでP波に対するSH波の確度をr（SHr（SHと定義する。r（SH

I

I

I

p）で表わし、

p）＝1−I

eβesH

I−l

aβasH

l

（3．5．31）p）＝1のとき最も確度が高い。これは、対応するSH波が充分に含まれている時には非常に有効であるがそうでないときはほとんど意味をなさないことに注意しなければならない。

以上、相の同定のための量および、P波、S波を表現する量を定義した。ここで定義した量は最終的な震源位置を決定する上で重要な評価量になる。そして、これらの量の信頼度は寄与率で評価される。また、これらはpolarizationfilterからの出力により計算することができ、これを用いてP波

信号あるいはS波信号の検出を行えることを付記しておく。

一93一

気象研究所技術報告

第16号

1985

参考文献
竹山一郎・柴田武男・佐藤

方位の推定

馨、1977：地震波形の処理に関する研究（IV）

ディジタル処理による震央、気象研究所報告、28、147−157。

中村豊・斎藤明夫・橋本民雄、1983二P波とS波の識別、地震学会予稿集、No．1。
Aki，K．and

P．G．Richards，1980：Quantitative

Seismology−Theory

andMethods，W．H．Freeman

andCompany，USA．
Ewing，W．M．，W．S．Jardetzky

New

and

F．Press，1957：Elastic

Wavesainusing

direction

LayeredMedia，Megraw−Hill，

York．

Flim．E．A．，1965：Singnal

analysis

rectilinearity

and

of

particle

motion，proc．L

E．E

E．，53，1874−1876．

Montalbetti，J．F．and
polarization

3．6

E．R．Kanasewich，1970：Enhancementofteleseismicbodyphaseswithafilter，Geophys．J．R．Aster．Soc．，21，119−129．

おわりに信号検出とその出現時刻の決定問題を、ベイズの識別手法として一般化して述べた。図3．6．1〜図3．6．5に若干の解析例を示す。本章では信号出現の場合の問題しか述べなかったが、信号がある期間だけ出現する現象においては、信号出現の判別と信号終了の判別を行わなければならない。しかし、

これら両者の問題は本質的に同等である。信号終了の判別の実際的手順としてはi）信号状態からノイズ状態への変化を識別する手法と、ii）ノイズ状態か否かを判断する手法がある。i）の手法は本章で述べたクラスNとクラスsの仮定を入れ替えることにより実行される。ii）の手法は、クラスNを既知としてクラスNか否かの識別を行えばよい。
地震波信号処理に特有な問題として震源決定の問題がある。この問題の最良解決手法は、3．5節の量を用いて、震源位置と地震波信号の比較判断を繰り返せばよい。このためのアルゴリズムのフ・・一を図3．6．6および3．6．7に示す。これまでの単に初動時刻を用いただけの震源計算では大きく誤っていた場合も、このアルゴリズムを適用することにより飛躍的に信頼度をますことが確認された。
現在、時系列解析は、時変システムの解析、多数情報源モデルの分離、非線型システムの解析等、新たな展開をしつつある。また、このような解析に不可欠な、母集団のパラメータの確率的構造を明らかにする手段に関しての研究も行われつつある。一方、地震の発生のような事象系列データの解析モデルについての研究も行われるようになった。今後、事象系列＋時系列の新しいデータ解析手法の開発により、更に効率のよい地震活動、地震波信号、地殼変動データを含めた総合的な処理も期待できるが、これには装置規模等の問題を含め未解決な部分があるので、ここではまとめない。

一94

気象研究所技術報告

第16号

1985

画

儲S

占

100

O

0

1sec，

I

｝

lサ

舶川口内
ヤ

嫡
1

サ

P

0

15ec．

図3．6．1 地震波信号に3．4．1項の手法1を適用した信号出現時刻の決定例。矢印P、Sはそれぞれ決定されたP波、S波の初動時刻を示す。

一95

気象研究所技術報告

第16号

1985

荊胃f
P置一model

AICB

fS

一一

ナP

N−model

㎞

500
100

0

0

1sec．

0

図3．6．2 地震波信号に3．4．1項の手法3−2を適用し、第1成分の事前確率に…とベイズ流の手順で決定した信号出現時刻の例。矢印P、Sは、それぞれ決定された
P波、S波の初動時刻を示す。

一『96一

気象研究所技術報告

第16号

1985

O Z一＞﹂﹂のZ田〇一く庄﹄鳳Q山氏の

oO

も

撃爵ず

も駕

↑
↑

ち

，

もψ

〜Y蒔

1

、

0（）
0

5SE（二

捨
うρ⊃

1午ーP

図3．6．3 地震波信号の信号検出とその出現時刻の決定に対比してスペクトルの変化を示す。矢印P、Sは、3．4．1項の手法2−2により決定されたP波S波の初動時
刻を示す。

一97

／

気象研究所技術報告

第16号

1985

榊一〜1
u
D
口．00

5．00

10．00

15．00

20．00

25．00

30．00

35．00

40．00

45．00

50．00

与5．00

60．00

〜0．00

25．00

ヨ0．00

ヨ5．00

40．00

45薗00

50．00

55．00

60．00

N
S
E

佃

W
u
D
0、00

5．00

10．00

15．00

図3。6．4 Polarization filterの適用例。

横軸は時間単位：秒である。上段は3成分の源波形を、下段はporalization
filterからの出力。ただし、固有値および固有ベクトルは1秒間の波形データの移動平均により求めた。

98

気象研究所技術報告

第16号

1985

v

E

孔

、1糊

W
u

妙瑠

D
oo

4．00

8．00

12．00

15．00

20．00

24。00

28．00

32．00

36．00

40．00

44．00

48．00

ゆ

P／S
o
面︒卜

oo

4，00

8．00

12．00

】5．00

〜0り00

24、00

23．00

32，00

36．00

40．00

44．00

48．00

8・00

12．00

15900

〜0，00

24．00

28，00

32。00

36．00

40，00

44．00

48．00

DE（S》

O︐O

OgOO

4。00

CONTRIBUTlON

ロ

0．00

0
⑩

4．00

8．00

12．00

15．00

20．00

24．00

28．00

32．00

36．00

40．00

44．00

48．00

12．00

15．OD

20．00

24。00

28．OG

32．00

36．00

40．00

44．GO

4B、00

AZIMUTHIPl

o①
0．00

4．00

8．00

図3．6．5 地震波信号の相の同定のための各種評価量の例。

横軸は、時間（単位：秒）。上3段は、細実細が3成分元波形、実線がpolari−zation filterからの出力。その下は順に、P／S波度、固有値依存なS波度、第1主成分の寄与率、そして最下段が方位角である。方位角は、P波として求めてある。したがって、S波以降乱れているのがわかる。各種量の中で黒く塗られている部分は、寄与率80％以上の部分に対応している。なお、固有値および固有ベクトルは2秒間のデータの移動平均により求めてある。それぞれ、必要な部分にフラックが立てられることがわかる。

一99

気象研究所技術報告

第16号

1985

各観測点毎に

START
next

仮P波検出
noise
START

P−time

検測

信頼区問，初動方向＊P／S

各観測点の検測情報の編集，結合

、

度

＊寄与率

観測点の区域分け

＊方位角

noise

振幅，周期
noise

震源計算に使用する観測点の選択，結合，分離（順方向サーチ）

震源計算

最大振幅，周期F一もime0−Cの検定S㊥

P＊誤差の評価F−timeの検定

最大振幅の検定（逆方向サーチ）

Ma8nitudeの検定

仮S波検出＊方位角の検定
nonext

exist？

地震の統合，分離yes
S−ti面eOK唖

検測

信頼区問，初動方向

＊P／S

＊再検測

度

＊寄与率
＊方位角

＊P波に対するS波の確度
振幅，周期

図3．6．7 観測点毎の解折結果を統合編集

し、震源を求める反復処理のフ

next

・一チャート。再検測は、必要に応じて、観測点毎の情報を再検測することを意味する。

図3．6．6 地震波信号の1観測点の信号に対する自動検測処理のフローチャート。＊印は、3．5節の量に対応する。

本章を閉じるにあたり、ていねいな指摘・議論等をいただいた地震火山研究部第1研究室伊藤秀美研究官に感謝いたします。

一100一