1977年3月

第17一号

国立防災科学技術センター研究報告550．34

強震言己録の 処理に関す る一考察
その（1）
木

下

繋

夫＊

国立防災科学技術セソター第2研究部○皿alM1ethodforProcessi皿gEarthq11akeof

Stmng−motion

the

（I）

A㏄e1emgram

By

ShigeoKi110shita

亙α舳岬α肋肋伽θθγ切Lα6・γαむ・〃，S…〃肋・θ舳ん刀伽・｛・〃N棚o伽ほ脇㈹ん0θ械θけoγ1）づ8α8乏〃P榊θ〃づo侃，No．棚9−1，
K〃仇α㈹，Sαん〃α一仇〃α，N舳α外σ舳，肋㈹肋一ん舳，300−32

AbstractThis

reportstrong−motion

deals

with

earthquake

and

its

the

autoregressi▽eamethodfora㏄elerogramapplications．This

methodspectrum

is

based

estimation

representing

by

using
on

a

both

the

characteristics

smal1number
themethod，assumingofprediction

such

a

of

the

parameters，

error刮terrecorded

anda㏄elera−tion is a sample function from a second order，covariance stationary stochastic
process．Parameters which this method requires are obtained by a fol1owing
procedure：

1）
2）

The

digitized

a㏄eleration

with

a

interval．

The

linearprediction

methoddigitized

data，and

this

data

of

the

uniform

isobtainedisappliedisbyinsamplizing
theexpandeddominant

by

the

a

recordedduration

uncorre1ated

wave

of

the

randomvariab1es．

3）

The

covariancematrixprediction

error

exPansion b㏄ome required Parameters．
There are two basic schemes to app1y these

and

the

matricesparamerters．One

is

of

the

this

app1i−

cation in the time domain．For this purpose，the prediction error filter whichis constructed directly by these parameters，is used．This丘1ter is driven by
uniform random numbers between zero and one，toa㏄eleration waves with the same sPectra1density

generate artificial
matrix rePresenting

earthquake
an origi−

nal

app1y

earthquake

parameterscoe舶cientsinacceleration

the

are

computedsPectral density matrix．
In one dimensional
theconcept

systemandrecord．

frequencyofvirtualthecalculation

the

other

hand，indomain，martices

whoseelements

from

0n

these

case，various

groundユayers，theofparameters乱nd

are

order

are

used

to

these

autoregressive−

to

the

aPP1icationsareidenti丘cationofamulti−degree−of−freedom

spectrumareimportant．Parameters，quasi−response＊耐震実験室一g1一considered．In

estimate

particular，

国立防災科学技術セソター研究報告
which

are

equivalentvirtualiscoe冊cientsoftocoe箭cients

reHection

groundanlayers．Supposing

constructedbyrandomexpansionthatgroundequivalentmulti−degree−of−freedom

equationa㏄elaration

which

appears

Artificial

acce1eration

waves

are

to

the

used

make

out

a

in

this

waves
the

which

by

scheme

is
is

are

are

of

generatedquasi−resPonse

scheme，i．e．

asasystemthatsystemisthepredicted
by

the

by

variables，areregardedestimated

denominatorrandommulti−reiectionispossib1e

whichsystem

1977年3月

uncorrelated

mathematicalvariables，itestimate
istic

in

第17号

at

to

site．Thethecharacter−

spectra1density

function．

theso1ving

each

which

used

prediction

error

filtersPectmm．

まえがき
ボーリング資料やPS検層資料に基づいて，種々の数値計算から得られる地震動はいわゆ
る

あるべき

地震動である．これに対して，強震計記録は実際にあった地震動である．

しかしながら，前者の各資料は利用され易い形で蓄積されているのに比較して，強震計記録

は大地震以外ほとんど整理されていない．中小地震においてはその記録紙が残されているだけのものが大部分である．

そこで，この報告では強震計記録をバラメータ化し，さらにこれらから原強震計記録の特徴を引き出す一方法について考察する．これが可能となれば，強震計記録はその特徴を表わ

す少数バラメータを数値化されたデータの代りに保存すれぼ良いことになり，特に中小地震の処理に有効である．もちろん，大地震における強震計記録の数値化されたデータは欠くことが出来ないが，強震計記録の今後の増加を考えると，このような処理法が遠からず必要となるであろう．

本報告では，なるべく簡単に強震計記録の特徴を把握することを目的として，記録の主要動部のみを部分的に定常と仮定して扱う．すなわち，強震計記録の特徴を時間領域及び周波数領域において耐震工学的に応用出来るようにバラメータ表現することを目的とする．

時間領域における敢扱いにおいて，強震計記録のパラメータ化はこの方法の最大の特徴を示す．すなわち，従来の人工地震波発生法は，記録波形のもつ連続スペクトノレを線スペクト

ノレの和で近似する三角関数系による展開法がほとんどである．したがって，展開項数はかな
り多くなる。これに対して，連続スペクトノレ構造をもつ予測誤差フィノレタをバラメータから構成して，直接利用する本方法ははるかに少い展開項数で人工地震波を発生出来る．さらに，一次元から相関のある三次元人工地震波までを系統的かつ容易に実現出来る．

1．強震計記録のパラメータ表現
ここでは強震計記録を三次元記録と見なして，その主要動部をパラメータ表現する．また，

これらのバラメータを使い，時間及び周波数領域において原記録の特徴を再現する一方法を
考察する．

一92一

強震記録の処理に関する一考察（その1）一木下

1．1

パラメータ表現

強震計記録の主要動部を標本化時間△Tで標本化したベクトル時系列を｛x（n・△T），n＝1，2，…，N｝とする．ベクトノレまたは行列の転置をプライム記号で表示すれば
x（n・△T）＝［x。（n・△T），x。（n・△T），x。（n・△T）1（1）であり，対応するベクトル値確率変数はX（n）＝［X1（n），X2（n），X3（n）］（2）

で表示される．以後の数学的取扱いを容易にするために，このベクトノレ値確率変数列に対していくつかの仮定を設ける．すなわち，j，k＝1，2，3；m，n＝1，2，…，Nに対して次の三条件を要請する．
（a）

E［Xj（n）］＝O（b）

E［Xj2（n）］＜∞（c）

E［X（m）・X（n）］＝R（m−n）

ここで，R（m−n）は共分散行列であり，E［X1は確率変数Xの期待値とする．（a），（b）は平均零の二次過程を，（C）は共分散定常をそれぞれ仮定している．

強震計記録を予測表現するために，部分空間皿、≡｛X（k−n），n＝1，2，…，1｝によるX（k）
及びX（k一（1＋1））の前向き及び後向き線形予測X（k1I1）＝ΣF（1，n）．X（k＿n）

（3）

n＝1

1

X（k一（1＋1）ll）＝ΣB（1，n）・X（k−n）

（4）

n＝1
を考える．予測行列｛F（1，n），B（1，n），n＝1，2，…，1｝は（3），（4）式の予測誤差を各々とX（kl1）＝X（k）一X（k

l

l）

（5）

X（k一（1＋1）l1）＝X（k一（1＋1））一X（k一（1＋1）11）（6）

すれぼ，下記の前向き及び後向きYule−Wa1kerの方程式の解である．
E［X（k

l

l）．X

1

（k−n）］＝R（n）一ΣF（1，m）・R（n−m）＝O，

n＝1，2，…

，1．

（7）

m＝1
E［X（k一（1＋1）l

l）・X

1

（k−n）］＝R（n一（1＋1））一ΣB（1，m）．R（n＿m）＝O，

m＝1
n＝1，2，…

，1．

（8）

ここで，0は三行三列の零行列とする．また，前向き予測誤差の共分散行列は
Σ（1）＝E［X（kl1）・X1（k＾＝E［（X（k）一X（kll）］l1））．X1（k）］＝R（O）一Σl

F（1，m）・R（一m）

（9）

In＝1

となる．また，皿。と皿1。、≡｛X（k−n），n＝1，2，…，1＋1｝において，Yu1e−Wa1kerの方程式から次の漸化式が得られる（P．Whitt1e，1963）．

一93一

国立防災科学技術セソター研究報告第17号

F（1＋1，n）＝F（1，n）＿F（1＋1，1＋1）．B（1，n），B（1＋1，n＋1）＝B（1，n）一B（1＋1，1）・F（1，n），

1977年3月

n＝1，2，…

，1．

（10）

n＝1，2，…

，1．

（11）

この報告においては，強震計記録を工学的に応用するための必要バラメータとして｛F（n，n），B（n，1），n＝1，2，…，p｝及びΣ（p）を選ぶことにする．すなわち，強震計記録主要動部は図1に示すようにフィノレタΦに入力
x（klP）を加えたときのフィノレタの出カであX（klρ）

fiけer亜るとする・ここで・1・（・，・）・・（…）・・一・
2．

・P｝がフィノレタΦの雛を・Σ（・）が

入カX（k

l

Fig．1．A

計記録主要動部の発生過程がモデル化され↑郵（k同

三（吋（帆・（｝入・斗・、・コ

p）の確率的特徴を各々記述している．すなわち，これらのパラメータで強震

f

X（k）

図1強震波モデル

mode1ofthe

Cele「atiOn

WaVestrong−motion

ac一る．また，このバラメータ抽出の方法ではnを増すごとにX（k）のもつスペクトノレ構造はフィノレタΦに逐次吸収され，最終的にはフィルタΦはX（k）のスペクトノレ包絡をその特性として持つこととなり，X（klp）㍉二は残りのスペクトノレ微細構造のみをその特性として持つこととなる．すなわち，X（klp）は，その共分散行列Σ（p）ともども三次元白色雑音に漸近し

てゆく．

1．2

時間領域における表現

1．1で定めたバラメータから，原強震計記録主要動部とスペタトル包絡特性を等しくする

時間波形を作るためには，フィルタΦを具体的に構成せねばならない．
Xllk10〕・X1｛k〕

X21k10〕＝X2㈹X3tklO〕・X3㈹DDD
榊一bl・nlメ1〕一fl・｛1〕。、｛、〕乳1，一f・1−11、ω払、．1〕執、1、叩
n20〕

fllω一hl｛1Wl〕一捌1〕ち・1㈱刀ω一f刀1〕一らメ1〕一f．m一ら。1〕イ。11〕
01ω

B7ω

1

十らl1〕十舳十n311〕十061〕十0411〕十D

11b12−P〕

一W

榊

D

D

1
1

1
一巧31P〕

榊D

1

1

一b231P〕

14・〕

日2−1〕

1

D
一f211P〕

一b311P〕

・。1・〕

榊

D
一ち21P〕

榊日1旧〕

n2｛P

−fl／p〕

bll｛P1他・bメp〕f2斐b2メp惚状ら外ち31P−b鶉［〕一ち1ρ

nl｛P〕Oメ・〕十ら1・〕州・1ρ〕十0・1・〕十0・一P〕十n・1・〕十XllklP〕

X21klP〕

X3｛k，P〕

図2三次元予測誤差フィルタ
Fig．2

Three

dimensionaI

prediction

一g4一

error刮tern・｛P）

強震記録の処理に関する一考察（その1）一木下

まず，Iを三行三列の単位行列として
F（n，O）＝一I，B（n，n＋1）＝一I，

n＝1，2，…

，p．

（12）を定義する．さらに，遅延作用素D：
DIX（n）＝X（n＿1），

1＝0，1，．．。

（13）を定義する．これらにより，（5），（6）式は各々（14），（15）式となる．

l

X（k11）＝一ΣF（1，n）・DnX（k）

（14）

n＝O
〜

1＋1

X（k一（1＋1）l1）＝＿ΣB（1，n）．D皿X（k）

（15）

n＝1

この両式と（10），（11）式から次の関係式が得られる．
〜

X（k

1＋1

l1＋1）＝＿Σ：F（1＋1，n）・D皿X（k）

n＝O
一一［・（1・…）・主［・（1，・）一・（1・・，1・・）・（1・・）1・ぴ・・（1・・，1・・）・刈・・（・）＝X（kl1）一F（1＋1，1＋1）・X（k一（1＋1）l1）

〜

（16）

1＋2

X（k一（1＋2）l1＋1）＝＿ΣB（1＋1，n）・DnX（k）

n＝1
一一・・［・（1・…）一之［・（1・・）一・（1・…）・・（1，・）1・㎜・・一1・（・）＝D［X（k一（1＋1）11）一B（1＋1，1）・X（k

l1）］

（17）

ここで，X（klO）＝X（k），x（k−110）＝DX（k）＝X（k−1）とする．
F（n，n），B（n，1）の（j，k）成分を各々（F（n，n））j，k＝fj，k（n），（B（n，1））』，止＝bj，止（n），

n＝1，2，…

，p．

j，k＝1，2，3．

（18）とすれば，（16），（17）式よりフィノレタΦは図2に示すようになる．このフィノレタヘの入カは三次元の独立なO−1間の一様乱数にΣ（p）の各主対角成分の平方根を乗じることにより近似的に得られる．すなわち入力列は下式で与えられる．
Xj（nlp）＝［（Σ（p））j，j］1／2・ηj（n），

E［ηj（m）］＝0，E［ηj（m）η。（n）］＝δj。・δ㎜，

j＝1，2，3；n＝1，2，…

（19）

j，k＝1，2，3；m，n＝1，2，…

ただし，（Σ（p））j，jはΣ（p）の（j，j）成分とし，δj，。等はクロネッカーのデノレタとする．

また，（16）式を1＝OからP−1まで辺々加えることにより下式を得る．
X（k）＝X（klO）
P〜＝ΣF（n，n）・X（k−nl〜n−1）十X（k

lP）

（20）

n＝1

（20）式は強震計記録の無相関なベクトル値確率変数列による有限フーリエ展開となる．
1．3

周波数領域における表現

一95一

国立防災科学技術セソター研究報告第17号

1977年3月

周波数領域におげる表現は｛F（n，n），B（n，1），n＝1，2，…，p｝から回帰行列｛F（p，n），n＝1，2，…，p｝を求めることに帰着する．これは（10），（11）式を逐次用いることによって求めら
れる．

そこで，

P
A（λ）＝I＿ΣF（p，n）．exp（＿iλn），

1川≦π

（21）

n＝1
とすれば，スペクトル密度行列は下式となる（E．Parzen，1968）．

1

〜

一

P（λ）＝一A■1（λ）・Σ（P）・［A1（2）］一1，

1川≦π

（22）

2π

ここで，λ＝2πf・△Tは基準化角周波数であり，AはAの複素供役とする．
1．4

数値計算における注意

実際の数値計算においては，共分散行列R（n）の代りに標本共分散行列
1N−n
N1＝1

C（n）＝＿Σx（（1＋n）．△T）．x（1．△T），

n＝O，1，…

，p．

（23）

を用い，（7），（8），及び（9）式を解いて必要なパラメータを求めることになる．ただし，｛x（n・△T），n＝1，2，…，N｝の加算平均は零ベクトルになっているものとする．

数値計算において決定せねばならないものとして，標本化時間△丁及び項数pがある．
この報告で用いている自己回帰法による推定スペクトノレ密度行列では項数不足のため近接するスペクトノレの山を分離出来ないことがある．近接する山の分離が真に意味を持つものならば，これには注意せねばならない．これを防ぐため，標本スペクトル密度行列I（λ）の特徴を把握しておく必要がある．I（λ）は標本フーリェ変換

N

Z（λ）＝Σ：x（n・△T）・exp（一iλn），

1λ1≦；π

（24）

n＝1
から

1
I（λ）＝一Z（λ）Z
2πN

（一λ），

1λ1≦π

（25）

として求まる．換言すれば，I（1）の特徴を目的に沿って表現出来るように△Tとpを決定せねばならない．

△TはI（λ）の形状から決定する．すなわち，I（λ）の形状から耐震工学的に意味を持つと兄なせる上限周波数f。（Hz）を決定する．しかる後に，△T（sec）を標本化定理からいくぶん余裕をもたせて（O・2〜O・3）・f；1程度に決定する．標本化の荒さによる折り返し効果はI（λ）と最終的に推定されるP（λ）との比較からその影響を判断する．ただし，前処理として低域通過フィルタが使えるような場合は別である．

バラメータの項数Pは情報量基準（AIC）に従って決定する（H．Akaike，1974）．
AIC＝一2・1n（最大尤度）十2・（バラメータ数）

一g6一

強震記録の処理に関する一考察（その1）一木下

情報量基準とはモデルの将来の値の分布を予測するための基準であり，真の分布の標本分布に関するエントロピーを最大にする（最大の情報量を得る）ためにモデノレのバラメータ白由

度を決定する方法である．この報告では二次までの統計量を扱っており，等価なガウス過程を想定出来る．したがって，実用的にはガウス分布の最大尤度から下式が求まる（R．H．
Jones，1974）．

AIC（n）＝N．ln（det（Σ（n）））十2M2n（26）

ここで，Mは次元の数である．項数pは（26）式を最小とするnとする．
2．一次元の場合

ここまでは強震計の三成分を同時に扱い，予測誤差フィノレタによる三次元地震動のシミュ

レーション法等の結果を得て来た．しかしながら，強震計記録の処理はその一成分のみに注目して行われているのが大部分である．また，実際の応用においても一次元での取扱いでほ
とんどの場合は問に合っているようである．そこで，以後は一次元の場合について詳論し，この方法の特徴や他の側面からの見方等を明確にする．

一次元の場合の記号はこれまでの太字記号から類推出来ると思われるので，添字記号は省賂し，なるべく見易いようにとした．
2．1

パラメータの特徴

一次元の場合，1で求めたパラメータは次のように三つの特徴が兄出される．まず，（3），（4）式に対応する前向き及び後向き予測を各々下式とする．
〈

lX（kl1）＝Σf（1，n）．X（k＿n）

（27）

n＝1
〈

1

■

X（k＿（1＋1）l1）＝Σb（1，n）．X（k＿n）

（28）

n＝1
各式の予測係数は（7），（8）式に対応して

1

R（n）＝Σf（1，m）．R（n＿m），

n＝1，2，…

，1．

（29）

m＝1

1

R（n一（1＋1））＝Σb（1，m）．R（n＿m），

n＝1，2，．．．，1．

（30）

m＝1
の解である．一次元の場合，共分散関数は偶関数であるから（30）式は下式となる．

1

R（n）＝Σb（1，1＿m＋1）．R（n＿m），
皿＝1

n＝1，2，…

，1．

（31）（29）式との比較により次の関係が得られる．
f（1，n）＝b（1，1−n＋1），

n＝1，2，…

，1．

（32）

故に，1．で定義したパラメータは一次元の場合には次のような関係を持っている．
f（1，1）＝b（1，1）≡r（1），1＝1，2，…

，p．

次に，前向き及び後向き予測誤差を（5），（6）式に対応して

一g7一

（33）

国立防災科学技術セソター研究報告第17号

1977年3月

X（kl1）＝X（k）一X（kl1）（34）

文（k一（1＋1）l1）＝X（k一（1＋1））一X（k一（1＋1）l1）（35）とすれば（34），（35）及び（32）式より
δ（1）＝E［（X（kl1））2］

1＝R（O）一Σf（1，n）・R（n）

n＝1
1

〜

＝R（O）一Σb（1，n）・R（1−n＋1）＝E［（X（k一（1＋1）l1））2］

（36）

n＝1

を得る．すなわち，前向き及び後向き予測誤差の自乗平均は等しい．さらに，この二つの予測誤差から次式が得られる．
E戊（k

l

1
l）・文（k一（1＋1）l1）］＝R（1＋1）＿Σ：f（1，n）．R（1−n＋1）

（37）

n＝1
最後に，（29），（36）及び（37）式を連立させることにより次の漸化式を得る．
1−1

R（1）＿Σ1f（1−1，n）・R（1＿n）
・（1）＝f（1，1）＝ n＋、

R（O）＿Σf（1−1，n）・R（n）

n＝1＿E［X（kI」1）X（k−111−1）］，1＿1，2，・，。．

（38）

E【（X（k＿111＿1））2］
f（1，n）＝f（1−1，n）一f（1，1）・f（1−1，1−n），

n＝1，2，…，1−1・（39）（39）式は（10）式に対応する．また，（38）式は（36）式により下式となる。

、（1）一一E［X（k11−1）聖k−111．1）］［E［（X（k11−1））2］・E［（X（k−l，1一・，2，。

（・・）l1−1））2］】1／2

ここで，X（k10）＝X（k），X（k−110）＝X（k−1）とする．Schwarzの不等式によればlr（1）1≦1となる．（40）式は偏相関係数の定義式である．

実際の数値計算においては，1＞PのときE［r2（1）1二N−1で評価される（G・P・Box

and

G・

皿．Jenkins，1970）から，N＝500〜1，OOO程度では偏相関係数の推定における有効桁はせい

ぜい三桁程度までである．また，（26）式は一次元の場合，偏相関係数の性質
6（n）＝R（O）・＝［圧［1−r2（1）］

（41）

1＝i

から，定数項を無視して下式となる．
AIC（n）＝NΣ1n【1−r2（1）］十2n
1＝1

2．2

（42）

仮想地盤

一次元の場合，バラメータ間には（33）の関係があり，予測誤差フィルタは図3とたる．
図3からも明らかであるが，（16），（17）式に対応する関係式を変形することにより次の二式が得られる．

一g8一

強震記録の処理に関する一考察（その1）一木下X（k11−1）＝X（kl1）十r（1）・X（k−1l1−1）
X（k一（1＋1）l1）＝D・［X（k−1l1−1）一r（1）・X（k（43）l1−1）］＝D・［一r（1）・X（kl1）十（1−r2（1））．X（k−111−1）］（44）

X（k411）X（klH）

X（klP）X（k−uユ→）X（k−H皿）

図3一次元予測誤差フィルタ
Fig．3

0ne

ror

dimensional

prediction

図4等価変換er−

Fig．4

Equiva1ent

transform

futer

U（工）

1−r（ユ）。

U（1−1）

X（）＝ヌ（klO）

U（O）
r（1）

十

1÷r（1）

1−r（1）

D（1）

D
DO

一r（1）

図5等価変換Fig．5

Equiv乱1ent

U（1）…

transform

D

lD（1）

U㈹

D

r（円

D（P→）

・1■

rQndomnumbers

1−r（P）

virtuαleαrth＿一r（P）
U（P）

ground・・α1・」αy…q・αk・
fαctorwαve

tiOnArtiicial

ea−rthquake

・1・

p
H p
岬一㈹岬一・（ΦX（klP）

図7人工波発生器
Fig．7

1・r（P）

accelera−

図6仮想地盤

generatOr

Fig．6

一99一

Yirtualground1ayers

国立防災科学技術セソター研究報告第17号

1977年3月（43），（44）式は図4のように表わされる．そこで，次の変換を行う．

X（k10）＝U（O），X（k−110）＝D（O）〜
1〜
l
X（kl1）＝1］1（1二r（n））・U（1），X（k一（1＋1）11）＝］］1（1−r（n））・D（1），
皿＝1

n＝1

1＝1，2，・．・，p．

（45）

この変換により（43），（44）式は各々次のようになる．
U（1−1）＝（1−r（I））・U（1）十r（1）・D（1−1）
D（1）＝D・［一r（1）・U（1）十（1＋r（1））・D（1−1）］

（46）
（47）

したがって，図4に対応して図5が得られる．また，図3は最終的に図6となる．

図6はU（p）を入射波とするp層の仮想な地盤における重複反射を示している．
1〜

1／■（1−r（1））はスケーノレファクタであるから，実際の入射波は分散6（p）の白色雑音X（k

l

p）

n＝1

である．また，r（n）は下降波の反射係数に対応し，L皿を換算層厚，V、を伝搬速度とすれば，標本化時間△Tとの間には△T＝2L皿／V皿の関係がこの仮想地盤の各層について言える．

さて，実際の地盤構成や波動の伝播機構の複雑さを考えると，地表での強震波を解析的に解くことは不可能に近いと思われる．そこで，この報告では図7に示すような強震波生成の確率的モデルを考える．このモデルでは，強震波の周波数特性を仮想地盤に，その大きさを
スケーノレ・ファクタにそれぞれ荷わせている．入力の一様乱数は観測強震波の位相特性が不規則であることから，これを積極的に利用出来る．

このモデノレは各震源一観測点毎に適用されるものであるが，個々のモデルにおいても仮想地盤の反射係数は観測波から推定されるもので，唯一的なものではない．したがって，各震源一観測点毎に多くの強震波が測定されれば，これらから確率変数としての反射係数を構成出来る．すなわち，仮想地盤を確率変数からなる系としてモデル化出来る．このとき，次の二点は工学的な応用に際して有用である．一つは全ての反射係数の絶対値が1よりも小さけ
れば，その系は安定であることである．他の一つは，この報告で考えている線形推定の範囲において，各反射係数は各々確率的に独立に処理が出来ることである（2．5参照）．最初は確率変数としての反射係数に適当な確率構造を入れてモデノレを行なうのも，一方法と思われる．
2．3

等価質点系

複素S平面において，振動のモードを表わす二次系の周波数応答関数は下式を含む．
（s2＋2hn・ωn・s＋ω宝）

1

（48）

z＝exp（s・△T）＝exp（伽）（49）

これに，変換を行えば（48）式は下式となる（E．I．Jury，1964）．

Z
H、（z）＝K，exp（一h皿ωn△T）．
（Z−Z。）（Zイ皿）
一ユOO一

（50）

強震記録の処理に関する一考察（その1）一木下

ここでz。＝exp（s皿．△T），

s皿＝一h、ω、十づω。（1−hζ）1／2

sm（ω、〉1−h芸

K皿＝

△T）

（51）

（52）

ω、〉1−h言とし，乏皿はZ皿の複素共役とする．したがって，二次系の伝達関数は下式を含む．
l

H皿（z）12＝H。（z）。H皿（z一工）＝K主． ZZ・（Z−Z皿）（Z−2；1）（Zイ。）（Z一乏；1）（53）

実際の地震動や構造物の振動はこの形の伝達関数の積の形をスペクトルに持つ時問波形として強震計に記録されている．したがって，強震計記録のパワースペクトノレから（53）式の形の因子の積を取り出すことが出来れぼ，（51）式によりそれぞれの振動モードの共振振動数と減衰定数が推定出来る．

一次元におけるスペクトル密度関数は（22）式に対応して下式となる．
P（λ）＝

6（P）
Q一ユ（z），
2π

1λ1≦π

（54）

ここで，特性多項式
Q（z）＝P
1一Σf（p，n）・z1

2
，

1z

l＝1

（55）

nヨ1
における予測係数｛f（P，n），n＝1，2，…，P｝は｛r（n），n＝1，2，…，P｝から（39）式を使うことにより求まる．そこで，Q（z）＝Oの根を｛z皿，n＝1，2，…，P｝とすれぼ，（54）式は下式となる．

6（P）
P（λ）＝一・1

2一血（・一州ド（一）P・δ（P）

2πII乞1
皿＝1

＝

。

・。

zp

II（z−zlユ）・（z−2≡1）
皿＝1

（56）

したがって，Q（z）＝Oの根にz皿及びその複素共役根乞皿が存在すれば（53）式の因子をP（λ）

から取り出し得る．すなわち，Q（z）を二次因子及び一次因子の積に分解したとき，二次因子から振動モードの共振振動数と減衰定数が決定出来る．

実際の数値計算において，高次代数方程式Q（z）＝0の解法はBairstow法（A．Ra1ston，
1965）が都合良い．また，この方法では共振振動数はかなり正確に推定されるが，減衰定数については共振振動数推定のような安定さと正確さは期待出来ない（木下繁夫・稲葉誠一，1976）．
2．4

擬似応答スペクトル

擬似応答スペクトノレの算出には人工地震波を用いる直接計算法と不規則振動論に基づく方法とが考えられる．ここでは時間領域における表現である予測誤差フィノレタを直接利用出来
る前者について考察する．
一101一

国立防災科学披術セソター研究報告第17号

1977年3月

まず，擬似相対変位応答スペクトノレの算出を考える．加速度入カ列に対する二次系の変位応答を｛W（n），n＝1，2，…，N｝とする．この二次系の周波数応答関数は形式的に（50）式と同じであるから，変位応答は下式より求まる．
W（z）＝一H。（z）．X（z）（57）

ここで，W（z），X（z）は各々｛W（n），n＝1，2，…，N｝，｛X（n），n＝1，2，…，N｝の片側z変換とする．系は静止状態から入力されたとすれぼ，（57）式は（50）式により下式となる．
W（n）＝［2cos（ω、、／1−h言△T）W（卜1）一exp（一h皿ω皿△T）W（n−2）

一K皿X（n−1）］・exp（一h皿ω、△T），

n＝2，3，…，N．

（58）

ただし，W（0）＝O，W（1）＝Oとする．（58）式は図8のようになる．図中のX（n）は図7における出カ加速度

である・

2cos（ωn（1−h書炉灯）

X（k）

擬似応答スペクトルは各固有振動数と減衰定数におげる｛W（n），n＝1，2，…，N｝との絶対最大値を求める

ことから算出される．さらに，図7における一様乱数

dn…exP←hnωn町）の組を換えることにより多くの擬似応答スペクトルを
Fig．8

W（k）

図8二次系
Secondorder

system求めることが出来る．したがって，それらに統計処理を行えば擬似応答スペクトノレの平均及び分散を推定出来る．また，この方法ではNと△Tの影響が含まれるので注意を要する．

微分演算は一階差分を取ることにより近似的に得られるので，これを利用すれぼ擬似速度
応答スペクトルも同様に算出出来る．
2．5

フーリェ展開確率変数列｛X（n），n＝1，…，N｝の線形結合である全ての確率変数からなるベクトノレ空間

をHとする．Hにおける二つの要素X（j），X（k）の内積を（X（j），X（k））＝E［X（j）・X（k）］とする．またX（j）のノルムを
l1X（j）い＝（X（j），X（j））1／2とし，X（j）とX（k）との距離をll

X（j）一X（k）llとする．

確率変数X（k）の無相関確率変数列によるフーリェ展開を求めるために，Hにおける閉部分空間M。＝｛X（k−n），n＝1，2，…，1｝を考える．M1と等しい無相関確率変数列｛V（n），n＝1，

2，…，1｝からなる部分空問はGram−Schmidt法により次式で得られる．
X（k−1）V（1）＝llX（k−1）ll
X（k−m l m−1）▽（m）＝ 〜，ll X（k−m lm−1）ll
一102一

m＝2，3，，1（59）

強震記録の処理に関する一考察（その1）一木下
ここで，文（k−mlm−1）＝X（k−m）一X（k−mlm−1）

（60）

m−1
文（k−mlm−1）＝Σ（X（k−m），V（n））・Y（n）

（61）

n＝1

とする．1を増して皿、＝｛V（n），n＝1，2，…，P｝を構成し，さらに部分空間｛叫1X（k）｝でこれを続けれぼ，
＾

X（k

l

P
p）＝Σ（X（k），V（n））・V（n）

（62）

n＝1

を得る．もちろん，X（klP）二X（k）一X（klP）

（63）

とすれば，

（文（klP），▽（・））＝O，・＝1，2，…，P．

（64）である．（64）式は文（k1p）がX（k）のM、の上への射影であり，距離11X（k）一X（klp川を最小にするという意味でM、によるX（k）の最良近似であることを示している（R・E・
Kalman，1960）．したがって，X（k）はP〜

X（k）＝Σ1（X（k），▽（n））．▽（n）十X（klp）

（65）

n＝1

とフーリェ展開される．
また，X（k）のM。．。＝｛▽（m），m＝1，2，…，n−1｝による近似を
A

X（k

l

n−1
n−1）＝Σ（X（k），V（m））V（m）

（66）

m＝1
とすれば，

X（kln−1）＝X（k）一X（kln−1）（67）としたとき，（文（kであるから，文（klln−1），V（m））＝0，

m＝・1，2，…

，n−1．

（68）

n−1）はM皿．、によるX（k）の最良近似である．したがって，（66）式より

（文（k

l

n−1），▽（m））＝0，

m＝n，n＋1，．・・，p．

（69）

であるから，（65）式の展開係数は（X（k），Y（n））＝（X（k）一X（kln一ユ），V（n））＝（X（k1仁1）lX（k一・I卜1））
ll

X（k−n

l

（70）

n−1）llとなる．故に，（38）式及び内積とノノレムの定義から
r（n）＝

〜

（X（k），V（n））（71）

llX（k−n1n−1）l12となるから，（65）式はP〜〜

X（k）＝Σr（n）・X（k−nln−1）十X（klp）
I1＝1

−103一

（72）

国立防災科学技術セソター研究報告第・・号・・。。年。月となる・すなわち・｛・（・）・・一・・・・・・…／は無相関確率変数列／文（・一、一ト。），、、1，、，．．．

p｝の展開係数となる．したがって，（「（n）X（k11卜・）・・（・）・文（・一・1・一））一・，。・。

（。。）

三あるから・各展開項は相関がな／確率的に独立に披／ことが出来る利点がある．また，二（k−n）I（卜1）・n＝1・a・・…1が無相関なべl／ル値確率変数111なるll／同様に言え

3．数値計算例
ここまで述べてきた方法の有効性を確認するため，実地震波について数値実馳行つた．
実土芒震波としては・S・・・・…i・・地震（・・・…．・・）の・・1・。。・、。、。、、。の記録を用いたが・選定の理由は特別にはな／・・数値計算は二次元の場合について行つた．すなわち，
苓1（k）をN・S成分・X・（・）を・一・成分として，・一…，・・一・・m、とした．（、6）式の規準ではn＝1・2パ・・・…の範囲において・一・・が求まつた．ただし，。、、である．
図9に二次元の予測誤差フィ／レタを，表・に推定された各パラメータを各々示す．実際の計算で

X1（k）

は仮数部37ビットの実数計算であるが，表1は

X、（k）

D

小数点以下・桁目を四捨五入したものである．ス柵岬一f！1）一・・1〕・〕一・。／l）D舳

ペクトル密度行列の推定に対して，この四捨五入

十はほとんど影響しない．

十

、

D

、
D

図10■1〜図10−4までに原記録の標本スペク／ l l ， ．
ノレ密度行列と推定スペク／ル醐了列を示す．横 I I l 1
軸は全て

（21

・T）である・密度行列は・…i・・一棚岬畏。l1一・、、舳，．妙，宅、、同行列であるから，対非角成分はその実数部と虚数部とで記述される．図・・は表・のバラメータを

十

十

十

・

用いた図9のフィルタによる合成二次元波を示している．図は最初の4秒間を示している．計算では△T＝20・・として・・・…サンプル準生させ

X1（kIPリヌ。（klPl

図9二次元予測誤差フィルタ

・i・・・・…i。。、、i。、、1、、、。i、．

ている・後半の…サンプルの標本スペクトノレ密tione「「o「釦t・・芒一列を図10−1〜図10−4に柑る推定スペク／／／密度行列と対比させて，図。、．。一図12．、
1」不す．

図9の入カ列は・η（k）をO−1問の一様乱数として，
1（・）一［1、；（k）11／2・・（・π1（…））（74）

1（・）一［、、；（、）丁／2…（・π・1（…））
（75）

一104一

強震記録の処理に関する一考察（その1）一木下

表1二次元予測誤差フィルタのパラメータ

Table1∫。。（仇）

〃

1
2
3
4
5
6
7
8
9

Parameters

of

the

∫。。（犯）∫。、（椛）

twodimensionaI

prediction

61。（〃）

6．1（犯）∫。。（物）

error

i1ter
b。。（刎）

b31（犯）

O．740

O．124

一〇．008

O．538

O．749

一〇．093

0．032

O．529

一〇．709

一0．100一〇．010一〇．633

一〇．7150．021

一〇．020一〇．626

O．357

O．043

O．008

O．405

O．359

一〇．O03

O．O06

O．401

一0．380

0．183

O．035

一〇．264＿O，387

O．120

O．072

一〇．265

一0．112

一〇．002

O．015

O．069

一〇．116O．044

一〇．007O．068

一〇．066

一〇．326

一〇．049一〇．019一0．069一〇．118＿O．127

一〇．O060．038

一〇．098

一〇．179

O．O09

O．031

一0．108

O，036

一0．072

一〇．149一〇．102

O．062

O．068

一〇、167O．181

一〇．044O．070

一〇．085O．148

一〇．015

一〇．182

一〇．089一〇．023O．071

一0．18110

0．016

一〇．114

一0．O09一〇．026O．019

一〇．026

一〇．043

一0．01711

一〇．023

O．117

0．011

一〇．077

一〇．027

0．033

O．051

一〇．080

12

0．163

O．303

一〇．017

一〇．171

O．183

一0．079

O．131＿O．181

5

10

Σ（・・）一［3；l111111

2

3
共10

−2

。18 gd・SeC

2−2

gd・SeC

lO一
10

oj

O．1

0

0

5

図10．1

Hz

0

15

スペクトル密度関数（記録波，W−S

図10．2

tion

15

スペクトル密度関数（記録波，σ一D

成分）

成分）

Fig．10．1

Hz

Power

spectral（recorded！V−Sdensity
comPonent）

func−

Fig．10．2 Power spectra1density func−tion（recordedσ一刀comPonent）

一105一

国立防災科学技術セソター研究報・告

10

2

10

−2

gql・SeC

・18

・18

一5

一5

一10

一10

5
図10．3
Fig．

10．3

l

1O

Hz

1977年3月

gQlく・eδ2

1
515

コ スペクトル（記録波）
Co−spectrum（recorded）

400

第17号

図10．4

1Hz
10

クウァッド・スペクトル（記録波）

Fig． 10．4 Quard−sPectrtm（recorded）

N−S

gα1

SeC

一400
U−D

200

3

1

gQl

1
15

O

SeC

一200
図11二次元合成波（0−4秒）

Fig． 11 A sample of two dimensiona1artiicial acceleration wave（O−4sec）

一106一

強震記録の処理に関する一考察（その1）一木下

2
−2
gdSeC

3
×10

。1ゴ

2 −2
gd・sec

10
10

O」

0．1

5

10

Hz

15

0

図12．1スペクトル密度関数（合成波，N−S

図12．2

Power

spectraltion（synthesized

W−S

densityfunc−

comPonent）

Hz

15

スベクトル密度関数（合成波，r刀Fig．12．2 Power spectral density func−tion（synthesizeφひ一D comPonent）・

10

10

gdろS・62。1♂

2−2
9Ql・SeC

刈03

一5

一5

一10

一10

0

10

成分）

成分）
Fig．12．1

5

l

5

l

10

Hz

図12．3 コ・スペクトル（合成波）
酬g．12．3 Co−spectrum（synthesized）

Hz
l

1O

15

1尋

l

図12．4 クウァッド・スベクトル（合成波）
Fig．12．4 Quad−sPectrum（synthesized）

一107一

国立防災科学技術セソター研究報告第17号

なる二次元の独立な平均零，分散1の白色ガウス列（B．Go1d

1977年3月

andC．M．Rader，1969）を用いて下式から構成した．
麦。（klp）＝［（Σ（p））。。］1／2・ζ（k）（76）

気（・1・）一1（Σ（・））rllll；llゴ・1（・）・1姜1；；；：：為（・1・）（・・）

また，表1からも知れるように各バラメータ問にはf・，。（n）竺b。，。（n），f。，。（・）！b。，。（n），

n＝1，2，…，P．

（78）の関係が見出される．すなわち，各成分におけるバラメータ間には（33）式の関係が近似的に成立する．これは二次元予測誤差フィルタにおいても仮想地盤の概念が拡張出来ることを意味している．

4．結論

強震計記録主要動部を定常確率過程とみなして，これをバラメータ化した．これらから，

周波数領域における表現と三次元までの模擬地震波合成法を導いた．この方法は従来の三角関数系によるフーリェ展開と異なり，無相関確率変数列によるフーリェ展開を基礎にしている．その利点は三角関数系と比較してはるかに少い展開項数で地震動の特徴を表現出来ることであり，多数の強震記録をファイルしておくのに適した方法と考える．

あとがき
この報告では，方法の確立に重点を置いた．そのため簡単な数値計算例を一つ述べたにすぎなかったが，次報告では耐震工学に使えるような型に適用例をまとめる予定である．また，

パラメータを時間的に変動させた場合の非定常的な取扱いについても考察する予定である．
謝

辞本研究の内容について，耐震実験室の諾兄から種々の御意見を頂いた．深く感謝致します．

また，一次元予測誤差フィルタの他分野への応用を扱った次の二文献は多くのことを参考にさせて頂いた．各著者の方々に感謝の意を表します。板倉文忠・斉藤収三（1969），D．C．Ri1ey
a−nd

J．P．Burg（1972）．

1）Akaike，H．（1974）：

19，

A

new

look

at

参

考

the

statistical

献
model

identiication，1亙亙亙仰α？z8．五．α，6，716−723．

2）Box，G・P・and G・M・Jenkins（1970）：
Holden−Day，65．
3）

文

Gold，B−and

C・ユV［・Rader（1969）：

Time

DigitalseriesanaIysis，forecastingprocessing

of

and

control，

signa1s，McGraw−Hill，144−146．

4）板倉文忠・斎藤収三（1969）：偏自已相関係数による音声分析合成系，日本音響学会講演論文集，
199．
5）Jones，R．H．（1974）：

Identi丘cation

anda・utoregressive

一108一

spectrumestimation，1亙亙亙〃α伽．

強震記録の処理に関する一考察（その1）一木下■4．0．，

19，6，894−898．

6）Jury，E．I．（1964）：Theory

andapplication

newapproachoftheZ−transform

linearmtering

method，John−Wiley＆Sons，298．

7）Kalman．R．E．（1960）：

A

to

and

prediction

problems

λ．8．〃．亙．，∫．3αs｛o亙刎伽θγ伽σ，82，35−45．

8）木下繁夫・稲葉誠一（1976）：電力用遮断器の振動実験，防災科学技術研究資料，第24号．
9）Parzen，E．（1968）：
ofMultiplestatistics，Stanford

10）Ralston，A．（1965）：

11）Riley，D．C，and

A

time

seriesmode1ing，Zθ6んっo｛ωZ㈹ρo材No．12，Department

University．

irst

coursein

J．P．Burg（1972）：

numericaI

Time

and

analysis，McGraw−Hin，372−378．

sPaceadaptive

deconvo1ution

刑ters，〃θ〃棚0グ伽42棚舳θ亡伽90∫伽S㏄｛物0グ肋ρ10㈹肋犯Gθ0ρ肋8械8．
12）Whittle，P．（1963）：On
canonical

factorization

the丘tting
of

a

of

spectralmultivariate
density

autoregressions，and

theaPProximate

matrix，別o伽θか肌α，50，129−134．
（1976年11月17日原稿受理）

一109一