地震や建築応答の時刻歴波形が加速度データで与えられているとき、それを速度や変位データに変換する場合、数値積分を行う必要があります。ここでは、書籍「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)のうち、「平成７年（1995年）兵庫県南部地震」の神戸中央区中山手(JMA KOBE)における加速度時刻歴データを変位データに変換してみます。もとの加速度時刻歴波形は下図の通りです。

これを、前述の積分計算プログラムを用いて2階積分を行い、変位時刻歴データを算出した結果は下図の通りとなります。長方形近似、台形近似およびシンプソン公式のいずれの方法でも、同じような波形となります。

x方向、z方向の最大変位は数メートルのオーダ、y方向に至っては50m以上の最大変位が発生していますが、現実の最大観測変位は25cm程度です。この大きな変位は、低振動数成分のノイズや誤差が存在することにより発生する、積分波形のドリフトによるものです。

このため適切な基線補正を施すことになりますが、補正の方法については別記事で改めて紹介します。