地震や建築応答の時刻歴波形が加速度データで与えられているとき、それを速度や変位データに変換する場合、数値積分を行う必要があります。ここでは、書籍「Excelではじめる数値解析」を参考に、近似的な数値積分をPythonを用いて求める方法を紹介します。
なお、本記事の最後に紹介する、気象庁|強震観測データの数値積分に用いたインプットとソースコードは、以下よりダウンロード可能です。
時々刻々の数値積分
通常、Pythonで数値積分をする際は、SciPyライブラリのモジュールであるintegrateを用います。たとえば、$y={\rm sin}x$を$0$から$\pi$まで定積分するプログラムは、以下の通り記述します。
import numpy as np
from scipy import integrate
y = lambda x: np.sin(x)
y_integral = integrate.quad(y, 0, np.pi)
print(y_integral)
実行結果は以下の通り。積分値と推定誤差が出力されます。
(2.0, 2.220446049250313e-14)
しかし、たとえば時刻歴加速度を積分して時刻歴速度を得たい場合は、積分計算を時々刻々で繰り返し行う必要がありますが、そのような関数はpythonには実装されていません(少なくとも、筆者は知りません)。
長方形近似
ここでは、元の離散的な時刻歴データを、下図のように長方形の短冊に分割し、その面積を時々刻々で足し合わせることで、積分後の時刻歴データを求める方法を紹介します。
被積分の時刻歴データを$f(t)$、積分後の時刻歴データを$F(t)$とすると、時刻$t_n$での積分後の時刻歴データ$F(t_n)$は、$t_n$までの長方形面積の和、すなわち
$$F(t_n)=\sum_{i=0}^{n}f(t_i)\Delta t$$
で表されます。これを、時間刻みdt(浮動小数点型)、各時刻の被積分データdx_dt(リスト型)を引数とするpythonモジュールの形で記述すると、以下の形となります。実行結果は後述します。
def time_hist_integrator_rect(dt, dx_dt):
xt = 0.
x = []
for dxt_dt in dx_dt:
xt += dxt_dt * dt
x.append(xt)
return x
台形近似
ここでは、下図のように台形の面積を時々刻々足し合わせることで、積分後の時刻歴データを求める方法を紹介します。
時刻$t_n$での積分後の時刻歴データ$F(t_n)$は、
$$F(t_n)=\sum_{i=0}^{n}\frac{f(t_{i+1})+f(t_i)}{2}\Delta t$$
で表されます。これを長方形近似と同様に、時間刻みdt(浮動小数点型)、各時刻の被積分データdx_dt(リスト型)を引数とするpythonモジュールの形で記述すると、以下の形となります。
def time_hist_integrator_trap(dt, dx_dt):
xt = 0.
x = []
i_x = len(dx_dt)-1
for i in range(i_x):
xt += (dx_dt[i] + dx_dt[i+1]) * dt / 2.
x.append(xt)
return x
シンプソン公式
シンプソン公式による積分は、2次曲線で離散データを近似し、その面積を求める方法です。
上図のように、3つの点を通る2次曲線と$t$軸に囲まれた部分の面積$S_i$は以下式により求まります(導出は、 書籍「Excelではじめる数値解析」 を参照)。
$$S_i=\frac{\Delta t(f_{i-2}+4f_{i-1}+f_{i})}{3}$$
これを、時間刻みdt(浮動小数点型)、各時刻の被積分データdx_dt(リスト型)を引数とするpythonモジュールの形で記述すると、以下の形となります。なお、このモジュールでは、最初の区間$t_0$~$t_1$の面積$S_1$は、定義点が2つしかなく2次曲線近似ができないので、直線近似による面積で求めています(”if i == 0″の部分)。
def time_hist_integrator_simp(dt, dx_dt):
xt1 = 0.
xt2 = 0.
x = []
i_x = len(dx_dt)-1
for i in range(0, i_x):
if i == 0:
xt1 = dt / 2. * (dx_dt[i] + dx_dt[i+1])
x.append(xt1)
elif i % 2 == 0:
xt2 += dt / 3. * (dx_dt[i-2] + 4*dx_dt[i-1] + dx_dt[i])
x.append(xt2)
else:
xt1 += dt / 3. * (dx_dt[i-2] + 4*dx_dt[i-1] + dx_dt[i])
x.append(xt1)
return x
積分結果
ここでは、気象庁|強震観測データ(https://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/index.html)のうち、「平成7年(1995年)兵庫県南部地震」の神戸中央区中山手(JMA KOBE)における加速度時刻歴データを変位データに変換してみます。もとの加速度時刻歴波形は下図の通りです。
これを、前述の積分計算プログラムを用いて2階積分を行い、変位時刻歴データを算出した結果は下図の通りとなります。長方形近似、台形近似およびシンプソン公式のいずれの方法でも、同じような波形となります。
x方向、z方向の最大変位は数メートルのオーダ、y方向に至っては50m以上の最大変位が発生していますが、現実の最大観測変位は25cm程度です。この大きな変位は、低振動数成分のノイズや誤差が存在することにより発生する、積分波形のドリフトによるものです。
このため適切な基線補正を施すことになりますが、補正の方法については別記事で改めて紹介します。