Random Decrement (RD法) によるシステム同定の理論と実装

概要
Random Decrement (RD法) とは
RD法の理論
RD法の実装
参考文献

概要

構造ヘルスモニタリングにおけるシステム同定の手法として用いられる、Random Decrement (RD法) の理論と実装例について整理する。

Random Decrement (RD法) とは

常時微動や風応答などランダムな外力を受ける構造物の振動データから、減衰自由振動波形を得るために利用される手法である。得られた減衰自由振動波形は、自由振動解とのカーブフィッティングや Eigensystem Realization Algorithm (ERA) によってモード特性を同定するために用いられる。 RD法はランダム応答の期待値がゼロであることを利用した手法で、アルゴリズムも直感的に理解しやすく、使い勝手の良い手法である。

RD法の理論

1質点系の運動方程式は、固有角振動数 $\omega_0$ と減衰比 $\zeta$ を用いて次のように表される。ここでは、期待値をゼロとするランダムな外力 $F(t)$ が作用するものとする。 $$ \ddot{X} + 2\zeta\omega_0\dot{X} + \omega_0^{2}X = F(t) \tag{1} $$

この場合、応答 $X$ の一般解は自由振動解 $D(t)$ とランダム応答 $R(t)$ で表される。 $$ X = D(t) + R(t) \tag{2} $$

応答 $X$ の期待値は次のようになる。 $$ E[X] = E[D(t)] + E[R(t)] \tag{3} $$

ランダム応答の期待値 $E[R(t)]$ はゼロであるため、次のようになる。 $$ E[X] = E[D(t)] \tag{4} $$

ここで、応答 $X$ について十分な数の標本 $X_i$ が得られ、各標本の極大値の位置を合わせて重ね合わせた時、ランダム応答の和 $\sum R_i(t)$ はゼロに近づくため、自由振動解の和 $\sum D_i(t)$ のみが残る。

$$ \sum X_i(t) \simeq \sum D_i(t) \tag{5} $$

以上により、十分な標本数の応答 $X$ の重ね合わせにより、系の減衰自由振動波形が得られる。実問題においてはバンドパスフィルタを用いて特定のモード成分が残るように前処理を行い、処理後の波形にRD法を適用することで、精度良く減衰自由振動波形を求めることができる。

また、自由振動解 $D(t)$ は下記のように表されるため、RD法により得られた減衰自由振動波形とカーブフィッティングを行うことで、系の固有角振動数 $\omega_0$ および減衰比 $\zeta$ を得ることができる。

$$ \begin{gather} D(t) = A \exp(-\zeta\omega_0t)\cos(\sqrt{1-\zeta^{2}}\omega_0t-\phi) \tag{6} \\ A,\phi: const \end{gather} $$

RD法の実装

PythonにてRD法の簡単な使用例を実装する。

環境

下記の外部ライブラリを使用する。

numpy
matplotlib
control
scipy

ソースコード

ソースコードの大まかな処理の流れとしては次の通り。

1質点系のモデルに乱数によるランダムな外力を入力して擬似的な観測波形を生成
観測波形の最大値の後の部分を切り出す
切り出した波形を重ね合わせ、減衰自由振動波形を生成
生成した減衰自由振動波形と自由振動解からカーブフィッティングを行い、固有角振動数と減衰比を推定
減衰自由振動波形の生成過程が確認できるように、波形の様子をgifアニメにプロットして出力

この実装ではバンドパスフィルタによる前処理は省略した。ノイズ成分も多数含まれる実際の計測波形に適用する場合には、バンドパスフィルタにより前処理を行うことで、少ない波形からでも精度良く減衰自由振動波形を得ることができる。Pythonでは scipy ライブラリの scipy.signal.butter を利用することで、バンドパスフィルタも簡単に実装できる。

ソースコードは下記に示す。

import numpy as np
import matplotlib.animation
import matplotlib.pyplot as plt
from control import matlab
from scipy.optimize import curve_fit

np.random.seed(0)


def free_vibration(t, omega, zeta):
    """
    カーブフィッティング用の自由振動関数
    """
    return np.exp(-zeta * omega * t) * np.cos(np.sqrt(1 - zeta**2) * omega * t)


def main():
    # 観測波形の個数
    obs_num = 100

    # 観測波形の長さ
    obs_len = 2000

    # 切り出し波形の長さ
    trim_len = 400

    # 時間間隔
    dt = 0.025

    # アニメーションの1フレームのミリ秒数
    animation_interval = 200

    # 1質点系モデルの生成
    omega = 2 * np.pi * 1.0  # 固有角振動数
    zeta = 0.02  # 減衰比
    m = 1.0
    k = omega**2 * m
    c = 2 * zeta * m * omega
    sys = matlab.tf([1], [m, c, k])

    # 疑似観測波形の生成
    t = np.arange(obs_len) * dt
    x_obs = np.zeros((obs_num, obs_len))
    for i in range(obs_num):
        x_obs_i, _, _ = matlab.lsim(sys, U=np.random.normal(size=obs_len), T=t)  # ランダム外力による応答
        x_obs[i] = x_obs_i + np.random.normal(scale=0.005, size=obs_len)  # 観測ノイズを付与

    # ピーク後の波形の切り出し
    max_indexes = x_obs[:, : obs_len - trim_len].argmax(axis=1)
    x_trim = np.zeros((obs_num, trim_len))
    for i in range(obs_num):
        x_trim[i] = x_obs[i, max_indexes[i] : max_indexes[i] + trim_len]

    # 切り出し波形を重ね合わせ、自由振動波形を生成
    x_sum = np.sum(x_trim, axis=0)
    x_sum /= x_sum[0]

    # カーブフィッティングにより固有角振動数と減衰比の推定
    t_trim = np.arange(trim_len) * dt
    (omega_hat, zeta_hat), _ = curve_fit(free_vibration, t_trim, x_sum, bounds=(0, (np.inf, 1.0)))
    print("固有角振動数")
    print(f" 理論値: {omega:.4f}")
    print(f" 推定値: {omega_hat:.4f}")
    print("減衰比")
    print(f" 理論値: {zeta:.4f}")
    print(f" 推定値: {zeta_hat:.4f}")

    # --- 以下、gifアニメ描画処理 ---
    fig, axes = plt.subplots(3, 1, figsize=(8, 6))
    axes[0].set_xlim(0, obs_len * dt)
    axes[1].set_xlim(0, trim_len * dt)
    axes[2].set_xlim(0, trim_len * dt)
    artists = []

    # 減衰自由振動の理論解を比較用にプロット
    free_vibration_wave, _ = matlab.initial(sys, X0=1.0, T=t_trim)
    axes[2].plot(t_trim, free_vibration_wave, color="0.7", label="Theoretical value")
    axes[2].legend()

    # 観測波形1つごとにアニメーションのフレームを作成
    for i in range(obs_num):

        # iまでの切り出し波形を重ね合わせ
        x_sum_i = np.sum(x_trim[: i + 1], axis=0)
        x_sum_i /= x_sum_i[0]

        # 各オブジェクトの描画
        trim_region = axes[0].axvspan(t[max_indexes[i]], t[max_indexes[i] + trim_len], fc="0.9", ec="0.7")
        (plot1,) = axes[0].plot(t, x_obs[i], color="dimgray")
        (plot2,) = axes[1].plot(t_trim, x_trim[i], color="dimgray")
        (plot3,) = axes[2].plot(t_trim, x_sum_i, color="tomato", lw=3.0)
        text1 = axes[0].text(0.0, 1.05, f"Observation wave: $X^{{obs}}_{{{i+1}}}$", transform=axes[0].transAxes, fontsize="large")
        text2 = axes[1].text(0.0, 1.05, f"Trimed wave: $X_{{{i+1}}}$", transform=axes[1].transAxes, fontsize="large")
        text3 = axes[2].text(0.0, 1.05, f"Free vibration wave: $\sum\,_{{i=1}}^{{{i+1}}} X_i$", transform=axes[2].transAxes, fontsize="large")
        artists.append([trim_region, plot1, plot2, plot3, text1, text2, text3])

    fig.tight_layout()
    animation = matplotlib.animation.ArtistAnimation(fig, artists, interval=animation_interval)
    animation.save("rd_animation.gif", writer="pillow")


if __name__ == "__main__":
    main()

出力結果

観測波形の個数 obs_num を変えた2つのパターンの結果を確認する。

`obs_num=10` の結果 (`animation_interval=1000`)

まずは10波形を重ね合わせた結果で出力内容を確認する。

Observation wave: 実際に観測される波形にあたるもの。グレーの枠の部分が次のTrimed waveになる。
Trimed wave: Observation waveから最大値の後の部分を切り出した波形である。当然、観測波形のピークから切り出しただけでは、きれいな減衰自由振動波形は得られないことがわかる。
Free vibration wave: Trimed waveを重ね合わせた減衰自由振動波形をオレンジで示している。初期値は1にスケーリングしてある。グレーで示したのは理論値である。重ね合わせの数が増えるにつれてランダム成分、ノイズ成分が除去されて滑らかな減衰自由振動波形に近づいている。

`obs_num=100` の結果 (`animation_interval=200`)

100波形を重ね合わせた結果を示す。50波形程度重ね合わせたあたりで減衰自由振動波形の変化が落ち着き、滑らかな減衰自由振動波形が得られているのが確認できる。理論値と比較すると、振幅の減衰についてはやや乖離があるものの、周期については精度良く一致した値が得られているように見える。

また、カーブフィッティングにより固有角振動数 $\omega_0$ と減衰比 $\zeta$ も得られた。グラフからも確認できたように、減衰比は理論値とやや乖離があるが、固有角振動数は理論値に近い値が得られている。

固有角振動数
 理論値: 6.2832
 推定値: 6.2680
減衰比
 理論値: 0.0200
 推定値: 0.0349

概要