建築物の振動測定結果解析プログラム Building Vibration Analysis System 常時微動測定の時刻歴から固有振動数・減衰定数・モード形状を同定（実稼働モード解析 OMA） [v260918_r2]

解析手法: ※ 手法ごとに同定結果（振動数・減衰）が異なるため、複数手法で照合することを推奨。

1. 測定情報

測定の概要
サンプリングと物理量

※ 同定される固有振動数・減衰定数・（正規化された）モード形状は、入力データの単位・スケールに依存しません。 単位は時刻歴グラフと計算書の軸ラベル・統計量の表記にのみ使われます。

2. 測定結果（時刻歴）の貼付

時刻歴データ

※ 空セル・非数値セルは欠測として扱い、前後の値で線形補間します（補間箇所は取込結果に表示）。
※ 想定規模は 約 10 万点 × 12 チャンネル。これを大きく超える場合は「3. 前処理」の間引き（デシメーション）を併用してください。

データ未取込
測定場所（センサ配置）
簡易（階・方向）
詳細（3D 座標）
使用 ch 名称 階 方向 高さ
GL+[m]		 平面 X
[m]		 平面 Y
[m]		 参照 入力

※ 同一階・同一方向に複数センサがある場合、平面 X / Y を入力するとねじれ（回転）成分を分離して表示・アニメーション化します。
※ 地表・基礎などで地動も測定している場合は、その ch の「入力」にチェックを入れると入出力同定（EMA）が行われます（ch 名称に GL・地盤・地表 を含む場合は自動でチェック）。

使用 ch 名称 x [m] y [m] z [m] dx dy dz 参照 入力

※ 骨格線は「センサ番号(1始まり)-センサ番号」をカンマ区切りで指定します。空欄の場合は、同じ高さ z のセンサを平面上で凸包順に結び、 さらに平面位置が近いセンサ同士を鉛直方向に結んだ骨格線を自動生成します。

3. 前処理

トレンド除去・間引き・帯域制限

※ 処理順序は ①トレンド除去 → ②帯域フィルタ → ③間引き とします。

4. 解析条件

スペクトル推定（FDD / EFDD 共通）
EFDD / FSDD
SSI
入出力同定（EMA）— 地動（入力）を測定している場合

センサ配置表で「入力」に指定した ch を地動、その他の使用 ch を建物応答として周波数応答関数（FRF）を求め、 有理関数モデルの当てはめで固有振動数・減衰定数・刺激関数 βφ を同定します。入力が無い場合は実行されません。

5. 対象モードの選択（卓越振動数）

まず「⚙ 解析実行」を押すとスペクトルが計算され、卓越振動数が自動検出されます。
その後「スペクトル」タブのグラフをクリックして追加、または下表で直接編集できます（編集後は再度「⚙ 解析実行」）。

No選択振動数 f [Hz] SSI
次数		 EMA
極		 備考削除

⚠ 入力エラー

    ⚠ 注意

      解析結果

      解析未実行 — 時刻歴を貼り付けて「⚙ 解析実行」を押してください。
      同定結果
      時刻歴
      データ診断
      スペクトル
      安定化線図
      入出力（FRF）
      モード形状
      アニメーション
      MAC 行列
      詳細計算ログ

      解析未実行です。

      データ未取込です。

      チャンネルごとの信号の性質を確認します。 常時微動の OMA は「定常・ガウス性の不規則振動」を前提とするため、確率密度が正規分布から大きく外れる、 スペクトログラムに時間的に局在した強い成分がある、といった場合は外乱（歩行・機械振動・衝撃）の混入を疑い、該当区間の除外を検討してください。

      解析未実行です（前処理後の波形を対象とします）。

      STFT スペクトログラム（短時間フーリエ変換の振幅）

      解析未実行です。

      グラフ上をクリックすると、その周波数を「対象モード」に追加します（Shift+クリックで最近傍のピークにスナップ）。

      解析未実行です。

      安定化線図（stabilization diagram）: 横軸=固有振動数、縦軸=モデル次数。 緑○ = Soft Criteria を満たす安定極、灰× = 不安定極。縦に並ぶ安定極の列が真のモードを示します。
      クリック = その振動数をモード追加／Shift+クリック = 最も近い極を選択（振動数と次数を記録。抽出方法「極を個別選択」で使用）。

      SSI 未実行です。

      入出力同定（EMA）: 地動（入力）に対する建物応答（出力）の周波数応答関数から極を求めます。
      安定化線図: 緑 × = 振動数・減衰とも安定な極、○ = 選択した極。クリック = その振動数をモード追加／ Shift+クリック = 最も近い安定極を選択（その場で LSFD を再計算し、結果に反映します）。

      解析未実行です。

      解析未実行です。

      解析未実行です。

      ※ モード形状は複素数のため、φj(t) = Re[ φj · eiωt ] として時間発展を描画します（節点位相差が表現されます）。 変形倍率は表示上の任意スケールであり、実際の振幅を表すものではありません。

      解析未実行です。

      
        
      

      
      
    
      

  
      

  
  
      
    
      使い方・解析ロジック詳細解説
      
    
      
      ▼ 展開して読む
      
      

        
      このアプリについて
      
        
      
          建築物の常時微動（または人力加振・地震観測）測定で得られた加速度・速度・変位の時刻歴から、
          固有振動数・減衰定数・モード形状を同定するアプリです。外力（入力）が未知のまま応答のみから
          同定するため 実稼働モード解析（OMA: Operational Modal Analysis）、あるいは
          出力のみ同定（output-only identification） と呼ばれます。
        
      
        
      
          同定手法は周波数領域の FDD / EFDD / FSDD と時間領域の SSI-cov です。
          各処理の内容は下記の解説に、実際に用いた条件と中間結果は「詳細計算ログ」に示しています。
        
      

        
      使い方（基本の流れ）
      
        
        
          
      1. 「1. 測定情報」に建物と測定の条件を入力します。サンプリング周波数 fs は必須です
            （時刻歴に時間列があれば自動推定できます）。
        2. 
          
      2. 「2. 測定結果（時刻歴）の貼付」の欄に、Excel の範囲をコピーしてそのまま貼り付け、
            「📥 データ取込」を押します。CSV/TXT ファイルからの読込も可能です。
            表形式で確認・編集したい場合は「📊 スプレッドシート入力」を押すと別ウィンドウの表が開きます。
            セルを選択して Ctrl+V で貼り付け（先頭行が見出しなら「名称」行に入ります）、必要なら修正して
            「📥 本体に反映（データ取込）」を押します。
        3. 
          
      3. 取込むとセンサ配置表が自動生成されます。ch 名称から階・方向を推定しますが、
            必ず「階」「方向」「高さ GL+」を確認・修正してください。ここが正しくないとモード形状図と
            アニメーションが正しく描かれません（同定される振動数・減衰は影響を受けません）。
        4. 
          
      4. 「3. 前処理」「4. 解析条件」を確認します。既定値は標準的な設定値です。
            まずは既定値のままで実行し、結果を見ながら調整するのが確実です。
        5. 
          
      5. 「⚙ 解析実行」を押します。卓越振動数が自動検出され、選択した手法（FDD／EFDD・FSDD／SSI-cov）で同定が行われます。
        6. 
          
      6. まず「データ診断」タブで信号の性質（自己相関・PSD・確率密度・正規確率プロット・スペクトログラム）を確認し、
            衝撃的な外乱や周期的な機械振動が混入していないかを見ます。
        7. 
          
      7. 右の「スペクトル」タブで CMIF（特異値スペクトル）を確認し、
            対象としたいピークをクリックして「5. 対象モードの選択」に追加・修正します
            （Shift+クリックで最近傍のピークにスナップ）。修正したら再度「⚙ 解析実行」。
        8. 
          
      8. 「安定化線図」タブで、緑○（安定極）が縦に並んでいる周波数が真のモードです。
            自動検出が拾えていないモードがあればここでクリックして追加します。
            Shift+クリックすると最も近い極が選ばれ、その次数が「5. 対象モードの選択」の表に記録されます
            （抽出方法が自動で「極を個別選択」に切り替わります）。
            「図の種類」で振動数－減衰図・合成スペクトルと実測の比較・Hankel 特異値にも切り替えられます。
        9. 
          
      9. 「モード形状」「アニメーション」タブで、同定されたモードが構造的に妥当な形状か
            （1 次は節なし、2 次は節 1 個、など）を目視で確認します。
            アニメーションは「🎬 動画出力」で MP4（または WebM）ファイルに保存でき、報告・説明資料に使えます
            （現在の表示方向・倍率で記録。記録には動画の長さと同じ時間がかかり、その間はタブを表示したままにします）。
        10. 
          
      10. 地動（地表・基礎の加速度）も測定している場合は、センサ配置表でその ch の「入力」にチェックを入れます
            （ch 名称に GL・地盤・地表を含む場合は自動でチェック）。解析を実行すると入出力同定（EMA）が追加で行われ、
            「入出力（FRF）」タブに周波数応答関数・安定化線図・刺激関数 βφ が表示されます。
            安定化線図を Shift+クリックすると極を選び直せます（その場で再計算）。
            動作確認にはツールバーの「サンプル（地動入力あり）」を使えます。
        11. 
          
      11. 「🖨 計算書出力」で A4 縦の計算書を新しいウィンドウに開き、ブラウザの印刷機能で PDF 保存します。
        12. 
          
      12. 入力一式は「💾 保存」で JSON に保存、「📂 開く」で読み込めます。
            結果の数値は「📤 CSV出力」で書き出せます。
        13. 
        
      
        
      
          動作を確認したい場合は、ツールバーの「サンプル読込」で
          既知の固有振動数・減衰定数を持つ合成データを読み込めます。
          解析結果がその真値に一致するかで、解析ロジックの妥当性を確認できます。
        
      

        
      入力のポイント
      
        
        
          
      • 単位は結果に影響しません。同定される固有振動数・減衰定数・（正規化された）モード形状は、
            入力データのスケールに依存しません。単位は時刻歴グラフと計算書の表記にのみ使われます。
            したがって生の count 値のままでも同定できます。
        • 
          
      • 記録長は減衰定数の精度を左右します。目安として、1 次固有周期の
            1000 周期以上（例: 1 Hz の建物なら 1000 秒 ≒ 17 分）が推奨されます。
            10 分程度では減衰定数に ±20〜30% のばらつきが生じます。
        • 
          
      • 2 ch 以上が必須です。FDD は第 1／第 2 特異値の比を用いるため、1 ch では動作しません。
        • 
          
      • 参照 ch（センサ配置表の「参照」列）は SSI で参照型（SSI-ref）として使う列です。
            応答が大きく S/N の良い ch（通常は最上階）を選ぶと計算が軽くなり精度も上がります。
            未選択の場合は全 ch を参照として扱います。
        • 
          
      • ねじれ成分を分離するには、同一階に平面位置の異なる同方向センサが 2 点以上必要です。
            センサ配置表の「平面 X / Y」に設置位置を入力してください。
        • 
        
      

        
      適用範囲・注意事項（先に確認）
      
        
        
          
      • 線形・時不変・比例減衰に近い系を前提とします。大振幅時の非線形挙動や、
            測定中に剛性が変化する場合は適用外です。
        • 
          
      • 出力のみの同定（FDD / EFDD / FSDD / SSI）は、入力が定常な広帯域（白色雑音に近い）外力であることを仮定します。
            特定周波数の機械振動（空調・ポンプ等）が卓越している場合、それを構造モードと誤認する危険があります。
            卓越方向・モード形状の形（節の数）と MPC/MPD で妥当性を判断してください。
            地動を測定している場合の入出力同定（EMA）はこの仮定を必要とせず、地震・強風など非白色な外力にも使えます
            （ただし地動以外の測定していない外力が大きいとコヒーレンスが下がり、精度が落ちます）。
        • 
          
      • 実装した手法は FDD / EFDD / FSDD と SSI-cov（cov, cov_R）です（データ診断図は v260917_r1 で追加）。
            pLSCF・クラスタリングによる自動モード抽出・不確かさ伝播・
            多セットアップ合成は実装していません。
            pLSCF は v260917_r1 で一度追加しましたが、減衰定数が系統的に偏る（コレログラム法では指数窓の影響が補正されず過大、
            ピリオドグラム法では過小）ことを確認したため、v260917_r2 で削除しました。
            また SSI のデータ駆動型（method = 'dat'）も未実装です
            （N×(r·q+l·(p+1)) の巨大行列の QR 分解が必要で、10⁵ 点規模ではブラウザで実行できないため）。
        • 
          
      • 減衰定数は同定手法・記録長・雑音レベルに敏感です。同一データでも手法間で
            30〜50% 異なることは珍しくありません。設計・診断に用いる場合は必ず複数手法で照合し、
            既往の統計値（日本建築学会『建築物の減衰』等）と比較してください。
        • 
          
      • モード形状図・アニメーションはセンサ配置の入力に依存します。
            各階を剛床と仮定してセンサ値から並進・ねじれ成分を最小二乗で推定しています。
            これは図示のための換算であり、同定結果そのものではありません。
        • 
          
      • 本アプリは設計者の判断を支援するツールであり、同定結果の妥当性の最終判断は
            技術者が行う必要があります。
        • 
        
      

        
      

        
      (1) 全体フロー
      
        
      解析は次の順序で実行されます。
      
        
      
          
      段階処理
      
          
            
      ①前処理: トレンド除去 → 帯域フィルタ → 間引き
      
            
      ②クロススペクトル密度行列 Syy(f) の推定
      
            
      ③各周波数で特異値分解 → CMIF
      
            
      ④卓越振動数（対象モード）の選択（自動検出＋グラフのクリックで修正）
      
            
      ⑤FDD: ピークピッキングで fn とモード形状
      
            
      ⑥EFDD/FSDD: SDOF ベル → 自己相関 → 対数減衰率で ξ
      
            
      ⑦SSI: Hankel 行列の構成
      
            
      ⑧SSI: 特異値分解 → 可観測行列 → 状態空間行列 A, C
      
            
      ⑨SSI: 次数を増やしながら固有値解析（+ Hard Criteria）
      
            
      ⑩Hard Criteria（MPC / MPD）と Soft Criteria（安定極判定）
      
            
      ⑪SSI: 安定極から対象モードを抽出
      
            
      ⑫データ診断図（自己相関・PSD・確率密度・正規確率・STFT）、合成スペクトル
      
            
      ⑬モード品質指標（MAC / MPC / MPD / MCF）の算出と手法間比較
      
            
      ⑭（地動 ch がある場合）入出力同定: FRF の H1 推定 → 極の推定（LSCF / LSCE / RFP）→ 安定化判定 → 極の選択 → LSFD → 刺激関数 βφ
      
                
        
      

        
      (2) FDD（周波数領域分解）
      
        
      
          定常白色雑音で加振された線形系の出力クロススペクトル密度行列 Syy(f) は、共振点近傍で
          そのモードのモード形状 φ による 1 階数（rank-1）行列に近づきます。したがって Syy(f) を
          特異値分解すると、共振点では第 1 特異値が卓越し、第 1 特異ベクトルがモード形状に比例します。
        
      
        
      
          $$ S_{yy}(f) \approx \sum_{j} \frac{d_j\,\varphi_j\,\varphi_j^{\mathsf{H}}}{i2\pi f - \lambda_j}
             + \frac{\bar{d_j}\,\bar{\varphi_j}\,\bar{\varphi_j}^{\mathsf{H}}}{i2\pi f - \bar{\lambda_j}} $$
        
      
        
      [Brincker, Zhang & Andersen (2001), Smart Mater. Struct. 10(3)]
      
        
      
          本アプリでは選択振動数 fsel ± DF1 の帯域で
          第 1 特異値 / 第 2 特異値の比が最大となるビンをピークとして採用します。
          この比が大きい点は「1 つのモードだけが卓越している点」であり、モード形状の推定精度が最も高くなります。
        
      
        
      
          $$ k^{*} = \arg\max_{k\,:\,|f_k - f_{sel}| \le DF_1} \frac{\sqrt{\sigma_1(f_k)}}{\sqrt{\sigma_2(f_k)}},
             \qquad f_n = f_{k^{*}},\quad \varphi = \frac{\overline{u_1(f_{k^{*}})}}{\left|\,\overline{u_1}\,\right|_{\max}} $$
        
      
                
      FDD は減衰定数を算出しません。減衰が必要な場合は EFDD/FSDD または SSI を使用してください。
      

        
      (3) EFDD / FSDD（減衰定数の同定）
      
        
      
          FDD で得たモード形状 φ を使い、そのモードだけを含む1 自由度系スペクトル（SDOF ベル）を切り出します。
          FSDD ではモード形状で空間的に濾波した拡張 PSD を、EFDD では第 csm 特異値そのものを用います。
          どちらも MAC(φ, ucsm(f)) > MAClim を満たす周波数ビンのみを採用することで、
          他モードの混入を防ぎます。
        
      
        
      
          $$ S_{\mathrm{sdof}}(f) = \begin{cases}
             \varphi^{\mathsf{H}} S_{yy}(f)\,\varphi & (\mathrm{FSDD}) \\
             \sqrt{\sigma_{csm}(f)} & (\mathrm{EFDD})
             \end{cases} \quad \text{ただし } \mathrm{MAC}(\varphi, u_{csm}(f)) > \mathrm{MAC}_{lim} $$
        
      
        
      [Zhang, Wang & Tamura (2010) MSSP 24(5) = FSDD ／ Brincker, Ventura & Andersen (2001) IMAC XIX = EFDD]
      
        
      
          この SDOF ベルを逆 Fourier 変換すると自由振動波形に相当する自己相関関数が得られます。
          その極値（正負の山）の減衰から対数減衰率を求め、減衰定数に換算します。
        
      
        
      
          $$ \tilde{R}(\tau) = \frac{\mathrm{Re}\left[\mathcal{F}^{-1}_{\mathrm{ortho}}(S_{\mathrm{sdof}})\right]}{\;\cdot\;\Big|_{\tau=0}},
             \qquad \delta_i = \ln\frac{|\tilde{R}(\tau_0)|}{|\tilde{R}(\tau_i)|},
             \qquad \lambda_0 = \frac{\sum_i i\,\delta_i}{\sum_i i^2} $$
        
      
        
      
          $$ \lambda = 2\lambda_0 \;\;(\text{method\_SD}=\mathrm{per}), \qquad
             \lambda = 2\lambda_0 - \frac{1}{\tau_w} \;\;(\text{method\_SD}=\mathrm{cor},\;\; \tau_w = \tfrac{-(n_{bins}-1)}{\ln 0.01}) $$
        
      
        
      
          $$ \xi = \frac{\lambda}{\sqrt{4\pi^2 + \lambda^2}}, \qquad
             f_d = \left[\overline{2\Delta\tau_{\mathrm{peak}}}\right]^{-1}, \qquad
             f_n = \frac{f_d}{\sqrt{1-\xi^2}} $$
        
      
        
      [極値は先頭 sppk 個を読み飛ばし np_max 個を採用する]
      
        
      
          ※ 先頭の極値を読み飛ばす（sppk）のは、逆変換直後の区間には窓関数やゼロ詰めの影響が残るためです。
          採用する極値の個数（np_max）を増やすと平均化されて安定しますが、後半は振幅が小さく雑音に埋もれるため、
          自己相関のグラフを見て適切な範囲を選んでください。
        
      

        
      (4) SSI（確率的部分空間同定）
      
        
      
          時間領域で状態空間モデルを同定する手法です。出力データから Hankel 行列（将来出力と過去出力の
          相互共分散）を作り、その特異値分解から可観測行列 O を取り出します。
        
      
        
      
          $$ H = \frac{1}{N}\,Y^{+}\left(Y^{-}\right)^{\mathsf{T}} = U\Sigma V^{\mathsf{T}},
             \qquad O_n = U_{[:,1:n]}\,\Sigma^{1/2}_{[1:n,1:n]} $$
        
      
        
      
          可観測行列は O = [C, CA, CA², …, CAᵖ]ᵀ というシフト構造を持つため、
          上下にずらした部分行列から状態行列 A を最小二乗で求められます。
        
      
        
      
          $$ O^{\uparrow} A = O^{\downarrow} \;\;\Rightarrow\;\; A_n = \left(O^{\uparrow}\right)^{+} O^{\downarrow},
             \qquad C_n = O_{[1:l,\,1:n]} $$
        
      
        
      [Van Overschee & De Moor (1996)。高速実装は Döhler & Mevel (2013) MSSP 38(2)]
      
        
      
          $$ A_n\psi_j = \mu_j \psi_j, \quad \lambda_j = \frac{\ln\mu_j}{\Delta t},
             \quad f_{n,j} = \frac{|\lambda_j|}{2\pi}, \quad \xi_j = -\frac{\mathrm{Re}(\lambda_j)}{|\lambda_j|},
             \quad \varphi_j = C_n\psi_j $$
        
      
        
      
          実際のモデル次数は未知なので、次数 n を ord_min から ord_max まで変えて同定を繰り返します。
          真の構造モードは次数を変えても同じ位置に現れますが、数学的・雑音由来の偽の極は次数ごとに移動します。
          この性質を可視化したものが安定化線図です。
        
      

        
      (5) 極の検証基準（Hard Criteria / Soft Criteria）
      
        
      
          同定された極から物理的に意味のあるものを選別するため、2 段階の基準を適用します。
        
      
        
      Hard Criteria（HC） — 各次数で個々の極を判定
      
        
        
          
      • (a) 共役対の存在: 実行列 A の固有値は必ず共役対で現れます。対になっていない極は
            数値的な偽の極なので除去し、対のうち片方（虚部が正の側）のみを残します。
            （実装注: 複素 QR 法では丸め誤差で固有値が厳密な共役にならないため、
            共役対を許容差でグループ化して対称化した後に、共役の有無を厳密一致で判定しています）
        • 
          
      • (b) 減衰定数の範囲: 0 < ξ < ξmax。負の減衰や非現実的に大きい減衰を除きます。
        • 
          
      • (c) MPC ≧ MPClim: モード形状の実部と虚部の直線性。低い極は物理モードらしくありません。
        • 
          
      • (d) MPD ≦ MPDlim: 位相のばらつき。大きい極は物理モードらしくありません。
        • 
        
      
        
      Soft Criteria（SC） — 隣接次数間の安定性を判定
      
        
      次数 n の各極について、次数 n−1 の中で振動数が最も近い極と比較し、3 条件すべてを満たせば「安定極」とします。
      
        
      
          $$ \frac{|f_n^{(k)}-f_n^{(k-1)}|}{f_n^{(k)}} < \mathrm{err}_{fn}, \qquad
             \frac{|\xi^{(k)}-\xi^{(k-1)}|}{\xi^{(k)}} < \mathrm{err}_{\xi}, \qquad
             1-\mathrm{MAC}\!\left(\varphi^{(k)},\varphi^{(k-1)}\right) < \mathrm{err}_{\varphi} $$
        
      
        
      [Reynders, Houbrechts & De Roeck (2012) MSSP 29]
      

        
      (5b) 入出力同定（EMA）— 地動を入力として測定している場合
      
        
      
          地表・基礎で測った地動 üg を入力、各階の応答を出力として、両者の比である
          周波数応答関数（FRF）を直接求めます。出力のみの同定と違い入力スペクトルの白色性を仮定しないため、
          地震観測記録などにも適用でき、さらにモードごとの刺激関数 βφ（刺激係数 × モード形状）が得られます。
        
      
        
      ① FRF の推定（H1） — Welch 法のクロススペクトル G から、入力が複数（X・Y の地動など）の場合も含めて
      
        
      
          $$ \mathbf{h}_p(f) = \mathbf{G}_{xx}(f)^{-1}\,\mathbf{g}_{x y_p}(f), \qquad
             \gamma_p^2 = \frac{\mathrm{Re}\big(\mathbf{g}_{x y_p}^{\mathsf{H}} \mathbf{h}_p\big)}{G_{y_p y_p}} $$
        
      
        
      
          床の加速度計は地動を含む絶対加速度を測るため、同方向の入力については Hrel = Habs − 1 として
          相対加速度 / 地動加速度の FRF に換算します（「応答は絶対量」をオフにすると換算しません）。
          γ² は多重コヒーレンスで、1 に近いほど出力が測定した入力で説明されている（FRF の精度が高い）ことを示します。
        
      
        
      ② 極の推定 — FRF に共通分母の有理関数モデルを当てはめ、分母多項式の根から極 λ を求めます。次数を 1 から nmax まで上げながら繰り返します。
      
        
      
          
      手法考え方特徴
      
          
            
      LSCFz 領域（z = e−iωΔt）の共通分母モデル。正規方程式
              [Σ(T − SᵀR⁻¹S)]a = 0 を解く（R, S, T は FRF の逆 FFT から作る Toeplitz 行列）
              高次でも数値的に安定で、安定化線図が明瞭。推奨
      
            
      LSCE帯域の FRF を逆変換したインパルス応答に自己回帰モデル Σβjh(t+j) = 0 を当てはめる
              時間領域の古典的手法。帯域内のビン数に比例して計算量が増える
      
            
      RFP有理分数多項式（Forsythe の直交多項式、分子次数 20）
              狭い帯域向き。広帯域では不安定になりやすい（次数 20 まで）
      
            
      RFP 分割帯域を 800 ビンごとに分割して RFP を適用ビン数が多い場合の RFP
      
                
        
      
        
      
          $$ f_n = \frac{|\lambda|}{2\pi}\,\mathrm{sign}(\mathrm{Im}\,\lambda), \qquad \xi = -\frac{\mathrm{Re}\,\lambda}{|\lambda|} $$
        
      
        
      ③ 安定化判定 — 次数 n の各極について、次数 n−1 の極の中に振動数の相対差 < errf（既定 0.001）となる極があり、
          かつ減衰の相対差 < errξ（既定 0.05）となる極があり、ξ > 0 のとき「安定極」とします（安定化線図の緑 ×）。
      
        
      ④ 極の選択 — 既定の「自動」では、対象モードの振動数 fsel ± rtol·fsel の範囲にある安定極の
          振動数の中央値を求め、それに最も近い極を採用します。「入出力（FRF）」タブの安定化線図を Shift+クリックすると、
          画面上で最も近い安定極を個別に選べます（選択はその場で反映されます）。
      
        
      ⑤ LSFD（留数と剰余項）と刺激関数 — 選んだ極を固定し、FRF を次式で最小二乗近似して留数 A と剰余項を求めます。
      
        
      
          $$ H_{pq}(\omega) = -\omega^2 \sum_{r}\left[\frac{A_{pq,r}}{i\omega - \lambda_r} + \frac{A^{*}_{pq,r}}{i\omega - \lambda^{*}_r}\right]
             + LR_{pq} - \omega^2\,UR_{pq} $$
        
      
        
      
          地動を受ける比例減衰系では、相対加速度 / 地動加速度の FRF が
          H(ω) = Σ ω²βrφr / (ωr² − ω² + 2iξrωrω) となることから、
        
      
        
      
          $$ \beta_r \varphi_{p,r} = -2i\,\omega_{d,r}\,A_{pq,r} \quad(\text{複素留数}),\qquad
             \beta_r \varphi_{p,r} = -A_{pq,r} \quad(\text{実留数・比例減衰}),\qquad
             \beta_r = \frac{\boldsymbol{\varphi}_r^{\mathsf{T}} \mathbf{M}\mathbf{1}}{\boldsymbol{\varphi}_r^{\mathsf{T}} \mathbf{M} \boldsymbol{\varphi}_r} $$
        
      
        
      
          ※ βφ は各モードで留数が最大の入力（参照入力）に対する値です。全モードを同定できていれば各 ch の Σβφ は 1 に近づきます。
          1 次モードの屋上の βφ は、せん断型の建物で 1.2〜1.4 程度が目安です。
      
          ※ 剰余項は解析帯域より下（LR: 定数）・上（UR: −ω² に比例）のモードの影響を近似するもので、既定で両方考慮します。
      
          ※ 実装上の修正点: LSCE で下限振動数 flo > 0 とした場合、切り出した帯域のインパルス応答の時間刻みは
          Δt′ = 1/(2(fup − flo)) であり、周波数の原点も flo だけずれます。本アプリはこの刻みとずれを用いて極を求めます
          （全帯域の刻み 1/(2fup) を用いると振動数が fup/(fup − flo) 倍に偏るため）。
        
      

        
      (6) モード品質指標の読み方
      
        
      
          
      指標意味良好の目安悪い場合に疑うこと
      
          
            
      MAC2 つのモード形状の相関（0〜1）
              対角 ≧ 0.8
      非対角 ≦ 0.3      		
              非対角が大 → 同一モードの重複抽出、近接モードの混在
      
            
      MPC実部と虚部の直線性（0〜1）≧ 0.7
              S/N 不足、近接モードの混在、非比例減衰
      
            
      MPD振幅重み付き位相ばらつき [rad]≦ 0.5
              同上。特に 0.8 を超える場合は物理モードでない可能性が高い
      
            
      MCFモード複雑度（0=実モード）≦ 0.3
              同上
      
            
      σ₁/σ₂ピークでのモード分離度大きいほど良い
      （目安 3 以上）      		
              小さい → 近接モードが重なっている、雑音が卓越
      
                
        
      

        
      (6b) データ診断図の読み方
      
        
      
          
      図見方
      
          
            
      正規化自己相関
              構造の卓越振動成分があると減衰振動波形になる。ゆっくり揺れる成分が支配的なら低周波ドリフトの混入を疑い、ハイパスフィルタを検討。
      
            
      PSD（Welch 法）
              1 ch ごとのパワースペクトル。鋭すぎる（幅のない）ピークは機械振動などの調和成分の可能性。
      
            
      確率密度
              常時微動は正規分布に近いのが普通。中央が尖り裾が重い（過剰尖度 > 0）なら衝撃的外乱、
                双峰形・平坦（過剰尖度 < 0）なら周期成分の卓越を疑う。
      
            
      正規確率プロット
              点列が直線 y = x に沿えば正規分布。両端が外れる場合は外れ値（外乱）がある。
                縦軸の基準分位点は、シードを 0 に固定した標準正規乱数 n 個を昇順に並べた値。
      
            
      STFT スペクトログラム
              構造モードは時間を通して連続する横縞。時間的に局在した縦縞・斑点は外乱。該当区間を除いた再解析を検討。
      
            
      振動数－減衰図
              安定極が (f, ξ) 平面で塊になっている位置が真のモード。塊の広がりは同定のばらつきの目安。
      
            
      合成スペクトル比較
              SSI で同定したモデルから合成したスペクトルと実測を比較。ピークの位置・形がよく合う次数ほどモデルが妥当。
      
            
      Hankel 特異値 √σ
              値が急に小さくなる位置がシステムの実効次数。ordmax はそれより大きく取る。
      
            
      MAC ヒートマップ・モード複素性
              MAC は (6) 参照。複素性図は矢印が 0°/180° の線上に並ぶほど実モード（良好）。
      
            
      時刻歴の RMS 線
              arms = √(Σx²/N)。ch 間の振幅レベル比較に使う。
      
                
        
      
        
      ※ 確率密度・正規確率プロットの参考線（標準正規分布・y = x）と、歪度・尖度による「要確認」の目安は、読み取りやすくするために追加した補助表示です。
      

        
      (7) 結果の妥当性確認（チェックリスト）
      
        
        
          
      1. 手法間で振動数が一致しているか。FDD／EFDD／SSI の固有振動数の差は通常数 % 以内です。
            「同定結果」タブの手法間比較表で「f の変動」を確認してください。
        2. 
          
      2. 手法間のクロス MAC が 0.8 以上か。下回る場合はモードの取り違えが疑われます。
        3. 
          
      3. モード形状の形が構造的に妥当か。1 次は節なし（片持ち梁状）、2 次は節 1 個、
            3 次は節 2 個…という順序になっているのが自然です。「モード形状」「アニメーション」タブで確認します。
        4. 
          
      4. 略算固有周期との比が妥当か。T = h(0.02 + 0.01α) との比は、常時微動レベルでは
            0.4〜0.8 程度になるのが一般的です（実剛性が設計想定より高いため）。
            1 を大きく超える場合は剛性低下（損傷）を検討します。
        5. 
          
      5. 減衰定数が常識的な範囲か。常時微動レベルの 1 次減衰定数の目安は
            RC 造 1〜3%、S 造 0.5〜2%、木造 2〜5% 程度です。大きく外れる場合は記録長不足・
            近接モードの混在・機械振動の誤認を疑ってください。
        6. 
          
      6. MPC / MPD が基準内か。「MAC 行列」タブの「モード品質指標」で「要確認」となったモードは、
            結果の採用前に必ず個別に検討してください。
        7. 
          
      7. 安定化線図で安定極が縦一列に揃っているか。揃っていない周波数を対象モードに
            入れている場合、そのモードは信頼できません。
        8. 
        
      

        
      (8) よくある失敗と対処
      
        
      
          
      症状原因対処
      
          
            
      「極値が N 個しか見つかりません」エラー
              SDOF ベルが狭すぎて自己相関の振動が足りない
              np_max を小さくする／DF2 を大きくする／MAC_lim を下げる
      
            
      「MAC > MAC_lim を満たすビンがありません」エラー
              MAC_lim が厳しすぎる、選択振動数がピークから外れている
              MAC_lim を 0.85 程度まで下げる／スペクトルを見て選択振動数を修正
      
            
      「DF1 が周波数分解能に対して小さすぎます」エラー
              Δf = fs/nxseg が粗い
              n_xseg を大きくする（2048 等）／DF1 を大きくする
      
            
      find_min で次数が見つからない
              全モードで同時に安定極が揃う次数が無い
              rtol を大きくする／ord_max を大きくする／SC を緩める／
                対象モードを減らす／「固定次数を指定」に切替えて安定化線図から次数を選ぶ
      
            
      減衰定数が異常に大きい／小さい
              記録長不足、近接モードの混在、機械振動の誤認
              記録長を確認（1 次周期の 1000 周期以上）／DF2 を狭めて近接モードを分離／
                機械振動の可能性がある周波数を除外
      
            
      解析に時間がかかりすぎる
              データ点数・ch 数・br・ord_max が大きい
              前処理の間引き q を上げる（例 q=2）／br を小さくする／
                step を 2 以上にする／使用 ch を絞る／参照 ch を指定する
      
            
      モード形状図が高さ方向に表現されない
              センサ配置の「高さ GL+」が未入力（すべて同じ）
              センサ配置表の高さを入力（「階高を自動入力」ボタンも利用可）
      
                
        
      

        
      (9) 数値計算の実装について
      
        
      
          本アプリは外部ライブラリの追加を避けるため、必要な数値計算をすべて自前実装しています
          （数式表示の KaTeX のみ CDN から読込）。下表の中欄は、同じ処理を行う代表的な数値計算ライブラリ（Python の NumPy / SciPy）の関数です。
        
      
        
      
          
      機能NumPy / SciPy本アプリの実装
      
          
            
      FFTnp.fft.*
              基数 2 Cooley-Tukey（2 の冪乗長）＋ Bluestein chirp-z（任意長）
      
            
      クロススペクトルscipy.signal.csd
              Welch 法を再現（Hann 窓・detrend='constant'・density スケーリング・片側 2 倍）
      
            
      特異値分解np.linalg.svd
              片側 Jacobi 法（実・複素の両方。特異値は降順）
      
            
      QR 分解np.linalg.qr
              Householder 変換（reduced モード）
      
            
      一般固有値問題scipy.linalg.eig
              Householder による Hessenberg 化 → 複素 Wilkinson シフト QR 法 →
                Schur 形の後退代入で固有ベクトル
      
            
      IIR フィルタ設計butter, cheby1
              アナログプロトタイプ → 周波数変換（lp2lp/hp/bp/bs）→ 双一次変換 → SOS
      
            
      零位相フィルタsosfiltfilt
              odd 拡張パディング ＋ 定常初期条件（sosfilt_zi 相当）＋ 前進後退
      
            
      間引きdecimate
              cheby1(8, 0.05, 0.8/q) ＋ 零位相 ＋ q 点間隔抽出
      
            
      極値検出find_peaks
              プラトー中央を採用する SciPy と同一のアルゴリズム
      
            
      曲線フィットcurve_fit
              δ = λ·i の場合は閉形式 λ = Σ(i·δ)/Σi²（厳密解）
      
            
      連立方程式・逆行列np.linalg.solve, inv
              部分ピボット付き LU 分解（実数・複素数）
      
            
      乱数np.random.seed, randn
              MT19937 と旧 API の極座標法ガウス乱数を再現（NumPy と同一の数列）
      
            
      補間interp1d
              linear・外挿あり（searchsorted と端点クリップを SciPy と同じ手順で）
      
            
      PSD・STFTwelch, stft
              Hann 窓・片側。stft は両端 nperseg/2 のゼロ詰めと padded=True を再現
      
            
      多項式の根np.roots
              コンパニオン行列を平衡化（対角スケーリング）し、複素シフト QR 法で固有値を求める。実係数では共役対を厳密に揃える
      
            
      最小二乗（最小ノルム解）np.linalg.lstsq
              Householder QR で直交変換した後、R の特異値分解で解く（rcond = eps·max(m, n)）
      
                
        
      
        
      
          実装上の工夫・留意点
      
          ① 共役対の判定: 上記 (5)(a) の通り、複素 QR 法の結果を対称化してから厳密一致判定を適用しています。
      
          ② A の算出: QR 分解後の inv(R)·S を、数学的に等価で数値的により安定な
            三角方程式の後退代入で解いています。
      
          ③ Hankel 行列（cov）の計算: 反対角線上でブロックの窓が 1 サンプルずつずれる性質を使い、
            漸化式で更新しています（結果は厳密に同一、計算量は約 br/2 分の 1）。
      
          ④ SSI の 'dat' 法・不確かさ伝播・クラスタリング・多セットアップは未実装です（上記「適用範囲」参照）。
      
          ⑥ 検証: FDD / EFDD / FSDD / SSI / データ診断の各計算は、同一データを Python の数値計算ライブラリで解析した結果と照合し、
            一致を確認しています（v260917_r1 時点）。入出力同定（LSCF / LSCE / RFP / RFP 分割・安定化判定・LSFD）も同様に照合しています（v260918_r2）。
      
          ⑤ Web Worker を使いません。file:// での動作を保証するため、解析は
            async 関数と setTimeout による協調的分割実行としています（進捗表示・中止が可能）。
        
      

        
      (参考文献)
      
        
        
          
      1. Brincker, R., Zhang, L., Andersen, P.
            『Modal identification of output-only systems using frequency domain decomposition』
            Smart Materials and Structures, Vol.10, No.3, pp.441-445, 2001.
        2. 
          
      2. Brincker, R., Ventura, C. E., Andersen, P.
            『Damping estimation by frequency domain decomposition』
            Proceedings of IMAC XIX, pp.698-703, 2001.
        3. 
          
      3. Zhang, L., Wang, T., Tamura, Y.
            『A frequency-spatial domain decomposition (FSDD) method for operational modal analysis』
            Mechanical Systems and Signal Processing, Vol.24, No.5, pp.1227-1239, 2010.
        4. 
          
      4. Van Overschee, P., De Moor, B.
            『Subspace Identification for Linear Systems: Theory, Implementation, Applications』
            Kluwer Academic Publishers, 1996.
        5. 
          
      5. Peeters, B., De Roeck, G.
            『Reference-based stochastic subspace identification for output-only modal analysis』
            Mechanical Systems and Signal Processing, Vol.13, No.6, pp.855-878, 1999.
        6. 
          
      6. Döhler, M., Mevel, L.
            『Efficient multi-order uncertainty computation for stochastic subspace identification』
            Mechanical Systems and Signal Processing, Vol.38, No.2, pp.346-366, 2013.
        7. 
          
      7. Reynders, E., Houbrechts, J., De Roeck, G.
            『Fully automated (operational) modal analysis』
            Mechanical Systems and Signal Processing, Vol.29, pp.228-250, 2012.
        8. 
          
      8. Allemang, R. J.
            『The Modal Assurance Criterion — Twenty Years of Use and Abuse』
            Sound and Vibration, Vol.37, No.8, pp.14-21, 2003.
        9. 
          
      9. Guillaume, P., Verboven, P., Vanlanduit, S., Van der Auweraer, H., Peeters, B.
            『A poly-reference implementation of the least-squares complex frequency-domain estimator』 Proceedings of IMAC XXI, 2003.
        10. 
          
      10. Verboven, P. 『Frequency-domain system identification for modal analysis』 Ph.D. thesis, Vrije Universiteit Brussel, 2002.
        11. 
          
      11. Brown, D. L., Allemang, R. J., Zimmerman, R., Mergeay, M. 『Parameter estimation techniques for modal analysis』 SAE Technical Paper 790221, 1979.
        12. 
          
      12. Richardson, M. H., Formenti, D. L. 『Parameter estimation from frequency response measurements using rational fraction polynomials』 Proceedings of IMAC I, 1982.
        13. 
          
      13. Maia, N. M. M., Silva, J. M. M. 『Theoretical and Experimental Modal Analysis』 Research Studies Press, 1997.
        14. 
          
      14. 日本建築学会『建築物の減衰』 丸善出版, 2000.
        15. 
          
      15. 日本建築学会『建築物の振動に関する居住性能評価指針・同解説』 第2版, 丸善出版, 2004.
        16. 
          
      16. 建築基準法施行令 第88条（地震力）、平成12年建設省告示第1461号
            — 設計用一次固有周期の略算式 T = h(0.02 + 0.01α)
        17. 
        
      

      
      
    
      
  
      

  
      
    建築物の振動測定結果解析プログラム（Building Vibration Analysis System） v260918_r2
      
    本プログラムの同定結果は工学的判断を伴う解析結果です。設計・診断への適用に際しては、複数手法での照合と技術者による検証を行ってください。
  
      


      












      
  
      
    
      解析中…
      
    
      
    
      0%