Web-Wave（模擬地震波作成プログラム）

概要

Lilhanand & Tseng (1988) による反復フーリエ振幅修正法で模擬地震波を生成します。初期波形のフーリエ振幅を目標応答スペクトルに合うよう繰り返し修正する手法です。

アルゴリズム

1. 継続時間 T、時刻刻み Δt を設定し、データ点数 M = T/Δt を決定する。FFT点数 N は M 以上の最小の2のべき乗。
2. 乱数位相 φ_k ∈ [0, 2π) を生成し、初期フーリエ係数を構成する。
   Ck = Ak · eiφk
3. IFFT して時刻歴波形 a(t) を得る。
4. 包絡線関数 e(t) を乗じる（Jennings 型）。
5. 再度 FFT し、フーリエ振幅 |C_k| を取り出す。
6. 各周期 T_j で応答スペクトル SA_j を計算し、目標値 S_0,j との比率 r_j = S_0,j / SA_j を求める。
7. フーリエ振幅を周波数方向に補間した r(f) で修正する：
   |Ck|new = |Ck|old × r(fk)
8. 修正した振幅と元の位相から時刻歴を再構成し、ベースライン補正を行う。
9. 3〜8 を収束するまで繰り返す（最大繰り返し回数・収束判定ε）。

収束判定

各照査周期で max |1 − SAj/S0,j|< ε（ε: 収束判定値[%]/100）を満たした時点で収束と判定します。

参考文献

Lilhanand, K. & Tseng, W.S. (1988): Development and application of realistic earthquake time histories compatible with multiple-damping design spectra. Proc. 9th WCEE, Vol.II, 819–824.

離散フーリエ変換（DFT）の定義

N 点の離散時刻歴 {x_n} に対するDFTとその逆変換（IDFT）：

X_k = Σ_n=0^N−1 x_n · e^{−i·2π·k·n/N}   （順変換）
x_n = (1/N) · Σ_k=0^N−1 X_k · e^{+i·2π·k·n/N}   （逆変換）

FFT の実装（Cooley-Tukey 基数2）

計算量を O(N²) → O(N log₂N) に削減します。N は2のべき乗。ビット逆順並べ替え後にバタフライ演算を繰り返します。

周波数分解能とナイキスト周波数

• 周波数刻み：Δf = 1 / (N·Δt) [Hz]
• ナイキスト周波数：f_Ny = 1 / (2·Δt) [Hz]
• k 番目のビン周波数：f_k = k·Δf (k = 0, 1, …, N/2)

フーリエ振幅と地盤増幅

地盤増幅を適用する場合、フーリエ振幅に伝達関数 |H(f)| を乗じます：
|Ck|amplified = |Ck| × |H(fk)|

1自由度系の運動方程式

ẍ(t) + 2h·ω₀·ẋ(t) + ω₀²·x(t) = −ü_g(t)

ここで h：減衰定数、ω₀ = 2π/T₀：固有円振動数、ü_g：地盤加速度。

Nigam & Jennings (1969) 逐次積分法

各時刻ステップ [t_n, t_n+1] で地盤加速度を線形補間し、解析解を用いて状態変数 {x, ẋ} を伝播します：

{x_n+1, ẋ_n+1} = A·{x_n, ẋ_n} + B·{ü_g,n, ü_g,n+1}

行列 A, B の各成分は ω₀, h, Δt の関数として解析的に与えられます（厳密な逐次解）。

応答スペクトル値の定義

• 加速度応答スペクトル：SA(T) = max|ẍ(t) + ü_g(t)| = max|ω₀²·x(t) + 2h·ω₀·ẋ(t)|
• 速度応答スペクトル（擬似）：pSV(T) = SA(T)·T/(2π)
• 変位応答スペクトル（擬似）：pSD(T) = SA(T)·(T/(2π))²

参考文献

Nigam, N.C. & Jennings, P.C. (1969): Calculation of response spectra from strong-motion earthquake records. Bull. Seism. Soc. Am., 59(2), 909–922.

告示スペクトル（平12建告第1461号）

設計用入力地震動の標準スペクトル（減衰5%）：

Z：地域係数（0.7〜1.0）。稀地震は上式を Z/5 倍。

簡略法 地盤増幅率 Gs

告示に基づく地盤種別ごとの Gs 倍率を目標スペクトルに乗じます。

• 第1種地盤 (式4a): T < Tg では線形増加、T ≥ Tg では一定
• 第2種地盤 (式4b, gv = 2.025): 特定周期帯で増幅
• 第3種地盤 (式4b, gv = 2.7): より長周期側で増幅

直接入力スペクトル

任意の周期−スペクトル値をテーブル入力し、照査周期でログ補間します。

Jennings 型包絡線

模擬地震波に時間的な振幅変動を与えるため、以下の包絡線 e(t) を乗じます：

e(t) = (t/t₁)²    (0 ≤ t < t₁)  【上昇部】
e(t) = 1           (t₁ ≤ t < t₂)  【定常部】
e(t) = exp(−c·(t−t₂))    (t₂ ≤ t)  【減衰部】

• t₁：上昇終了時刻
• t₂：定常終了時刻
• c：減衰係数（= −ln(e_end) / (T − t₂)、e_end：終端振幅比）

包絡線適用後に再FFTして振幅修正を継続するため、収束後の位相特性に包絡線の形状情報が織り込まれます。

乱数位相

線形合同法または xorshift 法による疑似乱数列を用いて位相 φ_k ∈ [0, 2π) を生成します。乱数シードを固定することで再現性のある波形を得られます。シード = 0 のとき現在時刻をシードとして自動設定します。

実地震波位相

観測地震波の時刻歴をFFTし、各フーリエ成分の位相角 φ_k = arg(C_k) を抽出して模擬波の位相として利用します。これにより実地震波に近い位相特性を持つ模擬波形が得られます。

処理手順（反復ループ内・加速度補正）

1. 速度・変位の線形加速度積分：補正前の加速度を台形則で積分し v(T)、x(T) を得る。
2. 拘束条件を満たす線形補正量 (c₀ + c₁t) を決定：
   v'(T) = 0、x'(T) = 0 を同時に満たす c₀, c₁ を連立方程式で解く。
3. 加速度補正：a'(t) = a(t) − (c₀ + c₁t)

ベースライン補正は各反復ステップで実施され、速度・変位の発散を防ぎます。

概要

加速度をフーリエ変換し、周波数ごとに iω で除算して速度・変位スペクトルを求め逆変換する方法。時間域の数値積分に生じる誤差累積がなく、各周波数成分を独立に処理できる。

計算式

A(ω) = FFT{ a(t) }
V(ω) = A(ω) / (iω)　　→　Re[V] = Im[A]/ω, Im[V] = −Re[A]/ω
D(ω) = A(ω) / (−ω²)　→　Re[D] = −Re[A]/ω², Im[D] = −Im[A]/ω²
v(t) = IFFT{ V(ω) },　d(t) = IFFT{ D(ω) }

k=0（直流成分）は 0 に設定して平均ドリフトを除去。負周波数成分（k>N/2）は符号付き ω=2π(k−N)Δf で処理し、IFFT 出力を実数に保つ。

参考文献

Trifunac, M.D. & Lee, V.W. (1973): Routine computer processing of strong-motion accelerograms. EERL 73-03, Caltech.
Boore, D.M. (2005): On pads and filters: Processing strong-motion data. Bull. Seism. Soc. Am., 95(2), 745–750.

概要

固有周期 T_p・減衰定数 h=1/√2 の減衰 SDOF（振り子）が地盤加速度を受けたときの相対変位・相対速度を Nigam 法（厳密な状態遷移行列解）で求め、地盤変位・速度に変換する方法。

運動方程式

ẍ_r + 2h·ω₀·ẋ_r + ω₀²·x_r = −ü_g(t)　　(h = 1/√2, ω₀ = 2π/T_p)

Nigam 法による厳密解（ステップ間で ü_g を線形補間）

ωd = ω₀·√(1−h²), E = exp(−h·ω₀·Δt)
C = cos(ωd·Δt), S = sin(ωd·Δt)
Δü_g = ü_g,n+1 − ü_g,n

x_n+1 = E·(C+h·ω₀·S/ωd)·x_n + E·S/ωd·ẋ_n
    + [(E·C−1)/ω₀² + E·h·S/(ω₀·ωd)]·ü_g,n
    + [2h·(1−E·C)/(ω₀³·Δt) + E·(1−2h²)·S/(ω₀²·ωd·Δt) − 1/ω₀²]·Δü_g

ẋ_n+1 = −E·ω₀²·S/ωd·x_n + E·(C−h·ω₀·S/ωd)·ẋ_n
    − E·S/ωd·ü_g,n
    + (a₁₁−1)·Δü_g/(ω₀²·Δt)

地盤速度・変位への変換

振り子の相対応答と地盤運動の関係（初期静止条件 x_r(0)=0, ẋ_r(0)=0）：

v(t) = −ẋ_r(t)
d(t) = −x_r(t)

T_p→∞（ω₀→0）の極限では台形則による直接積分に一致する。T_p は波形継続時間より十分長く（通常 40s 以上）設定すること。

参考文献

Nigam, N.C. & Jennings, P.C. (1969): Calculation of response spectra from strong-motion earthquake records. Bull. Seism. Soc. Am., 59(2), 909–922.

重複反射理論（SH波鉛直入射）

N 層地盤の複素剛性を用いた SH 波の伝達関数（地表加速度 / 工学的基盤入力加速度）を求めます。各層 m の複素せん断弾性係数：

G*_m = G_m · (1 + 2i·h_m)

複素せん断波速度：Vs*_m = √(G*_m / ρ_m）

各層の伝達行列を全層にわたって積算して地表−基盤間の伝達関数 H(f) を得ます：

H(f) = u_surface(f) / u_bedrock(f)

等価線形化の反復手順

1. 各層の初期せん断剛性 G₀ = ρ·Vs₀²、初期減衰 h₀ を設定する。
2. 入力スペクトルに H(f) を乗じて各層の応力フーリエスペクトルを計算する。
3. 各層の最大せん断歪 γmax を計算する（パワースペクトル積分の平方根）。
4. 有効歪 γeff = α · γmax（α：有効歪比率、通常 0.65）。
5. G/G₀ カーブおよび h−γ カーブから γeff に対応する G, h を読み取る。
6. 更新した G, h で次の反復を行う。
7. 全層の G, h の変化率が収束判定 ε 以下になるまで繰り返す。

G/G₀ モデル（Hardin-Drnevich 型）

G/G₀ = 1 / (1 + γ/γ_ref)

粘土系：γ_ref = 0.0005（基準歪）、砂系：γ_ref = 0.0003（目安値）

減衰 h−γ 関係（Hardin-Drnevich 型）

h(γ) = h_max · (1 − G/G₀)

h_max：最大減衰定数（粘土系 0.20、砂系 0.25 程度）

伝達倍率と設計への適用

収束後の |H(f)| を模擬波作成時のフーリエ振幅修正に用いることで、地盤増幅を直接反映した模擬地震波を生成します。設計スペクトルが地盤増幅込みの場合は地盤増幅をOFFにし、工学的基盤スペクトルに基づく場合は地盤増幅をONにします。

1. 入力パラメータ（Δt, T, h, 目標スペクトル, 包絡線, 位相設定）を確認
2. （地盤増幅 ON の場合）等価線形化法で伝達関数 |H(f)| を計算
3. 乱数位相または実地震波位相を生成
4. 初期フーリエ係数 = 一様振幅 × 位相
5. IFFT → 包絡線適用 → FFT（振幅抽出）
6. （地盤増幅 ON の場合）|H(f)| を振幅に乗算
7. 各照査周期で応答スペクトル計算 → 振幅比率 r(f) 計算
8. フーリエ振幅を r(f) で修正
9. IFFT → ベースライン補正 → 収束判定
10. 収束するまで 5〜9 を繰り返す
11. 最終波形を出力（加速度・速度・変位時刻歴、応答スペクトル）