B-4 2 SSCQE法

通信システム工学講座
平成17年度修論・卒論
B-4
MPEG-2 符号化動画像の No Reference 画質評価モデル
川除 佳和 (指導教授 村井 忠邦)
1
まえがき
3.1
近年, 高速 IP サービスやマルチメディア家電の普
及によって, デジタル動画像を利用する機会が増えて
きた. 動画像の符号化方式として, H.264 や MPEG2 等が挙げられるが, 特に MPEG-2 は DVD や BS・
CS, 地上デジタル放送の符号化に用いられ, 非常に重
要な役割を担っている. そのデジタル動画像を利用
する環境において, 動画像の品質の善し悪しは, サー
ビスの質を大きく左右する重要な要素である. 動画
像の品質を評価する方法として, 人間を用いた主観
評価法が用いられるが, 膨大な時間・労力や専用の
設備が必要となるため, 非現実的である. そのため,
符号化動画像の客観的な物理量 (符号化歪み等) か
ら主観品質を推定する技術 (客観評価法) の開発が望
まれている. 本研究では, 主観評価法の一つである,
SSCQE 法に対応し, かつ原動画像情報を一切使用し
ない No-Reference(NR) 型の客観評価法を提案する.
NR 型の客観評価法は, VQEG と呼ばれる標準化団
体が RRNR-TV テストプランとして議論を行って
おり, 提案する客観評価法は, このテストプランの要
求用件に従って作成した.
2
SSCQE 法
画像符号化により生じる画質劣化は, 動画像の内
容に強く依存するため, 多数の動画像について主観
評価を実施する必要がある. そこで, 長時間連続的に
評価が可能な主観評価法である SSCQE 法が提案さ
れた. SSCQE 法は, スライダ機構を持つ評価デバイ
スを用いて, 動画像品質の連続評価を行うことがで
きる (得られる評価値は, 0.5 秒毎である).
3
提案する客観評価法
提案法の処理の流れを図 1 に示す. 各特徴量はフ
レーム単位で求め, SSCQE 法のサンプリング間隔に
相当する 0.5 秒毎に処理を行う.
フレーム品質推定
MPEG-2 は DCT ベースの画像符号化であるた
め, フレームで単位みれば画質の劣化が JPEG 符号
化画像と似ている. JPEG 符号化画像は, 高圧縮に
なるにつれ高周波成分が削られていくという劣化を
伴うが, この符号化歪みに着目し, 高周波成分の有無
を輝度差の頻度情報から得ることで主観評価値を推
定しようと試みた. 輝度差の頻度情報を画像の特徴
量とし, 原フレームと, エッジ保存型フィルタ処理を
したフレームについて同じように求めた. これらと
主観評価値とで重回帰分析を行い, 回帰係数を求め,
フレーム品質推定モジュールの作成を行った.
3.2
フレーム品質補正
フレーム品質は, あらかじめ作成されたデータベー
スから補正値を得ることで行う. データベースには,
SI, TI ,FQ のそれぞれの平均値と標準偏差 (これら
6 つをクリティカリティ情報と定義) と, 補正値 α(主
観評価値とフレーム品質との差, 以下単に α と表記)
が関連付けられて格納されている. データベースは,
主観評価実験を行った動画像からクリティカリティ
情報と α を求め, これら二つを対にし, K-mean 法を
用いてクラスタ化され格納される. 補正を行う際に
は, 評価対象のクリティカリティ情報をデータベー
スに入力し, 得られた α をフレーム品質と加算する
ことで補正する.
3.3
3.3.1
時間方向へのフィルタ処理
Improved Assymmetric Tracking
これは, 品質が向上する場合とそれと同程度に品
質が低下する場合とでは, 低下する方をより強く評
価値に影響させる処理である. 同時に, フレーム品質
補正結果にその直前の結果の影響を考慮させる. 次
式により表される.
8
α
>
<
α[1 − (1 − gi ){φ·V QM (n−1)+ι} ]
gs =
> −β[1 − (1 + gi ){ψ·V QM (n−1)+κ} ]
:
−β
gi
V QM (n)
=
=
gi > 1
1 ≥ gi ≥ 0
0 > gi ≥ −1
gi < −1
preV QM (n) − V QM (n − 1) (1)
V QM (n − 1) + gs
(2)
ここで，ι, κ は曲線の傾きを決めるパラメータであ
り， φ, ψ は直前の品質予測値が曲線の傾きに与える
影響を表すパラメータである. また, V QM (n) は n
サンプル目における客観評価値, preV QM (n) は n
サンプル目におけるフレーム品質補正結果である.
3.3.2
平均値フィルタ
このフィルタは TAP 数が 4 の平均値フィルタで
ある. 前段の Improved Assymmetric Tracking で
取り除ききれなかった高周波成分をカットする目的
で使用する. 次式にフィルタの仕様を示す.
V QM (n) =
n
X
i=n−3
図 1: 提案する客観評価法の処理の流れ
(B-4)
V QM (i) × 0.2
(3)
通信システム工学講座
平成17年度修論・卒論
図 3: PS1,open
図 2: PS1,closed
表 2: open データに対する推定精度 (PS1)
表 1: closed データに対する推定精度 (PS1)
rate
3.0Mbps
4.0Mbps
5.0Mbps
8.0Mbps
4
最大誤差
0.12
0.10
0.08
0.07
平均誤差
0.02
0.02
0.02
0.01
rate
3.0Mbps
4.0Mbps
5.0Mbps
8.0Mbps
outlier 率
0.03
0.04
0.01
0.00
6
closed データに対する推定
フレーム品質補正のデータベースを, 全ての評価動
画像 (PS1∼PS8) から作成した. データベース作成時
に用いた評価動画像の主観評価値を推定した (図 2).
表 1 に主観評価値と提案法の推定値との, 相関係数,
outlier 率 (95 ％信頼区間からはずれた率), 最大誤
差, 平均誤差を示す. 結果, 非常に高い精度で主観評
価値を推定できた.
5
open データに対する推定
フレーム品質補正を行うデータベースに存在しな
い, 未知の動画像に対しての推定精度を調べるため
に, 主観評価実験を行った動画像 7 つのうち, 1 種類
を未知動画像とみなし, その動画像 (PS1) を除いて
データベースの作成を行った. このデータベースを
用いて, データベース作成時に除いた動画像 (PS1)
の主観評価値の推定を行った (図 3). 表 2 に主観評
価値と提案法の推定値との, 相関係数, outlier 率, 最
大誤差, 平均誤差を示す. 結果, closed データに対す
る推定精度と比較すると精度は大幅に落ちているが,
さほど大きな誤差は無かった.
最大誤差
0.33
0.27
0.30
0.24
平均誤差
0.10
0.09
0.09
0.12
outlier 率
0.44
0.44
0.44
0.74
まとめ
本論文では, SSCQE 法に対応した MPEG-2 符号
化動画像の No-Reference 画質評価モデルの構築を行
い, その推定精度を評価した. このモデルでは, 動画
像をフレーム単位で見た場合の静止画像品質をベー
スとし, その品質を画像特徴量によって補正する方法
である. 結果, 以下のことが明らかとなった.
• JPEG 符号化動画像の NR 画質評価モデルを
作成し, そのモデルを用いてフレーム品質の
推定を行った. 結果, JPEG 符号化画像で最
適化された評価モデルであるのにもかかわら
ず (MPEG-2 符号化動画像からみると, open
データに対する推定), 推定値は主観評価値と
高い相関を示した.
• フレーム品質をクリティカリティ情報と定義
した SI, TI, FQ のそれぞれの平均値, 標準偏
差を用いて補正を行った. 結果, closed データ
に対する推定においては非常に高い精度でフ
レーム品質の補正を行うことが出来た. また,
open データに対する推定においては, outlier
率が高くなったが, 大きな誤差も無く良好な補
正結果が得られた.
(B-4)
平成17年度修論・卒論
通信システム工学講座
(B-6) 準 3 次元有限要素モデルによる浅水長波流れ解析
滝澤 祐紀
1
はじめに
3
近似解と計算結果の比較
鉛直方向の速度を考慮に入れた準 3 次元モ
計算モデルとして長さ 10.0[m] の 1 次元モ
デル (レ イヤーモデル , レベルモデル ,σ 座標モ
デルを用いた. 分割数は 40 分割とした. 理想的
デルなど ) が提案されているが , 開境界の取り
な近似解は , 式 (5), 式 (6) を用いて求めた.
扱いなどに問題が残されている. 本論文では ,
固定メッシュを用いて鉛直流速を変数として
η = acosσ(t −
扱うレベルモデルを採用し , 浅水長波解析につ
いて開境界条件を提案し , 精度を検討した. ま
た, 浅水長波流れ解析によって求めた風の流れ
を利用して雷雲電荷流れ解析を行った.
u=
x
3 gσa2
x
)−
xsin2σ(t − )
c
4 c3
c
x
1 g 2 a2
x
ga
cosσ(t − ) −
cos2σ(t − )
c
c
8 c3
c
3 g 2 σa2
x
−
xsin2σ(t − )
4 c4
c
(1)
(2)
波高と流速の初期値を図２, 図３に , 近似解
2
準 3 次元有限要素モデル
と計算結果の比較結果を図４, 図５に示す.
解析領域の高さ方向を少ない層 (２∼３層)
に分け , それぞれの層については 2 次元的流れ
とする準 3 次元モデルとして解析する. 準 3 次
元モデルを用いることにより大まかな高さ方
向の速度成分を知ることができる. そこで , レ
ベルモデルの支配方程式 [1] を用いて計算を
行った.
支配方程式に対して, 平均流速の概念等を導
図 2: 波高 (初期値)
入し解析する.(図 1 参照)
図 3: 流速 (初期値)
図 1: レベルモデル
図１において,k は層番号,η は波高,h は水
深,u,v は水平方向の流速成分,w は鉛直方向の
流速成分を表している.
開境界条件としては , 境界の処理をしない (自
然境界) として計算を行った. 時間の離散化に
ついては , 前進差分法を適用した.
図 4: 波高の比較 (t=5.0[s])
平成17年度修論・卒論
通信システム工学講座
(a) 上空から見た点 (b)x 方向から見た点
電荷分布図
電荷分布図
図 5: 流速の比較 (t=5.0[s])
図 8: 電荷の移動の様子 (t=500.0[s])
4
雷雲電荷流れ解析
図６のような宝達山付近の地形を用いて北
西から風が流入すると仮定し , 領域全体の風の
流れを計算する. その結果を用いて雷雲電荷流
れ解析を行った. 結果として点電荷分布の時間
(a) 上空から見た点 (b)x 方向から見た点
電荷分布図
電荷分布図
的変化を図７∼図９に示す.
図 9: 電荷の移動の様子 (t=1100.0[s])
5
図 6: 宝達山付近の等高線図
まとめ
• 近似解と計算結果の比較については , 結果
が十分に一致しており, 提案した開境界条件
と作成したプログラムが妥当であることが分
かる.
• 雷雲電荷流れについては , 山に沿って点電
荷が上昇・下降している様子が分かる. しかし ,
雷雲を単純に風の流れで移動するとしている
が , 雷雲には , 様々な力が働いているので点電
荷間の力など 考慮に入れない場合でもどれく
らいの精度で推定出来るかに疑問が残る.
参考文献
(a) 上空から見た点 (b)x 方向から見た点
電荷分布図
電荷分布図
図 7: 電荷の移動の様子 (t=0.0[s](初期電荷分
布))
[1] T.Kodama and M.Kawahara:”Multiple Level
Finite Element Analysis for Tidal Current
Flow with Non-Reflective Open Boundary
Condition”,Structural Eng. /Earthquake Eng.
,vol.9,No.1,PP.77s-87s(1992)
通信システム工学講座
平成17年度修論・卒論
B-9
画像の特徴量を用いた符号化動画像の画質評価モデル
本田 昌志 (指導教授 村井 忠邦)
1
まえがき
デジタル動画像を利用する場合，その情報量の大
きさから非可逆符号化は必須とされており，その際
に生じる画質劣化を定量的に評価することは，伝送
レートの効率よい分配などに非常に有効である．動
画像の画質を評価する方法として，複数の人間によ
り画像品質を直接測定する主観評価法がある．しか
し，主観評価法は人間による評価であるため，時間
と労力を非常に要するという欠点がある．この問題
を解決しつつ，主観評価と同等の性能を有する評価
を行うために，客観評価法が古くから研究されてい
る．客観評価法は，主観評価時に人間が行っている
画像の物理量から主観評価値への変換処理を事前に
数値化しておくことにより，実際の評価を物理量測
定だけで行える様にする方法である．
本研究では，動画像の主観評価法の一つである SSCQE 法の評価値を推定対象として，RR 法と呼ばれ
る客観評価法を用いて画質評価モデルの構築を行う．
2
図 4: 提案法の流れ図
SSCQE 法
SSCQE 法では，30 分から 60 分の動画像をスラ
イダ機構を持つ評価デバイスを用いて，時間的に連
続して評価を行う．評価値は 0.5 秒おきにサンプリ
ングされる．
3
RR 法
客観評価法は，原動画像の取り扱い方により 3 種
類に大別される．原動画像をそのまま使う FR(Full
Reference) 法 (図 1) は，符号化動画像と原動画像そ
のものを比較できるため高精度な推定が可能である
が，原動画像の情報量が大きいため，用途が通信を
挟まないものに限られてしまう．一方で，原動画像の
情報を全く使わない NR(No reference) 法 (図 2) は，
用途を問わないが，高精度な推定が困難とされてい
る．そのため，原動画像から幾つかの情報を抽出す
る RR(Reduced Reference) 法 (図 3) は，原動画像
のリファレンス情報を極僅かに抑えることで NR 法
並の適用範囲を実現し，さらに NR 法より容易に精
度の向上を図れるバランスの取れた方法といえる．
4
4.1
図 2: NR (No Reference) 法
フレーム品質
原動画像と符号化動画像の各フレームに平坦部や
注視点を考慮したマスキングを行った後、そこから
テクスチャ特徴量やエッジ画像の特徴量を求め，その
差分値の線形重回帰式から動画像の各フレームにお
ける静止画像品質であるフレーム品質 F Q を推定す
る．さらに主観評価時の人間の反応を考慮して，過去
5 サンプルの F Q を基に連続的なフレーム品質 F Qt
を算出する．
4.2
グローバル動き量
次に，動画像中のカメラワークなどの大局的な動
きによる視覚の解像度低下を考慮するため，グロー
バル動き量 GM V を原動画像から求める．
4.3
図 1: FR (Full Reference) 法
提案する客観評価法
本モデルにおける，処理の流れを図 4 に示す．各
特徴量はフレーム単位で求め，0.5 秒ごとの平均値
を用いて処理を行う．以降，各処理の内容を簡単に
説明する．
動的劣化量
さらに，フレームを動画像として見た場合に検知
される動的劣化も考慮する．まず，人間の視覚特性
に注目して，フレーム内のエッジなどの高周波成分
にマスキングを行う．その後，フレーム内の DCT ブ
ロックサイズの輝度平均に対して，フレーム間の動
き差分を求め，その符号化誤差を動的劣化量 T Df l
とした．
4.4
フレーム品質の補正
式 (1) により，フレーム品質 F Qt をグローバル動
き量 GM V と動的劣化量 T Df l により補正し，動画
像品質 V Q を求める．
V Q = g0 − {g2 (T Df l + g3 )
(1)
√
+(1.0 − F Qt )(1.0 − ( GM V − 1.0)/g1 )}
図 3: RR (Reduced Reference) 法
(B-9)
通信システム工学講座
平成17年度修論・卒論
図 5: PS3 3Mbps (closed)
図 7: PS3 3Mbps (open)
図 6: PS8 3Mbps (closed)
図 8: PS8 3Mbps (open)
ここで，g0 ∼ g3 はパラメータである．このパラメー
タは，GA に似た最適化手法である PSO (Particle
Swarm Optimization) を用いて決定する．
4.5
伝送レート
RR 法を用いた画質評価モデルの場合，原動画像
のリファレンス情報の伝送レートは，モデルの性能
を決定する上で重要な要素である．そこで，今まで
述べた特徴量などを，想定している最大値で伝送し
た場合の情報量の総和を計算する．その内訳は以下
の通りである．
1.
2.
3.
4.
5.
フレーム品質計算用特徴量 (32bit) × 11 個
動的劣化量 T Df l (32bit)
グローバル動き量 GM V (32bit)
平坦部除外＋注視点マスク (726bit,33 × 22)
高周波部分除外マスク (726bit,33 × 22)
一般的な NTSC 方式の TV で規定されているフ
レームレート 29.97fps で計算すると，伝送レートは
合計で約 54.62Kbps となり，非常に僅かな値となっ
ている．
5
最大誤差の大きさが気になるものの，平均誤差と
Outlier 率は共に低い値を示している．ビットレート
別に見ると，Outlier 率，平均誤差共に高ビットレー
トになるほど小さくなっているが，やはり最大誤差
だけが大きくなっている．
図 5,6 に誤差の大きいものを示す．これらの評価
値が急激に変動している部分に注目すると，主観評
価値の変動と比べて，VQ の変動が遅れている点が
目に付く．これが最大誤差が大きい原因と思われる．
6
表 2: シーケンス別の推定精度 (open)
対象
PS1
PS2
PS3
PS4
PS6
PS7
PS8
ALL
クローズドデータの推定精度
一分間の動画像 (PS1∼PS4,PS6∼PS8) を 7 個作
成し，それぞれを 3/4/5/8Mbps のビットレートで
MPEG-2 符号化した．これらを SSDQE 法で評価
し，その結果を全てクローズドデータとして PSO に
よりパラメータの最適化を行った．表 1 に推定誤差
の値を示す．
表 1: 推定精度 (closed)
対象
PS1
PS2
PS3
PS4
PS6
PS7
PS8
3Mbps
4Mbps
5Mbps
8Mbps
ALL
最大誤差
0.18
0.28
0.38
0.30
0.29
0.31
0.29
0.32
0.35
0.35
0.38
0.38
平均誤差
0.060
0.037
0.053
0.059
0.054
0.066
0.068
0.068
0.059
0.052
0.047
0.057
Outlier 率
0.31
0.13
0.19
0.22
0.21
0.28
0.33
0.26
0.24
0.23
0.22
0.24
オープンデータの推定精度
本モデルの汎用性を調べるため，オープンデータ
に対して主観評価値の推定を行う必要がある．そこ
で，7 種類のシーケンスの内 6 種類で学習を行い，残
り 1 種類をオープンデータとして評価する処理を 7
種類分繰り返した．表 2 に推定誤差を示し，クロー
スドデータと同じ動画像の結果を図 7,8 に示す．
オープンデータの結果と比べて，PS6,PS4,PS3 で
やや誤差が大きくなるものの，全体的には，同等の
精度で推定出来ている．
7
最大誤差
0.19
0.28
0.38
0.31
0.30
0.30
0.29
0.38
平均誤差
0.061
0.037
0.055
0.061
0.058
0.066
0.069
0.058
Outlier 率
0.32
0.13
0.22
0.24
0.25
0.29
0.33
0.25
まとめ
本研究では，RR 法を用いて，SSCQE 法に対応
した MPEG-2 符号化動画像の客観評価モデルを構
築した．このモデルでは，テクスチャ特徴量やエッ
ジ画像の特量量から求めたフレーム品質を，グロー
バル動き量と動的劣化量で補正した値を評価値とし
ている．実際の主観評価値を推定した所，以下のこ
とが明らかになった．
• 主観評価値は高精度に推定出来たが，最大誤
差の大きさは改善の必要がある．
• オープンデータに対しても，クローズドデー
タを用いた場合と同等の推定精度を得ること
ができた．
(B-9)
平成17年度修論・卒論
通信システム工学講座
(B-10) マイクロフォンアレ イによる音源探査のシミュレーション
藪田 将貴
1
はじめに
到来電波や到来音波の波源を推定すること
は，雷放電に伴う電磁波からの雷雲位置の予
測，構造物や機器の振動部位もし くは故障箇
所の特定など 様々な応用が考えられることか
ら精度の良い推定法が望まれている．その中
で MUSIC 法が近年注目を集めている．しか
し ，到来波を平面波と見なすことが多く，音
源が信号観測部位に近い場合の報告がほとん
ど 見られず，種々の条件下でど の程度の推定
精度を有するかは十分明らかにされていない．
また MUSIC 法とは無関係に複数の音源の位
置等を推定する試みもなされているが 、有効
な音源推定法が見出せない状況である．
そこで，観測点と音源の相対位置など 様々
な条件を考慮した場合の MUSIC 法の推定精
度を数値シミュレーションにより検討した．
2
近傍音源探査のための MUSIC 法
MUSIC 法は，マイクロフォンアレイの出力
信号を操作して到来方向を推定する．本稿で
は，リニアアレ イの近傍に音源があるものと
して推定を行った．計算は図 1 の手順で行う．
入力信号ベクトルを次のように記す．
xm = [X1,m , · · · , XK,m ]T
(1)
ここで，Xk,m は k 番目のマイクから得られ
た信号に短時間フーリエ変換 (STFT) したも
のを表し，m は m 番目の周波数成分を意味し
ている．この Xk,m から空間相関行列 Rxx を
計算する．空間相関行列はエルゴ ード 性の仮
定より時間平均から求められ，周波数的に独
立なので，
Rxx m = E[xm xH
m]
(2)
となる (E[・] は期待値 (アンサンブル平均) を
示す)．次に Rxx を固有値分解する．
−1
Rxx m = Em Λm Em
(3)
ここでは，Em は固有ベクトル行列を表し ，
Rxx から計算した K 個の固有ベクトルを列ベ
クトルとする．Λm は，Rxx の固有値を対角
成分に大きいものから順に並べた対角行列で
ある．
さらに，得られた固有ベクトルを固有値の
大きい信号部分空間の基底ベクトルとなるも
の (E s ) と，その他の直交補空間の基底ベクト
ルとなるもの (E n ) に区別する．このとき，大
きい固有値の数を到来波の数とする．ロケー
ションベクトル am (r, θ)(式 (5) を用いて計算
n
との間には，
したもの) と Em
n
{a1,m , · · · , aL,m } ⊥ Em
·
al,k
図 1: 音源近傍における MUSIC 法アルゴ リズムの流れ図
原理は，アレイで得られる (空間情報を有す
る) 信号に固有値解析を実施し，比較的低周波
の成分から所望の (音源付近の) 空間特性を得
ようとするものである．
1 −jωm τl,k
1 −jωm τl,1
e
···
e
=
rl,1
rl,K
(4)
¸T
(5)
のような直交性がある．ここで，rl,k と τl,k は
l 番目の波源から k 番目のマイクまでの距離と
伝搬時間を示す．また，ωm は離散化された角
周波数を表す．これらの性質より MUSIC ス
ペクトル (式 (6) の値) がピークをとるように，
任意のロケーションベクトルと Ekn の値を用
いて式 (6) の計算を行うことにより音源位置
の推定を行う．その際，距離による影響をな
くすためにロケーションベクトルの正規化 (式
(7)) を行った．
P (r, θ, m) =
bm (r, θ) =
1
n 2
|bH
m (r, θ)Em |
am (r, θ)
||am (r, θ)||
(6)
(7)
平成17年度修論・卒論
3
3.1
通信システム工学講座
計算条件と数値例
計算条件
解析する空間は自由空間とし ，空間の構造
等による反射波は考慮しないものとした．マ
イクロフォンアレ イには図 2 の四素子の等間
隔リニアマイクロフォンアレ イを用いた．素
子間隔を様々に変えた場合の音源位置推定と，
マイクから得られる入力にノイズを考慮した
場合の分解能について考えた．音源の振幅波
形には音源１を式 (8)，音源２を式 (9) とした．
ここでノイズは平均値 0，標準偏差 σ のガウ
スノイズとしノイズの大きさを式 (10) のノイ
ズレベルＮＬで定義した．
図 2: マイクロフォンアレ イの模式図
F (t)1 = exp(−2.0t) ×
1
exp{−j2πf1 (t − τ1 )} (8)
r1
F (t)2 = exp(−1.0t) ×
1
exp{−j2πf2 (t − τ2 )} (9)
r2
N L = 10log10 (
σ2
)
1
0.2[m] 時の MUSIC スペ 2.0[m] 時の MUSIC スペ
クトル分布
クトル分布
ここでは音源 1 と音源 2 の角度差が 1[ ◦ ]
の時のノイズレベルの違いによる MUSIC ス
ペクトル分布を図７∼図 10 に示す。
図 7: NL=60[dB] 時の 図 8: NL=50[dB] 時の
MUSIC スペクトル分布 MUSIC スペクトル分布
(10)
ここで，f1 は音源周波数で 800[Hz] とし，f2
は 1700[Hz]，r1 は音源１までの距離，r2 は音
源 2 までの距離，τ1 ，τ2 は音源 1，音源 2 まで
のそれぞれの伝搬時間を示し，σ 2 はガウスノ
イズの平均電力を表す．
3.2
図 5: アレ イ素子間隔 図 6: アレ イ素子間隔
数値例
結果は，推定スペクトル強度 P の分布で表
示する．以下図 3∼図 6 に素子間隔を変えた
場合の結果を示す．ここでは音源 1 のみを用
いた．
図 9: NL=40[dB] 時の 図 10: NL=30[dB] 時の
MUSIC スペクトル分布 MUSIC スペクトル分布
4
まとめ
MUSIC 法を用いて音源の近傍における音
源探査を行った．その結果，以下のことがわ
かった．
• アレイ素子間隔が音源の波長の約 1/2 以
下では MUSIC スペクトルのピ ークは
1/2 以上のものと比べ鈍くなっているの
がわかる. だがすべての場合においてピー
クは音源位置を示す結果となっている．
図 3: アレ イ素子間隔 図 4: アレ イ素子間隔
0.002[m] 時の MUSIC ス 0.02[m] 時の MUSIC ス
ペクトル分布
ペクトル分布
• ノイズレベルが -60[dB] の時には二つの
鋭いピークが得られている．またノイズ
レべルを上げると徐々にピークは鈍くな
っていき-20[dB] の時においても二つの
ピークがみられるが音源位置からずれる
結果となっている．
平成17年度修論・卒論
通信システム工学講座
B − 1 フーリエ変換法を用いた BEM 非定常音場解析
稲見 正和 (通信システム) 1
はじめに
境界要素法（ BEM ）は任意形状の場の解析に
有効な手段である。しかし 、非定常場の解析では
必ずしも精度のよい結果が得られていない。そ
こで本研究では 、フーリエ変換と BEM を用い
た変換解法により厳密解が分かっている音響管
問題により検討する。
2
理論
二次元非定常音場の支配方程式は（ 1 ）で与え
られる。
1
∇2 p(t) − 2 p̈(t) = 0 in Ω (1)
c
½
p(t) = p̂(t) on Γ1
(2)
∂p(t)
Γ2
∂n ≡ q(t) = q̂(t) on ここで、p(t) は音圧、c は音速、q(t) はフラッ
クス、p̈(t) は時間の 2 階微分、( ˆ ) は既知量を
示す。
(1) 式にフーリエ変換を施すと、変換空間におけ
る支配方程式及び境界条件は以下のようになる。
∇2 P (ω) + k 2 P (ω) = 0 in Ω (3)
½
P (ω) = P̂ (ω) on Γ1
(4)
∂P (ω)
Γ2
∂n ≡ Q(ω) = Q̂(ω) on P(ω) は p(t) のフーリエ変換、k は波数で k
= ω/c である。変換後の (3) 式は、ヘルムホル
ツ方程式である。この方程式を、各周波数ごと
に解きある点での周波数成分を求め、逆フーリ
エ変換により周波数領域から時間領域に戻し過
渡波形を求める。
3
3.1

 P (ω) = P̂ (ω) on Γ1
∂P (ω)
= Q̂(ω) = 0 on Γ2 , Γ4
 ∂P∂n
(ω)
Γ3
∂n = kβP (ω) on 音響管モデルを図１に示す。
要素分割数は、x 方向：128 y 方向：8 であ
る。
また、Γ2 ,Γ4 は剛壁であり音の漏れは無いもの
とする。
Γ3 は開境界とし音の反射は無い。Γ3 を開境界と
する為に β なる音響アド ミッタンスを考える。
3.2
β の値の決定
今回扱う音響管問題を電気の分布定数回路と
考える。特性インピーダンス Z0 とすると
β=
1
Z0
(5)
なる関係がある。
電気回路と音波問題との間には、音圧 P は電圧、
フラックス ∂P
∂n は電流に対応している。この対
応に基づき、音響管の左端からある一定周波数
の音圧を入力し他端短絡の場合の βs 、他端開放
の場合の βf を求める。
β=
q
βs βf
(6)
この際、音響管が共鳴する周波数を除いて計
算した結果 |β| = 1.0 となった。
4
入力信号波形
計算条件
音響管モデル及び境界条件
図 2: 入力音圧波形
図 1: 音響管モデル
図 2 に示す音圧の信号を高速フーリエ変換
(FFT) して図 1 に示す境界 Γ1 に与える。
この時、サンプリング数 N = 1024 である。
(B − 1)
平成17年度修論・卒論
5
5.1
通信システム工学講座
6
計算結果
駆動面での結果
まとめ
結果を見て分かるように、厳密解と計算結果
の違いがほとんど 見られないほどの精度の良い
結果が得られた。誤差の原因には 、サンプ リン
グ点数、サンプ リング周期、また、分割した要
素の長さなどが関わってくる。
今回の研究で、音響管問題が変換解法によって
解析できることが示された。次に 、スピーカ−
などの任意の形状をした音響機器の過渡現象解
析への発展が期待できる。
参考文献
[1] ”パソコンによる境界要素法入門” 戸川隼
人/下関正義, サイエンス社
[2] ”信号処理入門” 雨宮好文/佐藤幸男, オー
ム社
図 3: 駆動面での厳密解と計算結果
5.2
[3] ”境界要素法による２次元非定常音場解析”
福田 正仁, 富山大学工学部修士学位論
文,1991
管の中央での結果
[4] ”騒音伝播の境界要素シュミレーション ” 藪
田 将貴, 富山大学工学部卒業論文,2003
図 4: 管中央での厳密解と計算結果
5.3
左端から 4 分の３での結果
図 5: 左端から 4 分の 3 での厳密解と計算結果
(B − 1)
平成17年度修論・卒論
通信システム工学講座
B-3 音声の可視化
大橋 理宏 (通信システム)
1
はじめに
音を見ることができれば様々な場面で役に立
つ。ホールや映画館などの音響施設を作る際、シ
ミュレーションにより音場分布を得ることによっ
てよりよい施設を建設できるはずである。本研
究は 、境界要素法と短時間フーリエ変換を用い
2 次元音響管内の音圧分布の解析を行った。
2
計算方法
1. 短時間フーリエ変換（ STFT ）をし 、各周波
数での音圧の実部と虚部、またその時間応
答を得る。
2. 周波数成分の実部と虚部を入力として音響
管内のヘルムホルツの問題を解き, それぞれ
の時間での音響管内の音圧分布を得る。
(a) ある周波数成分の実部と虚部で駆動す
る。まず音響管の境界値を求め、その
境界値で内点の値を求める。
(b) 全ての周波数成分で繰り返し 、音響管
内での音圧の和を得る。
(c) 各時間で繰り返し時間応答を得る。
図 2: 音声データ
4
解析結果
図 3 に音声信号に STFT を施したもの、図 4
に各時間での音響管内の音圧分布を示す。
以上の手順で音圧分布を求める。
3
計算条件
入力信号を音声データとし 、標本周期 512 、分
析窓長 0.032[s] 、時間ステップ数 10 、シフト長を
分析窓長の 1/2 、窓関数をハミング関数とした。
音響管モデルは図 1 である。横方向の分割数を
64 、縦方向の分割数を 4 、内点数は 441(63×7) と
し 、Γ1 を駆動面、Γ2 を剛壁としている。支配方
程式は領域 Ω 内で 52 P(ω)+k 2 P(ω)=0 となる。
ここで k は波数で k= ωc 、c は音速で 340[m/s],ω
は角周波数で今回は周波数を 280[Hz] で一定と
した。P(ω) は音圧にフーリエ変換を施したもの
である。また、要素は一定要素を用いている。
図 3: STFT の結果
図 4: 音響管内の音圧分布
5
図 1: 音響管モデル
まとめ
音響管内のそれぞれの時間での音圧分布を得
ることができた。問題点として波数計算の際に
一定周波数として計算しているため、周波数が
大きく変化するものには対応していない。周波
数が大きく変化する音声に対応できるものに改
善する必要がある。
(B-3)
平成17年度修論・卒論
通信システム工学講座
B-5 無限境界を有する電界問題の外挿法による計算精度改善
加納 正樹 (通信システム) 1
4
はじめに
雷雲による電界解析では 、雷雲電荷によって生ず
る地表付近の電界を正しく把握するが望まれている。
このような場の問題は３次元の場の問題であり、また
地表は必ずしも平坦ではないので数値解析が必要とな
り、境界要素法（ＢＥＭ）の適用される。また地表は
無限遠に及ぶので、半無限要素のような工夫が必要で
ある。本研究では本来無限に広がっている地表を有限
の境界要素領域として、地表電界を計算し 、無限遠の
影響をも考慮に入れた地表電界を計算する事を目的と
する。
2
解析結果
電荷の位置を (x,y)=(20,500) と (x,y)=(20,1000)
とし 、それぞれ 1[C] の点電荷がある時、外挿法によ
り求めた関数を表示し 、その結果を出力した。与えた
データはそれぞれ、４分割、８分割、１６分割、３２
分割の４つと、それに６４分割を加えた５つである。
電界の評価点は共に (x,y)=(-75,0)。
支配方程式及び境界条件
支配方程式と境界条件は以下のように表される。
52 φ +
Q
=0
ε
φ = φ̂ = 0
(in Ω)
(on Γ)
重み付き残差表現として次式が得られる。
Z
Z
Z
∂ψ
∂φ
Q
Ci φi +
φ
dΓ =
ψdΓ +
ψdΩ
∂n
∂n
Γ
Γ
Ω ε
(1)
(2)
(3)
この式を離散化し 、大地が接地 (式 (2)) していると考
えると、次のようになる。
G・E = b
(4)
ここで 、G は地形データにより決まる係数行列、E
は地上各節点の電界値ベクトル。b は各節点と点電荷
との関係（各節点への電荷による影響）から求まるベ
クトル。
3
表 1: 解析結果
電荷の位置 (x,y)
厳密解 [V/m]
計算値, データ 4 つ [V/m]
計算値, データ 5 つ [V/m]
解析モデル
以下の図１のように有限の領域に区切り、それぞれ
の領域においての電界を求める。求められた電界値が図
２のようなグラフになる場合 E = a(1/n)2 +b(1/n)+c
と近似できると仮定をし 、最小二乗近似を用いて、関
数を導く。
(20,500)
0.000314
0.000201
0.000292
(20,1000)
0.000159
0.000168
0.000226
ここで、電界の値は法線方向のベクトル値で、比較
のため定数である 1/ε は除いてある。
5
まとめ
外挿法により求まった結果では大まかな形ではその
特性は得られたが、与えたデータの数、電界の評価点
等によって誤差があらわれるので、これらの改善が必
要である。
参考文献
[1] C・A ブレビア著:『境界要素法入門』培風館、1980
[2] 中山正敏著:『電磁気学』培風館、2001
[3] 杉浦 洋著:『数値計算の基礎と応用』朝倉書店、1997
(B-5)
平成17年度修論・卒論
通信システム工学講座
B-7 コンクリート壁内欠陥探査の FDTD 音響解析
澤城 雅俊 (通信システム)
1
はじめに
２次元 FDTD 法を用いてコンクリートの欠陥部
を音によって検出できるかど うかを数値シミュレー
ションにより検討した。FDTD 法は時間領域差分方
のことで定式化が簡単で優れた精度を持っているの
でこの方法を用いた。
2
ぼ変化は無いが、B 点の観測点では穴による反射波
の影響が見られる。また反射波は音を送波してから
0.4[ms] にみられる。この時間は計算による理論値
とほぼ一致するため、穴による反射であると推測さ
れる。
計算条件
(a);A 点
(b);B 点
図 2: 穴なしの場合の音 図 3: 穴なしの場合の音
圧の時間変化
圧の時間変化
図 1: 解析モデル（ A は音源、B は受信点）
(a);A 点
解析領域とし て 、図１のよ うな 平面を 200 ×
120[cell] に分割し、空気中の一部をコンクリートにし
たものと、そのコンクリートに欠陥を仮定した穴を
開けたものの２種類を考えた。A 点の位置を (80,52
∼60) 、B 点の位置を (79,63) とする。セルサイズを
5.5 × 10−3 [m] 、時間ステップを δt = 5.5 × 10−7 [s]
とし 、いずれも courant の安定条件を満たすよう
に設定した。また空気の密度を 1.2[kg/m3 ] 、コン
クリートの密度を 2300[kg/m3 ] 、空気中での音速を
340[m/s] 、コンクリート中の音速を 3500[m/s] とし
た。モデルの境界面には Mur の１次吸収境界条件
を用いた。 音源は正弦波１周期分に相当するパル
ス音源とし 、粒子速度駆動で気体に垂直に進行する
平面波に近いものにした。
3
(b);B 点
図 4: 穴ありの場合の音 図 5: 穴ありの場合の音
圧の時間変化
圧の時間変化
4
まとめ
FDTD 法を用いて欠陥（コンクリート中の穴）の
探査シミュレーションを行った。ただ、空気とコン
クリートでは音速や媒質密度の差が非常に大きいた
め適切な courant の安定条件を選ぶことが困難だっ
た。また穴の大きさが徐々に小さくなるにつれて判
別しにくくなっていった。これらについてさらなる
検討が必要である。
参考文献
計算結果
図２∼５、に結果を示す。図２、図３は穴なしの
場合で、図４、図５には穴ありの場合の結果を示し
てある。いずれも図１に示す A 点及び B 点での音
圧の時間変化を示している。A 点は音源なのでほ
(B-7)
[1] 宇野亨著,”ＦＤＴＤ法による電磁界およびアン
テナ解析”, コロナ社,1998
[2] 奥 田 奈 生 著,”コン ク リート 壁 内 欠 陥 探 査の
FDTD 音響解析”, 富山大学卒業論文,2004
平成17年度修論・卒論
通信システム工学講座
B ー１０ 最近の雷害対策研究のデータベース作成
高堂 良平 （ 通信システム） 1
はじめに
ルとしてはリレーショナルデータベースが一覧の見
易さという点で適していると判断した。
毎年、雷害による停電や器物の破損が大きな被害
をもたらしている。そこで今後の雷害対策の発展の
ためにも、研究されている雷害対策について調査し 、
データベースを作成した。
2
データベースモデルの検討
データベースには個々にモデルが示されている。
そこでリレーショナルデータベースやハイアラキカ
ルデータベース、ネットワークデータベースの中か
ら、どのモデルを採用するかを検討した。図１に学
生の成績を用いたモデルの使用例を示す。
3
作成したデータベース
まず、データを蓄積した。データは論文誌から雷
害対策研究に関する該当論文を探し 、データとして
取り入れた。また、英語の論文や雷害対策に関する
公開特許からもデータを取り入れた。
図 2 に作成したデータベースを示す。
採用の検討に考慮した条件としては
・内容が見やすいデータベースである
・得たい情報に容易に到達できる
図 2: 作成したデータベース
図 1: 学生の成績を用いたモデルの使用例
それぞれの特徴として 、リレ ーショナルデータ
ベースモデルは主キーとなる言葉に関係付けられる
ものをデータベースシステムから表として取り出す
もので、ハイアラキカルデータベースやネットワー
クデータベースはポインタやチェイン、多重リンク
を使ったデータベースモデル（図１参照）で、複雑
にリンクが張られている。質問処理の速度やコス
トの点ではリレーショナルデータベースより優れて
いる。
4
まとめ
本研究では、主に過去 5 年間の雷害対策研究とそ
の問題点に関する調査をし論文 65 編、公開特許 74
件の計 139 件でデータベースを作成した。
データベースの特徴としては
・一画面で論文等の内容が見えるよう配慮した
・細かな項目におけるキーワード 検索への対応が不
十分である
前述の条件から、求めているデータベースのモデ
(B-１０)
平成17年度修論・卒論
通信システム工学講座
B-11 モバイル通信用通信プロトコルの開発
柘植 功蔵 (通信システム) 1
はじめに
現在、世界中でモバイル通信の研究は盛んに行
われている。特にモバイル通信用通信プロトコル
は発展段階にありより性能の良い物を提案しよ
うと日々研究されている。そこで、その前段階で
ある有線でのシミュレーションを The Network
Simulator Verstion 2 で行った。
2
１７の方向にルータのノード ８，９を介する。
・ 帯域、伝播遅延は図１の通り。
・ スタートしてからすぐに通信は開始され、５
秒後に通信は終了する。
3
解析結果
シミュレーションモデル
図 2: 障害のない場合の転送率
図 1: シミュレーションモデル
図１のようなモデルでシミュレーションした。
ノード ０∼７がデータを送る送信者。ノード １
０∼１７がデータを受ける受信者とし 、ノード
８，９をルータとする。イメージとしては送信
側の半分は常にデータを送り続けるヘビーユー
ザーで残り半分は一般的なユーザーがいるネッ
トワークで転送率を調べる。その他、細かい条
件は次に示す。
∼条件∼
・ 全ての接続は全二重接続とする。
・ 送信ノード０∼３は UDP プロトコルを使用。
・ 送信ノード ４∼７は TCP プロトコルを使用。
・ 受信ノード １０∼１３は NULL エージェント
とし 、受け取ったデータをすぐに破棄。
・ 受信ノード １４∼１７は SINK エージェント
とし 、データを受け取ったら受信確認を行う為
に送信側に返信。
・ノード０∼３のパケットジェネレータには CBR
トラヒックを利用。
・ノード ４∼７のパケットジェネレータには FTP
トラヒックを利用。
・ 通信はノード ０→１０、１→１１、２→１２、
３→１３、４→１４、５→１５、６→１６、７→
図 3: シミュレーション結果
4
まとめ
解析結果よりヘビーユーザーが帯域を圧迫す
るため一般ユーザ−は満足のいくスループット
を得られないということが分かった。この問題
を解決するためには通信の公平性が必要だと思
われる。有線でこの結果なので無線になったら
もっとこの問題は大きくなると考えられる。
今後は 、こういった問題を解決したモバイル
通信用通信プロトコルを開発していきたいと考
えている。
(B-11)
平成17年度修論・卒論
通信システム工学講座
B-16 道路沿道建物の騒音レベル上昇効果
橋本 博志（ 通信システム）
1
3
はじめに
私達の周りには少なからず騒音というものが
存在しており、特に道路の騒音は私達の生活に
も大きな影響を与えている。騒音伝搬を数値的
に解く方法としては一般的には差分法、有限要
素法、境界要素法などがあげられる。本研究で
は 、道路わきに高さの異なった建物を置きその
時の騒音伝搬の仕方を境界要素法を用いて解き、
その対応策などを考える。
2
解析結果
図 2: 太線部がすべて完全反射境界面
解析モデル
解析モデルは下図１のようなモデルとした。
騒音源は片側車線を走っている車を想定する。
実際の音は空気を伝搬するとき減衰するが本シ
ミュレーションではこれを考慮していない。
音速は 340[m]
音源 X からの音圧は 0.1[N/m2 ]
周波数は車の排気音を考えた 250[Hz]
要素は一定要素を用い要素数は 660, 内点は
10800 個をとった。
図 3: 左の障害物のみ吸収性境界面（吸音率 0.6 ）
境界積分方程式を以下に示す。
Z
Z
∂p(ξ)
∇2 p(ξ)ci +
Gij dΓ − p(ξ)Hij dΓ = P i (ξ, X)
∂n
P i (ξ, X) = P0
1 (2)
H (kr(ξ, X))
4j 0
ci は境界の形状によって決まる定数、 P0 は音
(2)
源からの音圧、H0 は第 0 次第 2 種ハンケル関
数となっている。
図 4: 右の障害物のみ吸収性境界面（吸音率 0.6 ）
4
図 1: 解析モデル
解析結果は相対的な表示の騒音レベルと呼ば
れるもので表示し 、騒音レベルは下式で定義さ
れている。
Ｌ = 20log10 (P/P0 )
[dB]
P は実効音圧、P0 は基準の実効音圧を表す。
P0 =0.00002
[N/m2 ]
P0 の値は若い人が両耳を使って聞くことがで
きる最小の音圧に近い値となっている。
結論
図 2 と図 3,4 より建物に吸収境界面を使用し
た時、騒音レベルが下がることがわかった。
次に 、建物の片方のみに吸収境界面を使用し
た時の騒音レベルを考察する。
図 3,4 を比べ、まず大小の建物間での騒音レ
ベルの影響を考える。そうすると、大きい建物
(図 4) を吸収境界面にした方が騒音を下げる効
果があることがわかる。
そして両建物の外側に与える騒音レベルの影響
を考えると、小さい建物 (図 3) を吸収境界面にし
た方が騒音を下げる効果があることがわかった。
(B-16)
平成17年度修論・卒論
通信システム工学講座
B-19 電磁波人命探査のＦＤＴＤシミュレーション
浜 智洋 (通信システム)
1
はじめに
近年、地震災害が多く発生している。その被害の
中では家屋の倒壊や土砂の滑落により人が生き埋め
にされてしまうという例が多くある。そのような場
面で生きた人を発見する手段として電磁波人命探査
装置が用いられる。
しかしながらこのシステムの詳細な内容について
は明らかにされていない。そこで電磁界解析により
人命探査のシミュレーションを行う研究がある。そ
の電磁界解析は FDTD 法を用いることで容易に行
うことができる。
FDTD 法においてセルサイズは、解析精度を決め
る重要な要素である。そこで２次元における２媒質
間の電磁界強度の変化をセルサイズを変化させて解
析し 、それによる精度の違いを比較検討した。
2
表 1: 大地の条件
比誘電率
乾地
4
湿地
20
導電率 [S/m]
0.00001
0.001
の sin 波を与えた。ここで t は時間とし 、周波数
f の値は 3[GHz] とした。
3
解析結果
解析結果を図２∼５に示す。
解析条件
解析モデルを図１に示す。計算条件は、解析空間を
1[m] × 1[m] の平面とし、セルサイズは、1[m] をそれ
ぞれ一辺のセル数で割った値とした。時間ステップ
は、今回解析を行った中での最小のセルサイズ ∆x =
3.125 × 10−3 [m] における Courant の安定条件を考
慮して、∆t = 2.00×10−12 [s] とした。媒質定数とし
て、空気部分では誘電率 ε0 = 8.854 × 10−12 [F/m] 、
導電率 σ0 = 0[S/m] とした。そして、領域の半分に
は大地を配置し 、電磁界強度変化の解析を行った。
図 2: 乾地での電界強度 図 3: 乾地での磁界強度
図 4: 湿地での電界強度 図 5: 湿地での磁界強度
4
まとめ
解析結果からもわかるようにセル数を多くするほ
ど細かい解析結果を得ることができることがわかる。
しかしながら、媒質によっても解析結果の精度に
変化が見られることから、媒質に合わせたセル分割
数を選択することが望ましい。
参考文献
[1] 宇野 亨著：
『ＦＤＴＤ法による電磁界およびア
ンテナ解析』コロナ社、1998
図 1: 解析モデル
また、大地の媒質定数は表１のように与えた。
波源としては
Ez = sin(2πf t)
(1)
(B-19)
平成17年度修論・卒論
通信システム工学講座
B-32 インピーダンス CT の有限要素シミュレーション
姜 東勲 (通信システム) 1
はじめに
インピーダンス CT では対象物に電流を注入し 、
その周囲の電圧を測定することにより、内部の導電
率分布を画像として再構築させるものである。特徴
として X 線 CT や MRI に比べ装置が簡単で安価に
構想可能、入力が電流であるので X 線による被爆の
危険性がなく人体への影響は非常に小さいが、再現
性の悪さから実用化はされていない。そこで SPMM
を用いて推定精度改善を検討した。
2
計算手順
各電極対に対するインピーダンス Z を求める。ま
ず各要素に初期伝導率を与え、
(n)
∆Zi
(n)
= Ẑi − Zi
=
M
X
(1)
j=1
(i = 1, 2, · · · ,
図 1: 求めたい導電率分布
(n)
Sij ∆σj
N (N − 1)
)
2
(n)
を解き、更新量 ∆σj が得られる。ここで、i は電
極対パターン、Sij は感度行列、j は要素番号、n は
反復回数、M は要素数、σ は導電率、N は対象物表
面の電極数 (節点数) である。SPMM ではこの更新
(n)
量 ∆σj にサブスペース法により先験情報を用い、
(n)
∆σj の要素が先験情報に ”適合”していればその要
素の値を修正し 、”不適合”はらその要素の値をその
まま用いる。そして、
(n+1)
σj
(n)
= σj
(n)
+ α∆σj
(j = 1, 2, · · · , M )
(2)
により反復計算を行い、
C=
M
X
(n)
|∆σj |
図 2: Marquardt 法 (推定誤差＝ 4.86 ％)
(3)
j=1
(n+1)
が十分小さくなれば収束したとみなし 、σj
を推
定値とする。ここで、α はスケール因子である。
3
結果
シミュレーション結果を図２、３に示す。SPMM
を用いることで誤差を少なくし精度を良くした。個
人差はあるが、臓器のだいたいの位置や形状は似て
いるため、先験情報が精度改善に繋がった。
4
まとめ
今回はＳＰＭＭ法を用いることで精度が改善した。
しかし 、実用化するためには更なる精度改善が必要
である。
(B-32)
図 3: ＳＰＭＭ (推定誤差＝ 1.72 ％)