アプリケーションに戻る

PRISM の使い方

はじめに

PRISM は，加速度の時刻歴から，時刻歴波形（加速度・速度・変位），フーリエスペクトル，地震応答スペクトルおよび波形の特徴（最大値，卓越振動数，計測震度）を求めるアプリケーションである。強震観測網 K-NET・KiK-net^[0-1] と気象庁^[0-2] の記録形式，および数値のみを並べた CSV を読み込む。計算はすべてブラウザ内で行い，データは送信しない。

• [0-1]防災科学技術研究所：強震観測網（K-NET，KiK-net）
• [0-2]気象庁：強震波形データ

1. データの読込

画面上部の「1. データの読込」にファイルを与える。複数のファイルを続けて読み込むと，一覧に追加される。読み込める形式は次のとおりである。

• K-NET・KiK-net の記録（拡張子 .NS，.EW，.UD など）：ヘッダの Scale Factor を乗じて cm/s² に換算し，NS・EW・UD を別々の波形として登録する。KiK-net の地中記録（Dir. が 1〜3）は除く。
• 気象庁の強震波形（.csv）：観測点コード，緯度経度，サンプリング周波数を読み取り，単位が m/s² の場合は 100 倍する。
• 常時微動計の CSV：先頭行に EPSON_X，EPSON_Y，EPSON_Z の列名があるファイル（コンマ区切りまたはタブ区切り）を見分け，この3列を加速度（単位 mG，1 mG = 0.980665 cm/s²）として読み込む。値のない行，温度（EPSON_T）および電池電圧（VBat）の列は読み飛ばす。読み込むときに確認の画面を表示し，各軸の成分（NS・EW・UD），向き（反転すると値に −1 を乗じる）およびサンプリング周波数を指定する。センサの X 軸は左向きが正である。サンプリング周波数は，UnixTime（秒）が同じ行の数を数えて求める（最初と最後の端数の秒は除く）。求められない場合は「数値CSVの既定の周波数」の値を初期値とする。名前はファイル名とし，前回の成分と向きの指定を次回の初期値とする。
• 数値のみの CSV・テキスト：1列を1波形として扱う。ヘッダ行は入れない。
• 上記を含む ZIP（.zip，.tar.gz，.tgz）：展開して読み込む。

読み込んだ波形は一覧に1行ずつ表示される。一覧の各列は次のとおりである。

• 名前・成分：強震記録は観測点コードと成分（NS・EW・UD）を入れる。数値のみの CSV はファイル名（複数列のときは ch 番号を付ける）を入れ，名前の末尾が「_NS」「 EW」などの場合は成分として読み取る。名前と成分は直せる。名前が同じ波形は1つの観測点として名前の欄をまとめ，成分を NS・EW・UD の順に並べる。オービット，水平合成および計測震度は，この名前と成分で組を作成する。
• 周波数・倍率・単位：波形ごとに直せる。サンプリング周波数は，強震記録ではヘッダの値を，数値のみの CSV では「数値CSVの既定の周波数」の値を入れる。単位では，読み込んだデータの量と単位を選択する（加速度：cm/s²・m/s²・G・mG，速度：cm/s・m/s）。計算は cm/s²（速度は cm/s）に換算して行う。倍率は単位とは別の補正の係数であり，読み込んだ値に乗じる。
• 点数・長さ：標本の数と記録の長さ（秒）を示す。
• 時刻歴・フーリエ・応答・特徴・減衰：波形ごとに，求めるものを選択する。読み込んだ時点ではすべてにチェックが入る。見出しの欄で列をまとめて，左端の欄で行をまとめて切り替える。右端の「×」で一覧から消す。
• CSV出力：一覧の下のボタンで，計算済みの結果を種類ごとのファイル（time_history.csv，fourier_spectrum.csv または transfer_function.csv，response_spectrum.csv，features.csv，damping.csv，オービットがあれば orbit.csv）に書き出す。

サンプリング周波数の異なる波形も重ねて描けるが，伝達関数，オービットおよび水平合成は，組になる波形の周波数がそろっている必要がある。

1.1 速度のデータの扱い

一覧の単位で cm/s または m/s を選択すると，読み込んだ値を速度として扱う。平均値の差し引き，倍率と単位の換算（cm/s に統一する）およびフィルタは，速度に対して行う。加速度は，フィルタ後の速度の各時刻における傾きとして，式(1.1)の中央差分で求める。前後1点ずつの差をとるため，その時刻のちょうど中央の傾きとなる。最初と最後の点では，隣の1点との差を Δt で割る。変位は，速度を台形則で積分したのち，同じフィルタを適用し直して求める（図1.1）。

(1.1)

ここで，

• a_i：時刻 i の加速度（cm/s²）
• v_i：時刻 i のフィルタ後の速度（cm/s）
• Δt：時間刻み（s）

中央差分は，高い振動数の成分ほど振幅を小さく見積もる。振動数 f の成分の振幅は sin(2πfΔt)/(2πfΔt) 倍となり，たとえば f_s = 100 Hz のとき，10 Hz の成分で約 6 % 小さくなる。建物の応答で重要な 10 Hz 以下では，実用上の影響は小さい。

各計算での扱いは次のとおりである。

• 時刻歴波形：速度はそのまま，加速度は式(1.1)，変位は台形則で積分して3段に描画する。
• フーリエスペクトル：速度のまま求める。縦軸の単位は cm/s である。
• 地震応答スペクトル：式(1.1)で求めた加速度を入力とする。
• 伝達関数：速度どうしの組だけを計算する。加速度と速度の組は省く。
• 波形の特徴：最大値は時刻歴波形と同じ値とし，卓越振動数は速度のフーリエ振幅から求める。計測震度は，フィルタを通さない速度（平均値を差し引いたもの）から式(1.1)で求めた加速度を用いる。

表1.1　読み込めないときのメッセージと対処

メッセージ	原因と対処
ヘッダを読み取れません	K-NET・KiK-net の記録として読もうとしたが，ヘッダの Station Code，Scale Factor または Sampling Freq(Hz) が見つからない。ヘッダを削除した記録や，表計算ソフトで編集した記録では起こる。取得した元のファイルを用いる。
成分（Dir.）を判別できません	K-NET・KiK-net のヘッダの Dir. が NS・EW・UD（KiK-net の地中記録では 1〜6）のいずれでもない。拡張子だけを変えたファイルでは起こる。
サンプリング周波数を読み取れません	気象庁 CSV として読もうとしたが，SAMPLING RATE の行が見つからない。
NS・EW・UDの見出し行が見つかりません	気象庁 CSV の列見出しの行が見つからない。見出しを削除したファイルでは起こる。
ZIPの形式を読み取れません	ZIP として読み取れない。ファイルが壊れているか，分割された ZIP である。展開してから読み込む。
数値の行が足りません	数値だけの CSV として読もうとしたが，数値として読める行が 64 行に満たない。ヘッダ行や時刻の列が含まれていないか確認する。
データが空です	ファイルに読み取れる行がない。
データが短すぎます	フーリエスペクトルを求めるには 64 点以上が必要である。区間の長さを短くするか，長い記録を用いる。

2. 求めるもの

一覧の「時刻歴」「フーリエ」「応答」「特徴」「減衰」の欄でチェックした波形と，「2. 計算」で作った伝達関数の組について，「計算する」を押すと次の6つを一度に求める。結果は右の欄のタブで切り替えて表示する。各計算の条件は，「2. 計算」の折りたたみで指定する。

• 時刻歴波形：加速度，速度および変位を縦に3段で描画する。速度は加速度を，変位は速度をそれぞれ台形則で積分し，積分のたびに同じフィルタを適用し直す。速度のデータの場合は1.1節による。観測点コードが同じ NS・EW の組があるときは，NS–EW のオービット（水平面内の軌跡）も描画する。「時間の範囲」に開始と終了の時刻（秒）を入れると，その範囲だけを描画する。空欄のときは全長を描画する。
• フーリエスペクトル：波形を区間に分割し，窓関数を適用して求めた各区間のスペクトルを平均する^[2-1]。求めるものは，フーリエ振幅スペクトル，フーリエ位相スペクトル，パワースペクトルおよびパワースペクトル密度から選択する。
• 伝達関数：「2. 計算」の「伝達関数」で，入力（基準）と出力の組を作成する。組は「＋組を追加」で1組ずつ作成するか，「まとめて作る」で基準の観測点を選んで「他の観測点の同じ成分と組にする」を押し，他の観測点の同じ成分（NS→NS など）との組をまとめて作る。組ごとに，区間ごとのクロススペクトルの平均から伝達関数を求め，振幅比，位相差（−180〜180度）およびコヒーレンスを描画する。区間・窓・平滑化はフーリエの条件を用いる。組になる2つの波形は，サンプリング周波数と量（加速度か速度か）がそろっている必要があり，そろわない組は省く。
• 地震応答スペクトル：固有周期の異なる多数の1自由度系に同じ波形を入力し，最大応答を周期ごとに並べる。変位 S_d，擬似速度 _pS_v または擬似加速度 _pS_a を求める。解き方と入力項目は4章に示す。
• 波形の特徴：最大加速度・最大速度・最大変位とその発生時刻，卓越振動数（周期）を一覧にする。最大値は時刻歴波形と同じ処理（①のフィルタ，積分のたびに適用し直し）で求める。卓越振動数は，②〜⑤の条件で求めたフーリエ振幅が最大となる振動数である。同じ観測点の NS・EW・UD がそろう場合は，気象庁の方法^[2-2]で計測震度と震度階級も求める。計測震度には，平均値を差し引いた3成分を用い，①のフィルタは通さない。
• 減衰：卓越振動数付近の減衰定数 h を，RD法，対数減衰率およびハーフパワー法から選択した方法で求める。対象の振動数 f₀ は成分（NS・EW・UD）ごとに指定する。空欄の成分と，成分を指定していない波形では，フーリエ振幅が最大となる振動数とする。求め方は下に示す。

フーリエ振幅は式(2.1)による。ここで w_n は窓関数，L は区間長，M は区間の数であり，窓の総和で正規化しているため，振幅の単位は元の波形と同じである。ただし，直流（k = 0）とナイキスト振動数（k = L/2）では，係数 2 を 1 とする。

(2.1)

ここで，

• A(fk)：振動数 f_k におけるフーリエ振幅。単位は元の波形と同じ
• f_k：k 番目の振動数。f_k = kf_s/L（Hz）
• L：1区間の点数。512〜16384 または全長に収まる2のべき乗
• M：平均に使った区間の数
• w_n：窓関数の値。区間の端を滑らかに 0 に近づける
• x_m,n：第 m 区間の n 番目の標本値
• i：虚数単位

フーリエ位相スペクトルは，各振動数成分の偏角 ∠X(f_k) であり，−180° 〜 180° で表す。「アンラップ」を選択すると，隣の振動数との差が 180° を超えないよう 360° を加減して連続させる。区間ごとの位相を平均しても意味がないため，位相を選んだときは区間の長さを全長で1回とし，区間の選別と平滑化は行わない。

パワースペクトルは，各振動数成分の平均パワーであり，振幅 A(f_k) から A²/2（直流とナイキスト振動数では A²）として求め，区間について平均する。単位は元の波形の単位の2乗である。窓を適用しない場合，すべての振動数について足し合わせると波形の二乗平均に等しい。パワースペクトル密度ではないため，値は区間の長さによって変わる。

パワースペクトル密度（PSD）は，1 Hz あたりのパワーであり，式(2.2)により求めて区間について平均する（Welch 法^[2-1]）。窓の2乗和で割ることで，窓による実効的な帯域幅の違いを補正している。単位は元の波形の単位の2乗を Hz で割ったものであり，振動数について積分すると波形の二乗平均に等しい（窓を適用した場合は，窓の重みを付けた二乗平均となる）。区間の長さやサンプリング周波数によらないため，条件の異なる記録どうしを比べるのに適する。

(2.2)

ここで，

• G(f_k)：片側のパワースペクトル密度。直流とナイキスト振動数では係数 2 を 1 とする
• X(f_k)：窓を適用した区間の離散フーリエ変換
• f_s：サンプリング周波数（Hz）

伝達関数は，組の入力 x と出力 y から，クロススペクトルの平均から式(2.3)^[2-3]の H₁ 推定量として求める。コヒーレンスは式(2.4)による。位相差は H の偏角であり，−180° 〜 180° で表す。

(2.3)

ここで，

• H(f)：伝達関数。入力に対する出力の応答倍率
• S_xy：入力 x と出力 y のクロススペクトル（区間平均）
• S_xx：入力 x のパワースペクトル（区間平均）

(2.4)

ここで，

• γ²(f)：コヒーレンス。0〜1 の値で，1 に近いほど2つの波形の関係が線形で安定していることを示す
• S_yy：出力 y のパワースペクトル（区間平均）

減衰は，対象の振動数 f₀ を中心に ±帯域（既定 10%）のバンドパスフィルタ（Butterworth 4次，ゼロ位相）を適用した波形 y から求める。ハーフパワー法だけは，フーリエの条件（②〜⑤）で求めたフーリエ振幅から求める。

RD法（ランダムデクリメント法）^[2-4]は，常時微動のような不規則な振動から自由振動に相当する波形を取り出す方法である。y が上向きにトリガの値（標準偏差の倍数，既定 1.0 倍）を横切った時刻 t_j から，決まった長さ（既定 10 周期分）を切り出して式(2.5)で平均する。不規則な成分は平均すると打ち消し合い，初期変位を与えた自由振動の波形 D(τ) が残る。平均した回数 N が多いほどばらつきが小さくなる。得られた D(τ) には，次の対数減衰率の求め方を適用する。

(2.5)

ここで，

• D(τ)：平均した波形（自由振動に相当する波形）
• τ：切り出した時刻からの経過時間（s）
• N：平均した回数
• t_j：j 番目にトリガの値を上向きに横切った時刻（s）
• y：バンドパス後の波形

対数減衰率は，自由振動の隣り合う山の振幅の比の対数である。指定した区間（「区間」の始め〜終わり，空欄は全長）のうち最大の山から順に山を拾い，最初の山の 5% を下回るまでの山の振幅 a_k の対数を時刻 t_k に対して直線で近似し，式(2.6)の傾き β から1周期あたりの対数減衰率 δ を求める。減衰定数は式(2.7)による。自由振動の部分を区間に指定して用いる。

(2.6)

ここで，

• a_k：k 番目の山の振幅
• t_k：k 番目の山の時刻（s）
• c：直線の切片
• β：直線の傾き（1/s）
• δ：1周期あたりの対数減衰率
• f₀：対象の振動数（Hz）

(2.7)

ここで，

• h：減衰定数

ハーフパワー法は，フーリエ振幅の山の高さ A₀ の 1/√2 となる振動数 f₁，f₂（f₁ < f₀ < f₂）を読み取り，式(2.8)で求める。直線補間で f₁，f₂ を求める。振動数の刻みが山の幅に比べて粗いと山が低く広く見え，細かすぎると区間の数が減ってばらつきが大きくなるため，結果は区間の長さによって変わる。

(2.8)

ここで，

• f₀：フーリエ振幅が最大となる振動数（Hz）
• f₁, f₂：フーリエ振幅が A₀/√2 となる低い側と高い側の振動数（Hz）
• A₀：山の高さ（フーリエ振幅の最大値）

「減衰」の折りたたみで指定する項目は次のとおりである。

• トリガ（RD法）：標準偏差の何倍の値を横切った時刻から切り出すかを指定する。小さくすると平均の回数が増えてばらつきが減るが，振幅の小さい揺れの減衰となる。1.0 倍前後がよく用いられる。
• 切り出し（RD法）：5〜10周期分を目安とする。短いと減衰が進まず傾きが決まりにくく，長いと末尾で打ち消し残りの雑音の割合が大きくなる。
• 帯域：RD法と対数減衰率ではバンドパスフィルタの幅として用いる。ハーフパワー法では山を探す範囲にだけ用い，f₁，f₂ の読み取りには影響しない。

• [2-1]Welch, P. D.: The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms, IEEE Transactions on Audio and Electroacoustics, Vol. AU-15, No. 2, pp. 70-73, 1967, https://doi.org/10.1109/TAU.1967.1161901
• [2-2]気象庁：計測震度の算出方法（平成8年気象庁告示第4号）
• [2-3]Bendat, J. S. and Piersol, A. G.: Random Data: Analysis and Measurement Procedures, 4th ed., Wiley, 2010, https://doi.org/10.1002/9781118032428
• [2-4]Cole, H. A.: On-line failure detection and damping measurement of aerospace structures by random decrement signatures, NASA CR-2205, 1973

3. 計算条件

「2. 計算」では，波形に施す順に，共通の条件（①平均値の差し引きとフィルタ）と，フーリエの条件（②区間の長さと重なり，③区間の選別，④窓関数，⑤平滑化）を指定する。フーリエスペクトルでは①〜⑤のすべてを，時刻歴波形と地震応答スペクトルでは①だけを用いる。波形の特徴では，最大値に①を，卓越振動数に①〜⑤を用いる。

3.1 波形の下ごしらえ（①）

読み込んだ値から，「平均値を差し引く」を選んだ場合は全長の平均値を差し引き，倍率を乗じたのちにフィルタを適用する。フィルタは，種類（ローパス・ハイパス・バンドパス），特性（Butterworth・Chebyshev I・Bessel），次数（2・4・6・8）および位相（ゼロ位相・因果）で指定する。既定はハイパス 0.05 Hz，Butterworth 4次，ゼロ位相である。

種類　ローパスは遮断周波数 f_c より高い成分を，ハイパスは f_c より低い成分を除去する。バンドパスは下限 f₁ と上限 f₂ の間の成分だけを残す。遮断周波数がナイキスト周波数 f_s/2 以上のときは，0.8f_s/2 に下げて計算し，その旨を表示する。3つの種類の働きを図3.1に示す。

特性　各特性は，遮断角振動数を 1 rad/s とするアナログのローパス原型で定義する。Butterworth の振幅特性は式(3.1)で表され，通過域が最も平坦である^[3-1]。

(3.1)

Chebyshev I の振幅特性は式(3.2)，式(3.3)で表される。通過域に高さ R_p のリップル（波打ち）を許すかわりに，同じ次数でも遮断が急になる。R_p は「リップル」欄で指定し，既定は 0.5 dB である。偶数次では直流の利得が 1 を下回るため，通過域の最大値が 1 となるよう全体を補正している。

(3.2) (3.3)

ここで，

• H(f)：フィルタの伝達関数
• f_c：遮断周波数（Hz）
• n：次数
• T_n：n 次の Chebyshev 多項式
• ε：リップルの大きさを表す係数
• R_p：通過域のリップルの高さ（dB）

Bessel は，群遅延（振動数ごとの時間の遅れ）が通過域でほぼ一定となるように定めた特性であり，波形の形を崩しにくいが，遮断はゆるい。本アプリケーションでは，遮断周波数で振幅比が 1/√2（−3 dB）となるよう規格化している。3つの特性の比較を図3.2に示す。

次数　次数 n を上げるほど遮断域の傾きは急になり，ローパス・ハイパスでは，遮断域で振動数が2倍（または1/2）になるごとに約 6n dB 減衰する（図3.3）。一方，次数が高いほど遮断周波数付近での位相の変化が大きく，因果フィルタでは波形のゆがみや遅れが増える。バンドパスは原型を帯域通過形に変換するため極の数が 2n となり，下限側と上限側がそれぞれ n 次の傾きをもつ。通常は 4次で十分である。

デジタル化　アナログ原型を式(3.5)，式(3.6)で目的の種類に変換し，双一次変換（式(3.7)）でデジタルフィルタとする^[3-2]。双一次変換では振動数軸がゆがむため，遮断周波数をあらかじめ式(3.4)で補正（プリワープ）しておく（図3.4）。式(3.5)の左はローパス，右はハイパス，式(3.6)はバンドパスへの変換である。得られたフィルタは2次ずつの部分に分けて順に通し，高次でも計算の誤差が大きくならないようにしている。

(3.4) (3.5) (3.6) (3.7)

ここで，

• Ω_c：プリワープした遮断角振動数（rad/s）
• f_s：サンプリング周波数（Hz）
• Ω₁，Ω₂：バンドパスの下限・上限をプリワープした角振動数（rad/s）
• s：ラプラス変換の複素変数
• z⁻¹：1標本分の時間の遅れを表す演算子

位相　因果は，波形を前向きに1回だけ通す方法である。実時間の処理と同じく過去の値だけを使用するが，位相が遅れてピークの時刻が後ろにずれる。ゼロ位相は，前向きに通したのちに時間を反転してもう1回通す方法であり，位相の遅れが打ち消されてピークの時刻が移動しない（図3.5）。その代わり振幅特性は2乗となり，遮断周波数における振幅比は 1/√2 ではなく 1/2（−6 dB）に，遮断域の傾きは2倍になる。また，応答が時間の前方にも広がるため，入力より前に小さな揺れが現れる。地震記録の処理では，位相を保つためにゼロ位相を用いるのが一般的である。

記録の始めと終わりでフィルタの過渡的な応答が生じないよう，波形の前後に 0 を補ってからフィルタをかけ，後で取り除く。補う長さは，ハイパス・バンドパスでは 1.5n/f₁ 秒（最大 120 秒）である。記録の両端で値が 0 から大きく離れている場合は，「平均値を差し引く」を併用するとよい。

選び方の目安　地震記録を積分して速度・変位を求める場合は，基線のずれによるドリフトを抑制するためにハイパス（0.05〜0.1 Hz 程度）を適用する。遮断周波数を高くするほどドリフトは減るが，長周期の成分も失われるため，対象とする周期の範囲より十分に低い振動数とする。高振動数のノイズを除く場合はローパスを，特定の帯域だけを見る場合はバンドパスを用いる。

3.2 区間に分割する（②）・使用する区間を選択する（③）

区間の長さ L は，全長で1回とするか，512〜16384 点から選択する。全長で1回とした場合は，全長に収まる最大の2のべき乗の点数を先頭から用いる。振動数の刻み（分解能）は式(3.8)で決まり，区間が長いほど細かくなる。たとえば f_s = 100 Hz，L = 4096 点では Δf ≈ 0.024 Hz であり，1区間の長さは約 41 秒である。

(3.8)

ここで，

• Δf：振動数の刻み（Hz）
• L：1区間の点数

一方，区間を短くして数を増やすと，平均によりスペクトルのばらつきが小さくなる（Welch 法^[2-1]）。分解能とばらつきは両立しないため，見たい振動数の間隔と記録の長さから区間の長さを決定する。重なりは隣の区間と共有する割合（0%・25%・50%・75%，または任意の 0〜95%）であり，窓によって区間の端の重みが小さくなる分を補う。窓を用いるときは50%とすることが多い。

③を選択すると，各区間の二乗平均平方根を求め，全区間の中央値に対する比が指定の範囲を外れる区間を除く（既定は 3 倍超と 0.1 倍未満）。常時微動から，人の通行などによる一時的な乱れを除くために用いる。すべての区間が除かれる場合は，選別せずに全区間を用いる。伝達関数では，入力の波形で選別する。

3.3 窓関数（④）

有限の区間を切り出してフーリエ変換すると，区間の両端で波形が不連続となり，本来ない振動数に振幅が漏れ出す（漏れ）。窓関数は区間の端をなめらかに 0 に近づけて，この漏れを抑制する^[3-3]。Hann を式(3.9)，Hamming を式(3.10)，Blackman を式(3.11)に示す。「なし」は全区間で w_i = 1 とする矩形窓である。各窓の形を図3.6に示す。

(3.9) (3.10) (3.11)

ここで，

• w_i：窓関数の値
• i：区間内の位置（0〜L−1）

窓関数の周波数特性を図3.7に，主な指標を表3.1に示す。中央の山（メインローブ）の幅が広いほど，近接した2つの振動数を分けにくい。周囲の小さな山（サイドローブ）が高いほど，離れた振動数からの漏れが大きい。矩形窓はメインローブが最も狭いが，サイドローブが高い。Hann は両者の釣り合いがよく，常時微動の解析で最もよく用いられる。Hamming は最初のサイドローブが低いが，遠方のサイドローブの減衰は遅い。Blackman はサイドローブが最も低いかわりに，メインローブが最も広い。

表3.1　窓関数の特性

窓関数	メインローブの半幅	最大のサイドローブ	最大の振幅の低下
なし（矩形）	Δf	−13 dB	3.9 dB
Hann	2Δf	−31 dB	1.4 dB
Hamming	2Δf	−43 dB	1.8 dB
Blackman	3Δf	−58 dB	1.1 dB

「最大の振幅の低下」は，正弦波の振動数が振動数の刻みのちょうど中間にあるときの振幅の目減りである。式(2.1)のとおり窓の総和で正規化しているため，振動数の刻みに一致する正弦波の振幅は，窓によらず正しく求まる。窓を適用すると実効的な区間長が短くなるため，分解能はやや下がる。地震記録のように区間の両端でもともと振幅が小さい過渡的な波形では，窓を適用せずに全長で1回とすることが多い。

3.4 スペクトルを滑らかにする（⑤）

「Parzen窓」を選んだ場合，帯域幅 b に対する重みを式(3.12)，式(3.13)に示す^[3-4]。重みの総和が 1 となるよう正規化しているため，平滑化の前後でスペクトルの面積は変わらない。

(3.12) (3.13)

ここで，

• W(f)：Parzen 窓の重み。この重みで振動数方向に移動平均をとる
• b：帯域幅。大きいほど滑らかになる（Hz）
• u：帯域幅から決まる係数

「Hann窓」を選んだ場合の重みは式(3.14)であり，b は重みが 0 となるまでの片側の幅である。

(3.14)

平滑化は見やすさのための処理であり，帯域幅を大きくすると鋭いピークは低く広がる。卓越振動数を読む場合は，平滑化の前後でピークの位置が大きく変わらないことを確認するとよい。

3.5 よく使う組合せ

「よく使う組合せ」の常時微動と地震記録は既定値の入口であり，押すと②〜④に典型的な値が入る（常時微動：4096点・50%重なり・Hann窓・選別あり，地震記録：全長で1回・窓なし・選別なし）。

• [3-1]Butterworth, S.: On the theory of filter amplifiers, Experimental Wireless and the Wireless Engineer, Vol. 7, pp. 536-541, 1930
• [3-2]Oppenheim, A. V. and Schafer, R. W.: Discrete-Time Signal Processing, 3rd ed., Pearson, 2010
• [3-3]Harris, F. J.: On the use of windows for harmonic analysis with the discrete Fourier transform, Proceedings of the IEEE, Vol. 66, No. 1, pp. 51-83, 1978, https://doi.org/10.1109/PROC.1978.10837
• [3-4]大崎順彦：新・地震動のスペクトル解析入門，鹿島出版会，1994

4. 地震応答スペクトル

地震応答スペクトルは，固有周期の異なる多数の1自由度系に同じ地震動を入力し，各系の最大応答を固有周期の関数として並べたものである。建物を1自由度系とみなしたとき，その地震動によってどの程度揺れるかを周期ごとに読み取るために用いる^[4-1]。

4.1 1自由度系の運動方程式

固有円振動数 ω，減衰定数 h の線形1自由度系に地動加速度 ÿ_g が作用するとき，地盤に対する相対変位 u の運動方程式は式(4.1)で表される。減衰は粘性減衰とする。

(4.1) (4.2)

ここで，

• u：地盤に対する相対変位（cm）
• ω：固有円振動数（rad/s）
• h：減衰定数
• ÿ_g：地動加速度（cm/s²）
• T：固有周期（s）

4.2 数値解法

式(4.1)を，地動加速度が各時間刻みの間で直線的に変化するとみなして厳密に解く^[4-2]。相対変位と相対速度を並べた状態量を用いると，時刻 j から j+1 への関係は式(4.3)となる。係数行列 A，B₁，B₂ は ω，h，Δt だけで決まり，行列指数関数から求める。この解法は時間刻みの大きさによらず安定であり，外力を直線とみなすこと以外の近似を含まない。

(4.3)

ここで，

• A：状態遷移行列
• B₁，B₂：時刻 j，j+1 の地動加速度に対する係数
• Δt：時間刻み（s）

固有周期が短い系では，時間刻みに対して応答の変化が速く，最大値を取りこぼすおそれがある。そこで，時間刻みが固有周期の 1/20 を超える場合は，式(4.4)の分割数 m で入力を線形補間して細かくしてから解く。たとえば Δt = 0.01 s の記録では，T < 0.2 s の周期で補間する。計算は記録の終わりまでとし，記録が終わった後の自由振動は含まない。

(4.4)

ここで，

• m：時間刻みの分割数（小数点以下を切り上げた整数）
• Δt'：補間後の時間刻み（s）

4.3 求める量

変位応答 S_d は相対変位の絶対値の最大値であり，擬似速度 _pS_v と擬似加速度 _pS_a は S_d から式(4.6)，式(4.7)で換算する。

(4.5) (4.6) (4.7)

ここで，

• S_d：変位応答スペクトル。1自由度系の相対変位の最大値（cm）
• _pS_v：擬似速度応答スペクトル（cm/s）
• _pS_a：擬似加速度応答スペクトル（cm/s²）

擬似速度と擬似加速度は，相対速度や絶対加速度の最大値そのものではない。減衰がない場合，_pS_a は絶対加速度の最大値に一致し，_pS_v は最大のひずみエネルギーに等価な速度となる。減衰定数が 5% 程度であれば，_pS_a は絶対加速度の最大値とほぼ等しい。一方，周期が長い範囲では，_pS_v は相対速度の最大値を下回る傾向がある。

4.4 入力項目

「2. 計算」の「地震応答スペクトルの条件」で，次の項目を指定する。

表4.1　地震応答スペクトルの入力項目

項目	内容	既定値
求める量	pSv，pSa，Sd から1つを選択する。	pSv
減衰定数 h	0〜99 %。建築では 5 % を用いることが多い。	5 %
周期の範囲	計算する最小周期と最大周期（0 < 最小 < 最大）。	0.02〜5 s
周期の刻み	計算する周期の点数（10〜2000）と，等間隔・対数等間隔の別。	200 点，等間隔
水平2成分の合成	NS・EW を同時に解いた合成（4.5 節）も描画する。	なし
トリパータイト	3つの量を1枚で読む図で描画する（4.6 節）。	なし

対数等間隔では，周期を式(4.8)で与える。周期の範囲が広い場合や横軸を対数で見る場合は，短周期側の点が粗くならないよう対数等間隔とするとよい。点数を増やすとスペクトルは滑らかになるが，計算に時間がかかる。

(4.8)

ここで，

• T_min，T_max：周期の範囲の下限と上限（s）
• N：周期の点数

4.5 水平合成

「水平2成分（NS・EW）の合成も作成する」を選択すると，観測点コードが同じ NS・EW の組を探し，同じ周期の1自由度系に2成分を同時に入力して，時刻ごとの相対変位のベクトル和の最大値を式(4.9)で求める。_pS_v，_pS_a は，この値から式(4.6)，式(4.7)で換算する。各成分の最大値を別々に求めて二乗和平方根をとる方法とは異なり，最大値の生じる時刻の違いを考慮した値であり，水平面内のあらゆる方向の応答の最大値に一致する。

(4.9)

ここで，

• S_d,H：水平合成の変位応答スペクトル（cm）
• u_NS，u_EW：NS・EW 成分に対する相対変位（cm）

4.6 トリパータイト

式(4.6)，式(4.7)より _pS_v = ωS_d = _pS_a/ω であるから，横軸を周期，縦軸を _pS_v とした両対数の図では，右下がりの斜めの直線上で S_d が，右上がりの斜めの直線上で _pS_a が一定となる。これらの直線を目盛りとして書き込むと，1枚の図から3つの量を同時に読み取れる（図5.3）。

4.7 前処理

地震応答スペクトルでは，3章の①（平均値の差し引き，倍率およびフィルタ）だけを波形に施し，区間の分割・窓・平滑化は用いない。記録の全長を入力として計算する。

• [4-1]柴田明徳：最新 耐震構造解析 第3版，森北出版，2014
• [4-2]Nigam, N. C. and Jennings, P. C.: Calculation of response spectra from strong-motion earthquake records, Bulletin of the Seismological Society of America, Vol. 59, No. 2, pp. 909-922, 1969, https://doi.org/10.1785/BSSA0590020909

5. 例題

El Centro 1940 NS（Δt = 0.02 s，2687点）を1列の CSV として読み込み，サンプリング周波数に 50 を入力した例を示す。条件は初期値（ハイパス 0.05 Hz，Butterworth 4次，ゼロ位相）のままである。

表5.1　最大値

量	最大値	備考
加速度	341.2 cm/s2	フィルタ後
速度	34.71 cm/s	加速度を1回積分
変位	10.70 cm	加速度を2回積分

「区間の長さ」を全長で1回とすると，全長 2687 点に収まる最大の2のべき乗である 2048 点が用いられる。窓なし・選別なしで求めたフーリエ振幅を，Parzen 窓（帯域幅 0.3 Hz）で平滑化した結果を図3.8に示した。平滑化後の卓越振動数は 1.15 Hz（周期 0.87 秒）である。

応答スペクトルは，周期を 0.02〜10 s の対数等間隔 600 点として求めた。周期の点が粗いと，鋭いピークを取りこぼして最大値を低く見積もることがある（たとえば既定の等間隔 200 点では，h = 2% のピークは 0.7 % 程度低くなる）。

表5.2　擬似速度応答スペクトルのピーク

減衰定数 h	ピーク周期 (s)	pSv (cm/s)
2%	0.86	118.5
5%	0.58	80.4
10%	0.58	66.1

減衰定数を大きくするとピークは低くなり，スペクトルの凹凸もなだらかになる。同じスペクトル（h = 5%）をトリパータイトで描いたものを図5.3に示す。

トリパータイトでは，1枚の図から S_d，_pS_v，_pS_a の3つを同時に読み取れる。短周期側では _pS_a 一定（加速度が支配的），長周期側では S_d 一定（変位が支配的）の傾向が読み取れる。

6. 結果の表示と出力

• 軸：図の上の欄で，横軸・縦軸の上限と下限，および線形・対数を指定する。空欄の項目はデータに合わせて自動で決定する。指定はタブごとに保たれる。
• ピークの印：図に付ける印を，表示しない・1位だけ・1位だけ（数値も）・上位5つ から選択する。表には常に上位5つを示す。時刻歴波形では，加速度・速度・変位の各段で絶対値が最大の点に印を付け，「数値も」では値と時刻を書き込む。「上位5つ」では絶対値の大きい山を5つ（近すぎる山は除く）示し，最大のものに値と時刻を書き込む。
• 位相：フーリエ位相スペクトルと伝達関数の位相差の縦軸を，−180〜180° と 0〜180°（絶対値）から選択する。目盛りは 90° 刻みとする。
• 図の保存：各図の下のボタンから PNG または SVG で保存する。図の文字の大きさと縦横比は画面上部で変更する。
• CSV：一覧の下の「CSV出力」で，計算した結果をすべて書き出す。時刻歴は時間と各波形の加速度・速度・変位，フーリエは振動数と選んだスペクトル，応答スペクトルは周期と各成分の値，伝達関数は組ごとの振幅比・位相差・コヒーレンス（transfer_function.csv），波形の特徴は表の値と計測震度を，減衰は波形ごとの方法・振動数・減衰定数（damping.csv）を出力する。横軸（時間・振動数）が波形によって異なる場合は，波形ごとに横軸の列を付ける。オービットは別のボタンで出力する。
• 設定：画面上部から，入力条件を .prsm 形式で保存および読込できる。一覧（名前・成分・周波数・倍率・単位・選択欄）と波形データも含めて保存し，開くとそのまま戻る。波形データは倍精度のまま圧縮して入れるため，記録の長さと数によってはファイルが数 MB になる。Web 上で開いた Chrome または Edge では，「上書き保存」により同一のファイルに書き戻すことができる。その他の環境では，保存のたびにファイルをダウンロードする。

7. 前提条件と適用範囲

• フーリエ振幅は，窓の総和で正規化した各調和成分の振幅であり，単位は元の波形と同じである。パワースペクトルは各成分の平均パワーであり，パワースペクトル密度とは Δf で割るかどうかが異なる。
• 「計算する」を押すたびに，読み込んだ値から平均値の差し引き，倍率，フィルタの順で処理し直す。遮断周波数を変えて押し直しても，前のフィルタは残らない。
• 速度・変位は積分のたびに同じフィルタを適用し直すため，変位には同じフィルタが3回適用される。低周波側のずれを抑制するための処理である。
• コヒーレンスは，区間が1つしか取れない場合は常に 1 になる。意味のある値を得るには②で区間を分割する必要がある。
• 伝達関数は，組の入力に対する出力の推定量である。区間の取捨は入力の波形の二乗平均平方根で決定する。
• 応答スペクトルは線形1自由度系に対するものであり，減衰は粘性減衰とする。
• 減衰は，1つの振動数の成分を粘性減衰の1自由度系とみなして求める。近い振動数に別の山があると正しく求められない。減衰定数が帯域の半分（既定 5%）に近い場合は，フィルタが減衰の波形を変えるため帯域を広げる。
• RD法の結果は平均した回数によってばらつき，数百回以上を目安とする。対数減衰率は自由振動の区間に用いるものであり，常時微動にそのまま用いると意味のない値となる。ハーフパワー法は区間の長さによって結果が変わる。

用語索引

本文に出てくる主な用語をまとめた。英数字で始まる語をABC順に，続けて日本語の語を五十音順に並べている。日本語の語には対応する英語を括弧で添えた。説明の末尾の括弧は，その用語を扱っている章である。

Parzen 窓（Parzen window）

スペクトルを振動数方向に平滑化するための重み。帯域幅で滑らかさを決定する。（3. 計算条件）

RD法（random decrement method）

不規則な振動を，同じ条件を満たす時刻からそろえて平均し，自由振動に相当する波形を取り出す方法。（2. 求めるもの）

Welch 法（Welch method）

波形を区間に分割し，窓を適用して求めたスペクトルを平均する方法。ばらつきが小さくなる。（2. 求めるもの）

アンラップ（phase unwrapping）

位相スペクトルが ±180° をまたぐときに生じる不連続をつなぎ，連続した値に直すこと。位相の傾きから時間の遅れを読み取るときに用いる。（2. 求めるもの）

因果フィルタ（causal filter）

波形を前向きに1回だけ通すフィルタ。過去の値だけを用いるため実時間の処理と同じであるが，位相が遅れて波形のピークの時刻がずれる。（3. 計算条件）

オービット（orbit）

NS と EW を同時に描いた水平面内の軌跡。（2. 求めるもの）

擬似加速度 _pS_a（pseudo acceleration）

変位応答スペクトルに円振動数の2乗を掛けた量。減衰が小さければ絶対加速度の最大値にほぼ等しい。（4. 地震応答スペクトル）

擬似速度 _pS_v（pseudo velocity）

変位応答スペクトルに円振動数を掛けた量。エネルギーに対応する指標として用いられる。（4. 地震応答スペクトル）

計測震度（JMA seismic intensity）

3成分の加速度に周期特性のフィルタをかけ，ベクトル和から求める気象庁の震度の指標。（2. 求めるもの）

減衰定数（damping ratio）

1自由度系の減衰の大きさを，臨界減衰に対する比で表した値。応答スペクトルでは 5% とすることが多い。（4. 地震応答スペクトル）

コヒーレンス（coherence）

2つの波形の関係が線形で安定しているかを 0〜1 で表す指標。（2. 求めるもの）

サンプリング周波数（sampling frequency）

1秒間に記録した標本の数（Hz）。時間刻み Δt の逆数である。（1. データの読込）

地震応答スペクトル（earthquake response spectrum）

固有周期を変えた多数の1自由度系について，最大応答を周期の関数として並べたもの。（4. 地震応答スペクトル）

状態遷移行列（state transition matrix）

1自由度系の変位と速度を1時間刻み先へ進める行列。応答を厳密に求めるために用いる。（4. 地震応答スペクトル）

水平合成（horizontal vector resultant）

NS と EW の同じ時刻のベクトル和をとり，その最大値を用いる方法。（4. 地震応答スペクトル）

ゼロ位相フィルタ（zero-phase filter）

前向きと後ろ向きに1回ずつ適用し，位相のずれを無くしたフィルタ。（3. 計算条件）

双一次変換（bilinear transform）

アナログフィルタをデジタルフィルタに置き換える変換。振動数軸のゆがみはプリワープで補正する。（3. 計算条件）

帯域幅（band width）

平滑化で平均する振動数の幅。広いほど滑らかになるが，鋭いピークは低く，幅広くなる。（3. 計算条件）

対数減衰率（logarithmic decrement）

自由振動の隣り合う山の振幅の比の自然対数。減衰定数に換算できる。（2. 求めるもの）

台形則（trapezoidal rule）

数値積分の方法の一つ。隣り合う2点の値を直線で結び，その下にできる台形の面積を順に足し合わせて積分を近似する。加速度から速度を，速度から変位を求めるときに用いる（図1.1）。（1. データの読込）

卓越振動数（predominant frequency）

フーリエ振幅が最大となる振動数。地盤や建物のゆれやすい振動数の目安となる。（2. 求めるもの）

中央差分（central difference）

数値微分の方法の一つ。ある時刻の前後1点ずつの値の差を，その間の時間で割って傾きとする。前後を対称に用いるため，片側だけを用いる差分より誤差が小さい。（1. データの読込）

通過域のリップル（passband ripple）

Chebyshev I の通過域に生じる振幅の波打ち。（3. 計算条件）

伝達関数（transfer function）

入力と出力の組について，入力に対する出力の応答倍率と位相差。（2. 求めるもの）

トリパータイト表示（tripartite plot）

変位・擬似速度・擬似加速度を1枚の対数グラフで同時に読めるようにした表示。（4. 地震応答スペクトル）

ドリフト（drift）

積分を重ねることで速度や変位が一方向にずれていく現象。加速度に含まれるわずかな低振動数の成分が原因であり，ハイパスフィルタで抑制する。（3. 計算条件）

ナイキスト振動数（Nyquist frequency）

サンプリング周波数の半分の振動数。これより高い成分は正しく表せないため，フーリエスペクトルはこの振動数までを示す。（2. 求めるもの）

二乗平均平方根（root mean square）

値の2乗の平均の平方根。波形の大きさの目安であり，本アプリケーションでは区間の選別に用いる。（3. 計算条件）

ハーフパワー法（half-power bandwidth method）

スペクトルの山の高さの 1/√2（パワーで半分）となる2つの振動数の幅から減衰定数を求める方法。（2. 求めるもの）

ハイパスフィルタ（high-pass filter）

指定した振動数より低い成分を除去するフィルタ。積分によるドリフトを抑制する。（3. 計算条件）

パワースペクトル密度（power spectral density）

1 Hz あたりのパワーを表したもの。振動数について積分すると波形の二乗平均に等しい。区間の長さが変わっても値の水準が変わらない。（2. 求めるもの）

フィルタ（filter）

波形から特定の振動数の成分を取り除く，または残す処理。本アプリケーションでは，ローパス（遮断周波数より高い成分を除去），ハイパス（低い成分を除去）およびバンドパス（ある帯域だけを残す）を用いる（図3.1）。（3. 計算条件）

フィルタの次数（filter order）

フィルタの極の数。大きいほど遮断が急になるが，位相の変化も大きくなる。（3. 計算条件）

フーリエスペクトル（Fourier spectrum）

波形を種々の振動数の正弦波の重ね合わせとみなし，その振幅を振動数ごとに示したもの。（2. 求めるもの）

プリワープ（prewarping）

双一次変換では振動数軸が高振動数側で圧縮されるため，設計したとおりの遮断周波数にならない。これを避けるため，変換の前に遮断周波数をあらかじめずらしておく操作（図3.2）。（3. 計算条件）

分解能（resolution）

近い2つの振動数を別のものとして見分けられる細かさ。区間の長さの逆数で決まり，区間を長くするほど細かくなる。（3. 計算条件）

平滑化（smoothing）

フーリエスペクトルの細かい凹凸をならす処理。隣り合う振動数の値に重みを付けて平均する。（3. 計算条件）

窓関数（window function）

区間の端を滑らかに 0 に近づける重み。切り出しによる漏れを抑制する。（3. 計算条件）

メインローブとサイドローブ（main lobe and side lobes）

窓関数の周波数特性のうち，中央の山と周囲の小さな山。前者は分解能に，後者は漏れに関係する。（3. 計算条件）

漏れ（leakage）

波形を有限の長さで切り出すことにより，本来その振動数にない成分が周囲の振動数に広がる現象。窓関数で抑制する。（3. 計算条件）