計測自動制御学会東北支部 第 281 回研究集会 (2013.6.21)
資料番号 281-11

ガイド波の速度分散特性に基づく多孔質丸棒の密度推定
Estimation of porous rod density based on ultrasonic guided
wave s dispersion○斉藤匠 ∗ ，齊藤玄敏 ∗○ Takumi Saito∗ , Hiroyuki Saito∗*弘前大学大学院*Hirosaki Univ.
キーワード : ガイド波 (Guided wave)，アレイ信号処理 (Array signal processing)，速度分散(Dispersion)，多孔質体 (Porous solid)，センブランス解析 (Semblance analysis)，最尤推定法 (Maximumlikelihood estimation)

連絡先 : 〒 036-8561

青森県弘前市文京町 3 番地 弘前大学大学院 理工学研究科 齊藤研究室

齊藤玄敏，Tel./Fax.: (0172)-39-3694，E-mail: saitoh@cc.hirosaki-u.ac.jp

1.

諸言捉えることが困難であり，また，実体波 (P 波，近年，平均寿命が延びたことや生活習慣の変化により，骨粗鬆症の患者は増加傾向にある 1) ．
骨粗鬆症は日常生活程度の負荷によって骨折を引き起こし，高齢者が寝たきりとなる原因の一つとなっている 2) ．骨粗鬆症を予防するためには，まめに検査を行い，症状の進行度を早期に把握する必要がある．しかし，現在主流となっている X 線を用いた検査法は放射線被曝を受けることから，患者のストレスの増加や気軽に何度も検査を受けることができないことが問題点として挙げられる 2) ．それに対し，超音波骨密度計測法は放射線被曝の心配は無く，軽量，安価で操作が容易である 3) ．従来の超音波骨密度計測法は踵骨に超音波パルスを透過させ，骨内伝播速度 (SOS) を計測し，それを指標に骨密度を計算する．しかし，従来の計測法では取得する波形が 1 つのみで，波の立ち上がりを正確にS 波) は高密度の部分を選択的に伝播するため
4, 5) ，骨密度の低下を正しく捉えることができず，推定速度に大きな誤差が生じてしまう．
本研究室では，ガイド波を用いて骨の実体波速度を推定し，骨密度を評価する新たな計測方法の開発に取り組んでいる．ガイド波は物理的境界によって形成された導波路に沿って長手方向に減衰することなく伝播するという特性を持っている 6) ．また，境界から 1 波長程度の深さの力学的性質を反映しているので，骨密度計測において，ガイド波を使用することにより骨の内部に存在する海綿骨の特性を評価することができる可能性がある．しかし，ガイド波には速度分散性があり，これを扱うには速度分散性を正しく理解しなければ，正確に計測することができない．7) また，複数の波形を取得するためにアレイ計測を用い，センサを展開する長さを確保するため計測部位は脛骨を想定している．

–1–

鎌田，齊藤はガイド波を利用した実体波速度

Table 1 試料物性値．Physical property.
specimen density[kg/m3 ] porosity[％]推定法を構築し，骨を均一な中実丸棒と仮定した場合の実体波速度の推定を行い，本手法の有効性を確認した 8) ．しかし，実際の骨は内部にNo.1
No.2No.3

海綿骨という多孔構造をもつ部分が存在するの800
640
48038.550.863.1で，この影響を明らかにする必要がある．
本研究では密度 (空隙率) の異なる同一素材の多孔質丸棒を 3 種類用意し，ガイド波を伝播させ，密度による速度分散性の違いを調べた．またTable 2 実体波速度参照値．Reference valueof body wave velocity.
specimen Cp [m/s] Cs [m/s]No.1
No.2No.3

速度分散特性から実体波速度の推定を行い，実体波速度から密度を推定した．これを実際の密度の値と比較し，推定結果の正確さを調べ，本1809
16481529

手法の多孔質構造に対する有効性を調べた．

E = E02.

多孔質体の理論計算式と参照値

G = G0の算出√密度，弾性係数から多孔質体の力学的性質を求める計算式 9) を用いて P 波・S 波速度を求める

Cp =√必要がある．しかし本研究で試料として用いた

Cs =ソリッドリジッドポリウレタンフォーム (Solid1−P0 +5E0 )(9K0 −E0 )P
1 + (11K4K
0 (27K0 +5E0 )1−P0 +2G0 )P
1 + 6(K
9K0 +8G0(4)
(5)

1−P0P
1 + 3K
4G0(6)

K + 34 G
ρ

(7)

G
ρ

(8)

K = K0

評価には P 波・S 波速度を用いるため，材料の811
785760

Rigid Polyurethane Form : SRPUF) は正確な

ここに，弾性係数が分かっていないため，P 波・S 波速ρ0 : 中実材料の密度, ρ : 多孔質体の密度,度を求めることができなかった．そこで，密度P : 空隙率, E0 : 中実材料の縦弾性係数,

800[kg/m3 ](試料 No.1)，640[kg/m3 ](試料 No.2)，

G0 : 中実材料の剛性率,

480[kg/m3 ](試料 No.3) の 3 種類の試料から各 6

K0 : 中実材料の体積弾性率,本の長さの異なる角材を作成し，パルス透過法

E : 多孔質体の縦弾性係数,による弾性波計測実験によって取得した波形か

G : 多孔質体の剛性率,

ら P 波速度を推定した．各 P 波速度と密度の値

K : 多孔質体の体積弾性率,から，以下に示す多孔質体の力学的性質を求め

Cp : P 波速度, Cs : S 波速度,る計算式を用いてα=0.34，E0 =3.2[GPa] と推である．

定し，これらから S 波速度の参照値を算出した．
結果を Table 2 に示す．

G0 = αE0
G0 E 0K0 =3(3G0 − E0 )
ρP = 1−ρ03.

実験

3.1

実験方法

(1)
(2)
(3)

試料には，直径 30[mm]，長さ 400[mm]，空隙率 38.5 ％ (No.1)，50.8 ％ (No.2)，63.1 ％ (No.3)の 3 種類の SRPUF 丸棒を用いた．パルス幅–2–

Fig. 1弾性係数推定結果．

Fig. 3ガイド波計測波形 (空隙率 50.8 ％)．

関が最大になるときの速度を求める方法である．
)2(∑τ +pzk +T ∑N −1

k=0 x(t, zk )t=τ +pz

S=N∑τ +pzk +T ∑N −1(9)2k=0 x(t, zk )t=τ +pz

ここに，Fig. 2実験の様子．

S : センブランス値,τ : 計算開始時刻,

10[µs] の方形パルスをトランスデューサで電気信号を超音波に変換してから試料表面に伝播させ，発振点から 100[mm] 離れた位置で信号を取得した．さらに受信信号を取得するセンサの位置を 10[mm] ずつずらして，各試料につき 6 本の波形を取得した．Fig. 3 に実験で取得した波

p = 1/v : スロウネス, zk : 伝播距離, t : 時刻,
k : 波形番号, T : ウィンドウ幅, N : 波形数,x(t, zk ) : 信号データ,である．フィルタバンクにはガウシアンフィルタを用いた．ガウシアンフィルタは計測波形に対してフーリエ変換を行い，計測波形のスペクトルに対して，ガウス関数 Hn (ω) を窓関数とし形を示す．

てかけ，得られたスペクトルを逆フーリエ変換3.2することで，設定した中心周波数付近の成分を解析方法と実体波速度推定方法抽出した波形を得るというものである。

本研究では速度解析に周波数センブランス法{と最尤推定法を用いた．

ω − ωnHn (ω) = exp −αωn周波数センブランス法は，フィルタバンクを用い，周波数ごとの波形データ群に対して，位相()2 }(10)(for ωl,n ≤ ω ≤ ωu,n )

速度をパラメータとして相関を取る方法で，加算後の平均したエネルギーに対する加算前の各波形の平均エネルギー比を求めることで，周波数ごとの波形群に対して同時に相関をとり，相Hn (ω) = 0(for ω < ωl,n or ω > ωu,n )–3–(11)α = β/BAN D2(12)ωu,n = (1 + BAN D)ωn(13)ωl,n = (1 − BAN D)ωn(14)ここに，Hn (ω) : ガウス関数, ωn : 中心周波数,α, β : 減衰率を表すパラメータ,
BAN D : バンド幅,ωu,n : 上限周波数, ωl,n : 下限周波数,である．
最尤推定法は，式 (15) で示される尤度 δ が最大になるような速度を推定する方法で，尤度δ をFig. 4 周波数センブランス解析結果 (空隙率50.8 ％)．

角周波数ω ，計算開始時刻τ ，スロウネス p の関これらの解析から得られた速度分散と，実体数として評価する (+ はエルミート行列を表す)．

波速度をパラメータとして算出した理論曲線を

δ(ω, τ, p) = [e+ (ω, p)R−1 (ω, τ, p)e(ω, p)]−1 (15)

最小二乗法で誤差評価し，最も残差が少ない速

ここに，度を実体波速度とする．さらに実体波速度から

x1 (ω, τ, p)..
X(ω, τ, p) = .


多孔質弾性体の理論式を用いて密度を推定する．

(16)

3.3

xN (ω, τ, p)R(ω, τ, p) = X(ω, τ, p)X + (ω, τ, p) (17)
を Fig. 4-8 に示す．周波数センブランス解析で(18)

ejωpzNxk (ω, τ, p) = e−iωzk p∫ τ +Tτ

は，ウィンドウ幅 T を 1024[point](=102.4[µs])，
速度分解能を 10[m/s]，時間分解能を 0.1[µs]，計x(τ + zk p, zx )W (t − τ )e−iωt dt

算開始時刻 τ を 1500 に設定し，周波数 0-200[kHz](19)の範囲を 1[kHz] 刻みで解析を行った．また，フィ(20)

ルタバンクに組み込まれたガウシアンフィルタtπ
W (t) = sin2 ( )Tの係数を BAND=0.2，β=5 と設定した．そのδ : 尤度, W (t) : 窓関数, xk (ω, τ, p) : 信号データを短時間窓フーリエ変換したスペクトル,X(ω, τ, p) : ベクトル行列,R(ω, τ, p) : 共分散行列, e(ω, p) : 単位行列,t : 時刻, T : ウィンドウ幅, ω : 角周波数,
τ : 窓のずらし幅, p = 1/v : スロウネス,k : 波形番号, x(t, zk ) : 信号データ,zk : 伝播距離, N : 波形数,である．

周波数センブランス解析と最尤推定法の結果

ejωpz1..
e(ω, p) = .


速度解析の結果と実体波速度推定

結果，周波数 10〜20[kHz] 付近でたわみ波 1 次モード，20〜80[kHz] 付近でたわみ波 2 次モードと考えられるピークが抽出できた．
最尤推定法ではウィンドウ幅 T を 1024[point](=102.4[µs])，速度分解能を 2[m/s]，時間分解能を0.2[µs] とし， 周波数分解能は 2.441[kHz] で解析を行った．周波数 10〜40[kHz] 付近でたわみ波 1 次モード，20〜80[kHz] 付近でたわみ波 2 次モード，80〜120[kHz] 付近で P 波と考えられる–4–

Fig. 5 P 波のピークの例 (空隙率 50.8 ％，周波数 95.199[kHz])．
Fig. 7 たわみ波 2 次モードのピークの例 (空隙率 50.8 ％，周波数 48.82[kHz])．

Fig. 6 たわみ波 1 次モードのピークの例 (空隙率 50.8 ％，周波数 24.41[kHz])．
ピークが抽出できた．そのため本研究では誤差評価時の P 波速度の値を発振周波数付近の帯域

Fig. 8％)．

最尤推定法解析結果ピーク (空隙率 50.8で検出できた P 波速度の平均値で設定し，S 波4.

速度の推定を行った．

結言

さらに，推定した S 波速度の値から密度を計

本研究では，実験から多孔質丸棒を伝播する算し，誤差を調べた．その結果，空隙率 38.5 ％，ガイド波の伝播波形を計測し，その波形を周波50.8 ％の試料では実際の値に近い S 波速度が推数センブランス法と最尤推定法により，速度分定でき，密度の誤差も Table 6 に示すように，散性を評価した．その結果から得られた実際の1.48 ％，0.675 ％と小さいものだったが，空隙密度と解析値の誤差から，空隙率 50 ％付近まで率 63.1 ％の試料では 83.20 ％もの誤差が生じるはこの方法で密度の推定が可能であるが，今回ことが分かった．

用いた理論計算式では 50〜60 ％を境に大きな誤差が生じることが分かった．

–5–

Table 3 P 波ピークの平均値．Average of P
wave velocity peak.
Cp [m/s]
specimenRef average
No.1
No.2
No.3

1809
1648
1529

1781
1648
1509

811
785
760

812
757
686

810
787
709

810
786
705

specimen

Ref

MLM

Sem

ALL

No.1
No.2
No.3

800
640
480

806.1
455.5
-127.6

794.6
655.8
80.6

794.6
649.5
45.8

Table 6
ror.

No.1
No.2
No.3

9) 近藤連一: 多孔材料-性質と利用-，272/276，技報堂出版 (1973)
10) 森田真史，山本真，内田貞夫: 海綿骨の機械的性質とその骨梁構造依存性について (1984)

density error[％]
MLM
Sem ALL
0.763
28.80
1270.675
2.4783.20

6) 工藤圭，廣瀬壮一: SH 板波を用いた平板の減
肉欠損の逆解析，土木学会第 64 回年次学術講演会，I-127（2009）8) 鎌田正洋，齊藤玄敏: フィルタバンク法による骨を伝播するガイド波の速度分散特性の評価，
計測自動制御学会東北支部第 265 回研究集会，
265-10(2011)

実際の密度に対する誤差．Density er-specimen5) 松川真美，水野勝紀，長谷芳樹，細川篤: 超音波による定量的な骨評価に向けて-海綿骨中の縦波超音波伝搬挙動の解明-，Fundamentals Review，Vol.3，No.4，pp.47-52(2010)7) 大谷靖弘: 配管壁貫通部の腐食診断技術，JFE
技報 No.25，41/46，(2010)

密度計算結果．Calculated density.
density[kg/m3 ]Table 53) 大谷隆彦: QUS の標準化における基本的な考え
方-工学側面からみた標準化とその問題点，第10 回日本骨粗鬆症学会イブニングワークショッ
プ，Osteoporosis Japan，17-2(2009)4) 長岡洋樹，須崎琢而，大谷隆彦，真野功，堀井薫:
超音波による骨強度の定量測定，REABOUT，
29(2004)

Table 4 S 波速度推定結果．Estimated S wavevelocity.
Cs [m/s]
specimen Ref MLM Sem ALL
No.1
No.2
No.3

2) 小野田敏行: DEXA 法による骨密度の集団検診に関する基礎的研究

0.675
1.48
90.50

一般的な海綿骨の空隙率は 50 ％以上であり
10) ，骨粗鬆症の場合は空隙の割合がより大きくなるため，今回用いた理論計算式で骨密度を計測することは難しい．今後は，高い空隙率にも対応する理論計算式について調査する必要がある．

参考文献
1) Dr Ambrish MIthal，Dr Vibha Dhingra，DrEdith Lau: The Asian Audit Epidemiology，costs and burden of osteoporosis in Asia
2009，International Osteoporosis Foundation，
30/32(2009)–6–