JP2020034379A - 振動騒音伝達経路推定システム及び伝達率推定方法 - Google Patents

Abstract

【課題】振動騒音伝達経路推定システム及び伝達率推定方法を提供する。【解決手段】複数の振動騒音センサーからの信号をデシベル値として測定時刻ごとに記憶し、同一測定時刻における複数の振動騒音センサーからの信号を入出力データ１４ａとし、入出力データと伝達率を用いて振動騒音現象を表現するモデル式に適用し、同一測定時刻の入出力データに対応する伝達率のデシベル値に対する尤度１６１〜１６３を計算し、計算した尤度を事前に想定している伝達率の事前確率分布に乗算することで伝達率の推定値の事後確率分布を得、同一測定時刻を逐次変更して事後確率分布を繰り返し求めることで振動騒音の伝達経路を推定することを特徴とする振動騒音伝達経路推定装置。【選択図】図５

Description

本発明は、振動騒音問題の発生防止に向けた効果的な事前施策実施の選定に供する振動騒音伝達経路推定システム及び伝達率推定方法に係り、しいては、その施策の効果について確率論的に扱うことで、コスト期待値を判断基準とした費用対効果の定量的判断方法を可能とする振動騒音伝達経路推定システム及び伝達率推定方法に関するものである。

製品から発生する振動や騒音は製品使用時の快適性を悪化させる要因として広く認識されている。そのため、振動や騒音の発生源と伝達経路を把握し、その結果を元に現状の構造案に対する変更とその効果の事前予測が必要となる。

一方で、一般に振動騒音を抑止するために実施される防振・防音施策は、コストや重量を増大させるばかりか、他の性能とトレードオフの関係にあることが多い。そのため、振動騒音を抑止するためのみを目的とした施策を、むやみやたらに実施することは忌避される。例えば、あえて事前の防振・防音施策は適用せず、問題発生の後に事後対策を行なうという戦略も、それが費用対効果の観点から最適であると判断させるのであれば成立し得る。また逆に、問題発生後に実施した事後対策の費用が、事前に施策を実施した場合の想定費用よりも高額になることも、当然ながら避けなければならない。

このように、究極的には要求される振動騒音に関する仕様を必要十分なぎりぎりのレベルで達成し、且つ施策費用や重量、他性能への影響なども含めた総合的な広義のコストを最小限にして実現できる施策の適用が求められる。しかしながら、その判断は実際には極めて難しい。なぜならば、そもそも現行製品の騒音がどこから発生し、どのように伝達しているかを示す伝達経路と、その割合を示す寄与の判定が難しく、結果、現状構造に対する構造変更の効果を示しにくいためである。

このような振動騒音の事前対策に関する判断支援の技術に関しては特許文献１がある。特許文献１では既存の製品の実稼働時の振動騒音測定データをもとに、未知である伝達関数を同定し、さらに各部位からの騒音伝達の寄与を再構築する技術が提供されている。

特許第４０２８５６２号

製品が発生する振動や騒音を抑制するには、まず音源や加振源および伝達経路の把握が必要であり、その実現には一般には現行製品の実測データを用いた音源特性や伝達特性の同定が行なわれ、さらには支配的な伝達経路と各経路からの寄与の大きさを示す寄与率の判定が行なわれる。

背景技術で説明したとおり、一般に振動騒音を抑止するために実施される防振・防音施策は、コストや重量を増大させるばかりか、他の性能とトレードオフの関係にあることが多い。そのため振動騒音に関する仕様を必要十分なぎりぎりのレベルで達成し、且つ施策実施費用や重量、他性能への影響なども含めた総合的な広義のコストを最小限にして実現できる施策の適用が求められるが、そのためには、まず初めに、現行製品の騒音の寄与を把握する必要がある。

ここではまず、図１１から図１７を用いて従来における騒音伝達解明手法の一例を説明するとともに、本発明において対応すべき課題について具体的に説明する。

図１１は鉄道車両の車内騒音を例に、車内騒音に対する伝達経路とその寄与率の例を示している。

図１１の左には車内騒音に対する伝達経路の例を示しており、車内の特定個所に設けたセンサーにより、車内各所（屋根、荷棚、窓、側、床など）の振動や騒音を測定する。

騒音の測定結果は解析により図１１の右に示すような寄与率の例として把握され、モニタなどに表示される。図１１の右の例における表示では、横軸に騒音をスペクトル解析した時の周波数を、縦軸に車内各所ごとの騒音寄与率を示している。

このようにして、次期車両の開発にあたっては現行車両の徹底した計測による実測データをもとに、図１１右に示すような寄与率を判定し、次期車両において静音化設計上重要な個所を把握する。たとえば、図１１の右の解析結果の場合、低周波では１点破線で示した床からの寄与が大きく、一方で高周波では破線で示す屋根からの寄与が大きい。そのため、次期車両では、床構造は高周波の遮音性を犠牲にしてでも低周波における遮音性を重視させ、逆に屋根に関しては低周波の遮音性はある程度妥協してでも高周波の遮音性を確保することが有効などと判断される。

図１１は鉄道車両の場合の例を示したものであるが、他の製品においてもこのような指針が得られることが次期製品の開発において極めて重要となる。

図１２は伝達率や寄与の理解を容易にするために導入する「伝達経路モデル」の概念を示す図である。図１２は、図１１の左に示した鉄道車両車内騒音の伝達経路について、屋根、荷棚、窓、側、床など内装各部位の振動が音として放射され、それぞれある伝達率で車室内へ伝わったものが車内騒音へのそれぞれの寄与として表現されることを示している。ここでは、図１１に示す伝達経路について、入力信号と伝達率の線形一次結合として出力信号が求まると仮定した伝達経路モデルで表現した例を示している。

具体的には（１）式に示すように、ｉ番目の内装の振動加速度、ないしは振動速度をｘ_ｉ ^（ｎ）としたとき、車内騒音ｙ^（ｎ）は、前記ｘ_ｉ ^（ｎ）とそれぞれの伝達率ａ_ｉを一次結合させた形で表現できるものと考える。なお、右肩の添え字（ｎ）は条件番号を表しており、たとえば走行区間や走行時刻などにあたる。

ここで、ｘ_ｉ ^（ｎ）およびｙ^（ｎ）は既存車両の走行試験などの実稼働データから得られ、伝達率ａ_ｉが未知である。ここで伝達率ａ_ｉを知ることができれば、前記ｘ_ｉ ^（ｎ）とａ_ｉの一次結合としてあらわされるｙ^（ｎ）の各項の大きさとして各部位からの寄与を表されると考えることができる。このため、ここでは図１１の屋根、荷棚、窓、側、床など内装各部位に設置したセンサーによりｉ番目の内装の振動加速度、ないしは振動速度ｘ_ｉ ^（ｎ）を得、また図１１の車内騒音の部位に設置したセンサーにより車内騒音ｙ^（ｎ）を得るものである。なお各種センサーから得られる情報のうちｘ_ｉ ^（ｎ）は振動原因側要因であり、ｙ^（ｎ）は結果側要因の騒音であることからｘ_ｉ ^（ｎ）を入力データと称し、ｙ^（ｎ）を出力データと称することがある。またｘ_ｉ ^（ｎ）とｙ^（ｎ）を総称して入出力データと称するものとする。なお，ｘ_ｉ ^（ｎ），ｙ^（ｎ），ａ_ｉともに周波数依存性があり，各周波数のデータそれぞれについても先述と同様な線形一時結合性が成立していると考える。以下ある周波数に対応する入出力データをｘ_ｉ ^（ｎ），ｙ^（ｎ），入出力データをａ_ｉと表現する。

伝達率ａ_ｉの同定方法としては特許文献１などで適用されている最小二乗法による線形重回帰分析があげられる。具体的には、（２）式に示すように、上記右肩の添え字（ｎ）で示した複数の条件下における入出力データｘ_ｉ ^（ｎ）およびｙ^（ｎ）を取得しておけば、条件に依存しない伝達率ａ_ｉに対する逆同定問題として解くことができ、得られた伝達率を再び上記一次結合表現に戻せば、各部位からの寄与が各項の大きさとして求めることができる。なお、特許文献１では重回帰分析の前に主成分分析ないしは特異値分解などを行うことで、伝達率算出の安定化やノイズの除去などを行っている。

なお（２）式において、ｙ^（ｎ）は複数の多点同時計測結果であり、伝達率ａ_ｉを一定と仮定したものである。また伝達率ａ_ｉに対する逆同定問題は、擬似逆行列を用いた伝達関数の逆計算であり、入出力データｘ_ｉ ^（ｎ）およびｙ^（ｎ）から伝達率ａ_ｉを求め、求めた伝達率ａ_ｉを（１）式に代入してｙ^（ｎ）を求めることで伝達経路を求めるものである。この伝達率ａ_ｉに対する逆同定問題による手法は、最小二乗法による線形重回帰分析と同義である。

ここで説明をさらに簡単にするため、以後、図１３に示す２入力１出力モデルを用いて説明する。図１３に示す２入力１出力モデルＭでは、二つの入力ｘ_１とｘ_２に対してそれぞれ伝達率ａ_１、ａ_２が掛け合わされたものの合算値として出力ｙが得られるとする。このモデルＭは、図１２における屋根や床がｘ_１、ｘ_２に、車内騒音がｙに相当すると考えてよい。ただし、実際の測定データはこれらにノイズｅ_ｘ１、ｅ_ｘ２、ｅ_ｙが混入された情報を観測することとなる。そのため、上記モデルｘ_１、ｘ_２、ｙと区別した観測系Ｏとしてｘ_ｍ１、ｘ_ｍ２、ｙ_ｍと表記する。

なお以降の説明においては、観測系Ｏにおいて計測されたｘ_ｍ１、ｘ_ｍ２、ｙ_ｍのことを入出力データ１４ａとして、表記、説明する。また観測系Ｏにおいて計測された入力データ１４ａであるｘ_ｍ１、ｘ_ｍ２は、図１１の屋根、荷棚、窓、側、床など内装各部位に設置したセンサーにより検知したｉ番目の内装の振動加速度、ないしは振動速度ｘ_ｉ ^（ｎ）に対応し、観測系Ｏにおいて計測された出力データ１４ａであるｙ_ｍは、図１１の車内騒音の部位に設置したセンサーにより検知された車内騒音ｙ^（ｎ）に対応している。

観測系Ｏ（センサー）において計測されたこれらの入出力データ１４ａは、伝達経路解析のために一般的には計算機で構成される伝達経路解析システム内のデータベースＤＢに取り込まれる。

ところで振動騒音の分野では各種の運転稼動条件における振動加速度や振動速度、音圧などの物理量の変動値の二乗平均値をパワーとしてとらえ、このパワーの常用対数に１０を乗じた、いわゆるデシベル値で表現することが多い。なお「デシベル値」は，「レベル」と称されることもある。以下特に断りがない限り「デシベル値」で統一する。このことからたとえば、図１３下部に記載のデータベースＤＢには、計測時刻ごとの入出力データ１４ａ（ｘ_ｍ１、ｘ_ｍ２、ｙ_ｍ）が、デシベル値１３ａ（Ｘ_ｍ１、Ｘ_ｍ２、Ｙ_ｍ）として記憶される。なお大文字はデシベル値を、小文字はパワー値を表している。

１３ｂで示した（３）式は、モデルＭを表現するモデル式であり、モデル式は解析の対象とした２入力１出力モデルＭのモデル式と、パワー値をデシベル値に変換するときの変換式を含む例を示している。

この状態、すなわち、入出力データ１４ａのデシベル値の時系列データとして取得されたデータベースＤＢのデータ１３ａと、想定される（３）式のモデル式１３ｂから伝達率ａｉないしそのデシベル値であるＡｉが算出されれば、前述の手法により寄与が求まることとなる。図１３下の右側に示すグラフ１４ｄは、伝達率ａｉのデシベル値であるＡｉ（ここではＡ_１とＡ_２）を縦横軸として伝達率Ａｉの推定結果を表示したものであり、推定結果の領域内に伝達率Ａ_１、Ａ_２が存在していると考えることができる。なお、ここでは本発明により伝達率ａｉのデシベル値である（Ａ_１、Ａ_２）が確率推論されることを示しているが、詳細は本発明の実施例において後述する。

図１４は特許文献１などで適用されている従来技術である最小二乗法による線形重回帰分析を用いた伝達率同定手法について図１３の２入力１出力モデルＭに適用した場合について視覚的に説明した図である。ここで図の各軸（ｘ_ｍ１、ｘ_ｍ２、ｙ_ｍ）は入出力データ１４ａのパワー値表現である。これは重回帰分析を用いるうえでの前提条件であるモデル式１３ｂの線形性がパワー表現で達成されるためである。

図１３に示した２入力１出力モデルＭの場合、データベースＤＢに含まれる各入出力データ１４ａのひとつひとつは、図１４の３次元空間（ｘ_ｍ１、ｘ_ｍ２、ｙ_ｍ）に配置された点の座標として表現され、最小二乗法による線形重回帰分析を用いた分析では、二乗残差の合計が最小となるように、すなわち、３次元空間上に配置された各データ点群のほぼ中央を通るように近似平面ＰＬが選択され、この近似平面ＰＬのそれぞれの軸に対する傾きが、すなわち図１４におけるａ_１、ａ_２が伝達率ａｉとして算出される。

ただし、この手法には以下に示す３つの観点で課題がある。

一つ目の課題は算出された伝達率ａ_１、ａ_２が負値になる場合があることであり、非負値となることを想定していた設計者がその結果の解釈に悩むことがしばしばある点である。たとえば図１３に示した２入力１出力モデルＭの場合、設計者は暗黙のうちに伝達率ａ_１、ａ_２が非負値になることを想定しており、ｙからｘ_１ないしはｙからｘ_２へのパワーの逆流は想定していない。しかしながら、上記に示す重回帰分析では算出される伝達率ａ_１、ａ_２が非負値になることを保証しない。

図１５はこのように暗黙の想定と違って伝達率ａ_１、ａ_２が負値になる場合を説明する図であり、各データ点群をもとに二乗残差が最小となるように近似平面ＰＬを定義すると、伝達率ａ_１、ａ_２のうちのひとつａ_１が負値となってしまった状況を表している。この場合、図１３におけるパワーの流れを考えると、ｘ_２から入力したパワーは一部ｘ_１に吸い取られ、その残りとしてｙに伝達していることと解釈される。しかしながらこのような結果は経験的に知っている現象と容易に折り合いがつかないことがしばしば発生し、これを回避するために、伝達率算出結果を修正して所望のように非負値の範囲で探査するために、幾らか非合理的な制約や仮定を設けざるを得ない場合が発生する。理想的には原理上非負の範囲で探査することのできる手法の適用が求められる。

二つ目の課題は機器が高負荷で稼動している状態で取得されたデータと低負荷で稼動している状態で取得されたデータのバラツキの取扱が等価でないことであり、算出された伝達率ａ_１、ａ_２はほとんどが高負荷稼動状態でのデータをもとに決定されてしまい、低負荷稼動状態のデータに伝達率ａ_１、ａ_２を同定する重要な特性があったとしてもそれがほとんど無視されてしまうことである。これを説明したのが図１６である。

図１６では説明を容易に行う為に、入力をｘ_ｍ、出力をｙ_ｍとした１入力１出力モデルＭでのデータ点群の例を示しており、これらデータ点群のうち回帰直線ｙ＝ａｘからいくらか偏差をもったデータ１（ｘ_１、ａｘ_１＋ｅ_１）とデータ２（ｘ_２、ａｘ_２−ｅ_２）について、どちらの方が回帰直線ｙ＝ａｘの傾きａの算出結果に対して大きく影響を与えたかを議論する。なお、ｘ_１＞ｘ_２であり、データ１は高負荷稼動時のデータを、データ２は低負荷稼動時のデータを模擬している。

また、回帰直線ｙ＝ａｘからの偏差ｅ_１、ｅ_２について、入力に対する比率（ここでは便宜上「比偏差」と呼ぶこととする）をそれぞれｒ_１＝ｅ_１／ｘ_１、ｒ_２＝ｅ_２／ｘ_２とし、今、ｒ_１≒ｒ_２であったとする。この場合、入力量で規格化した偏差である比偏差は、すなわち回帰直線ｙ＝ａｘの傾きに対するズレの感度を表しており、これはどちらもほぼ同じであることとなる。しかしながら、最小二乗法による線形重回帰分析を用いた場合では、入力の量で規格化していない偏差の総和を最小化するため、データ１のほうが圧倒的に回帰直線の傾きａの算出結果に対して影響が大きい。

すなわち高負荷稼動状態におけるデータのバラツキの影響が、低負荷稼動状態におけるデータのバラツキよりも大きく影響してしまい、先述のように低負荷稼動状態のデータに重要な特性が表れていたとしてもそれが無視されてしまうことが懸念される。また先述のように、振動音響分野では、物理量の変動値の二乗平均値をパワー値としてとらえ、このパワー値の常用対数に１０を乗じた、いわゆるデシベル値で表現することが多いため、どちらかというと誤差・バラツキはデシベル値と親和性の高い、上記「比偏差」をもとに議論するべきである。

しかしながら従来手法の線形重回帰分析ではモデル式１３ｂの線形性が前提条件となっているため、デシベル値そのままでの探査はできず、線形重回帰分析を用いる以上はどうしてもパワー値での探査にしかならざるを得ず、前述のような高負荷稼動状態のデータのバラツキの影響を大きく受けることとなる。

三つ目の課題は重回帰分析を用いる際の前提条件である入力データの独立性が、実際の振動騒音計測結果では成立しにくいことである。図１７はこれを示したものであり、たとえば図１１、図１２で挙げた鉄道車両の車内騒音などの場合、データの大部分が走行速度に依存する（走行速度のべき乗に比例するなど）として説明されてしまうため、入出力データ１４ａを空間上にプロットすると、ほとんど平面状にはならず、広がった線状ないしは曲線状に分布し、ここから近似平面を切り出すこと自体が原理的に不可能となる。

にもかかわらず、特許文献１などに記載の方法では、あくまでデータは平面上に分布しているとの仮定の下、そこからの偏差はノイズとして認識され、除去する処理がなされる。この場合、データの従属性によって引き起こされる伝達率同定結果の信頼性の低さは、その結果に保持されない。

一方で設計者は、入出力データ１４ａの独立性の低さにはなかなか気づきにくく、結果的に決定論的に逆同定された伝達率を特に疑いなく扱うこととなってしまう。また、たとえ原理的に困難があったとしても、設計上何某かの結論が求められる場合も多々ある。このような場合、データから無理に画一的な数値として表される伝達率を算出するのではなく、伝達率の推定結果について確率論的な幅を持たせた取り扱いをすることにより、その確率が結果の信頼性を意味するように保持させておくべきである。

本発明は上記の従来技術の「伝達率が意図せず負値になる」「バラツキの取扱が高負荷稼動状態と低負荷稼動状態で等価でない」「バラツキをノイズとして切り捨ててしまうため不確実性ないしは結果の信頼性を議論できない」という３つの観点における課題を解決することができる振動騒音伝達経路推定システム及び伝達率推定方法を提供するものである。

以上のことから本発明においては、「複数の振動騒音センサーからの信号を格納する格納装置と、格納装置に格納されたデータを用いて振動騒音を演算する演算装置と、演算装置の演算結果を元に、振動騒音の伝達経路推定結果を表示する表示装置からなる振動騒音伝達経路推定装置であって、格納装置は、複数の振動騒音センサーからの信号をデシベル値として測定時刻ごとに記憶し、演算装置は、同一測定時刻における複数の振動騒音センサーからの信号を入出力データとし、入出力データと伝達率を用いて振動騒音現象を表現するモデル式に適用し、同一測定時刻の入出力データに対応する伝達率のデシベル値に対する尤度を計算し、計算した尤度を事前に想定している伝達率の事前確率分布に乗算することで伝達率の推定値の事後確率分布を得、同一測定時刻を逐次変更して事後確率分布を繰り返し求めることで振動騒音の伝達経路を推定することを特徴とする振動騒音伝達経路推定装置」のようにしたものである。

また本発明は、「複数の振動騒音センサーからの信号をデシベル値として測定時刻ごとに記憶し、同一測定時刻における複数の振動騒音センサーからの信号を入出力データとし、入出力データと伝達率を用いて振動騒音現象を表現するモデル式に適用し、同一測定時刻の入出力データに対応する伝達率のデシベル値に対する尤度を計算し、計算した尤度を事前に想定している伝達率の事前確率分布に乗算することで伝達率の推定値の事後確率分布を得ることを特徴とする伝達率推定方法」のようにしたものである。

本発明により、上記の従来技術の３つの観点における課題を解決し、設計者の意図通り非負値の伝達率が得られ、しかもバラツキの影響を稼動状態に依存せず等価に扱いながら、バラツキを不確実性として受け入れながら結果の信頼性を確率という形で保持させた状態での伝達経路推定結果の取得を実現することができる。

デシベル値１３ａのうちＸ_ｍ、Ｙ_ｍの関係を表記したグラフを示す図。 各周波数で得られた伝達率のデシベル値の確率分布を示す周波数−伝達率特性グラフを示す図。 確率推論の一手法であるベイズ統計の原理を示す図。 ベイズ統計のなかで重要な概念であるベイズ更新を示すために導入した「コイントス問題」を説明する図。 図３に示したコイントス問題において、得られた結果からその結果の尤度を求め、これを事前確率分布に乗算することで事後確率分布を算出するベイズ更新の概念を説明する図。 図１３に示した２入力１出力系伝達経路モデルにおいてベイズ更新を適用した場合について説明する図。 本発明で採用するベイズ統計に基づく伝達率推定方法について、従来方法と比較してその優位性を示す図。 本発明の振動騒音伝達経路推定システムの機器構成を説明する図。 振動騒音伝達経路推定システムの画面表示方法の例を示す図。 任意の入出力データのデシベル値が与えた時にそこから求まる伝達率のデシベル値の尤度の確率分布を示す図。 入出力データのデシベル値が与えた時にそこから求まる伝達率のデシベル値の尤度の確率分布を算出する際に、計算速度を向上するための効率化について説明する図。 鉄道車両の車室内騒音に対する伝達経路とその寄与の例を示す図。 図１１に示す伝達経路について、入力信号と伝達率の一次結合として出力信号が求まると仮定した伝達経路モデルで表現した例を示す図。 以降の説明を簡単にするための２入力１出力系伝達経路モデル例を示す図。 図１３に示した２入力１出力系伝達経路モデルについて、従来技術である最小二乗法による線形重回帰分析を用いた伝達経路分析技術による伝達率同定手法を示す図。 従来技術の伝達経路分析技術の１つ目の課題である、伝達率が負値となる例を示す図。 従来技術の伝達経路分析技術の２つ目の課題である、バラツキの影響が稼動状態の負荷の大小に依存することを示す図。 従来技術の伝達経路分析技術の３つ目の課題である、入力信号同士の独立性が低いデータから重回帰分析で伝達率を求める際の不確定性を示す図。

以下、本発明の実施形態について図を用いて説明する。

本発明の実施例１について、図１ａ、図１ｂから図８を用いて説明する。

まず図１ａ、図１ｂを用いて、本発明の振動騒音伝達経路推定システム及び伝達率推定方法により、バラツキをもつ入出力データ１４ａから伝達率をデシベル値のまま確率推論することができることを説明する。

特許文献１に示すような従来技術ではバラツキを持ったデータからノイズを除去してしまうため、伝達率が画一的な数値として得られてしまい、得られた同定結果に対する信頼性の議論ができなかった。しかも意図せず伝達率が負値になるばかりか、入力データの独立性が十分でない場合は、そもそも手法適用の妥当性が低かった。

これに対して、本発明では、バラツキをもつ入出力データ１４ａから伝達率を物理量やパワー値で取り扱わず、デシベル値のまま確率推論する。図１ａは、図１３のデータベースＤＢに蓄積されたパワー値のデータ１３ａのうちｘ_ｍ、ｙ_ｍの関係を表記したグラフであり、ここに示すように、入出力データ１４ａをパワー値表現したときに、各データが多次元空間上の座標として表される入出力データの点群から、その関係性を示す伝達率を確率分布として推定・算出する。ただし，伝達率の確率推論の際，伝達率の物理量やパワー値ではなく，伝達率に関してはデシベル値を推論するように工夫する。このようにして各周波数で得られた伝達率のデシベル値の確率分布は図１ｂに示すように周波数−デシベル表現伝達率特性グラフ上で、ある幅を持った滲んだ曲線として表される。

これにより、パワー値表現の伝達率について非負値探査が可能で、データのバラツキの取り扱いが稼動条件で等価となり、さらに従属性の強いデータからの同定に対してもその不確実性に関する尺度を保持することが可能となる。さらには、得られた伝達率推定結果を、その後で寄与率計算結果に対して用いた場合も、結果が確率論的な取り扱いをすることができる。すなわち次期製品の騒音推定結果についても、その推定結果に対する信頼性の尺度を結果に対する確率として付与することができる。

バラツキを持ったデータからそれを結論として得られる原因について確率推論する手法としては、近年普及しつつある方法としてベイズ統計がある。ベイズ統計ではモデル式の線形性を必要としない。つまり図１３のようにデシベル値で表現されたモデル式についても、推定することが可能である。そこで所望とする伝達率のデシベル値に関する確率推論に対して、ベイズ統計を適用することを考えた。

その実施例について図２から図５を用いて説明する。ベイズ統計はベイズの定理を基礎として、因果関係を形成する原因と結果について、結果として得られたデータをもとに原因の発生確率を確率分布として推定する方法である。

図２はベイズの定理を説明した図で、横軸に仮説Ｈ、縦軸にデータＤを示したグラフである。このグラフで、実線で示す範囲は仮説Ｈが成立する範囲、点線で示す範囲はデータＤが得られる範囲である。また（４）式はベイズの定理を示す式である。

ここでは、仮説Ｈが成立しさらにデータＤが得られる同時確率Ｐ（Ｈ∩Ｄ）が、Ｐ（Ｈ｜Ｄ）Ｐ（Ｄ）とも、Ｐ（Ｄ｜Ｈ）Ｐ（Ｈ）ともどちらでも表せることを利用したものである。この結果、データＤが得られたときの仮設Ｈに関する条件付き確率Ｐ（Ｈ｜Ｄ）について、仮説Ｈが成立していると考える事前確率分布Ｐ（Ｈ）と、仮説Ｈが成立しているときにデータＤが得られる条件付き確率であるＰ（Ｄ｜Ｈ）の積で表せることを示している。なお，この条件付き確率Ｐ（Ｄ｜Ｈ）は尤度とも呼ばれる。

（４）式に示した事前確率分布Ｐ（Ｈ）と尤度Ｐ（Ｄ｜Ｈ）の積として定まる事後確率Ｐ（Ｈ｜Ｄ）の関係式において、新たなデータＤ´が得られた時の新たな事後確率Ｐ（Ｈ｜Ｄ´）を求める式が（５）式である。

ここでは、データＤが得られたときの条件付き確率Ｐ（Ｈ｜Ｄ）、つまり事後確率分布について、これを再び事前確率分布として活用することで、次にデータＤ´が得られたときに、そのデータに対応する尤度Ｐ（Ｄ´｜Ｈ）を再び事前確率分布に乗算することで事後確率分布を更新でき、以下これを逐次繰り返すことによってデータが生み出され得る仮説に関する確率分布を真値に漸近させることができる特徴がある。

なお前記データによる事後確率分布の更新はベイズ更新と呼ばれる。ベイズ更新について図３、および図４でコイントス問題を例により詳しく説明する。図３はベイズ統計を説明する際によく引き合いに出されるコイントス問題を表したものであり、コインを放り投げた時、地面に落ちたコイン表が上か、あるいは裏が上かの情報を逐次得ることにより、コインの表が上になる確率分布の推論結果を徐々に真値に漸近させる問題である。

図４に図３に示したコイントス問題におけるベイズ更新による事後確率分布の変化について具体的な例を挙げて示す。図４において１５１、１５２、１５３、１５４は表が出る確率θと確率密度分布Ｐ（θ）の関係を示すグラフであり、１５１は初期状態、つまり事前情報がない無情報の状態を表しており、このときには確率密度分布Ｐ（θ）は表が出る確率θに対して一定である。また図４において１６１、１６２、１６３はコイントスの結果（表、裏、表の順に発生したものとする）であり、尤度であらわすことができる。ここでは、コイントスの結果１６１、１６２、１６３を受けて表が出る確率θと確率密度分布Ｐ（θ）の関係を示すグラフ１５２、１５３、１５４が順次変化していく様子を示している。

最初、コインの表が出る確率θについて事前情報がない（無情報）であるとすれば、表の出る事前確率分布は一様であると仮定する。これを示したのが図４の１５１である。次にコインを投げた時、表が出たとする。この時の尤度とは、「表が出る確率がθであるコインを投げて表が出る確率」であるから、この確率分布はθであり図４の１６１のようになる。これをもとの事前確率分布１５１に掛け合わせると、１回目表が出た時の事後確率分布がグラフ１５２のように得られる。この事後確率分布を次の推定で事前確率分布として使用する。次にコインを投げたとき、裏が出た１６２とする。この時の尤度は、「表が出る確率がθであるコインを投げて裏が出る（＝表が出なかった）確率」であるから、この確率分布は１−θで１６２のようになる。これを先ほどの事前確率分布１５２に掛け合わせると、２回目裏が出た時の事後確率分布が１５３ように得られる。このようにして得られたデータに対応する尤度を事前確率分布に逐次かけていくことで、データが得られた後の事後確率分布を更新していき、真値に近づけていく。

図５は図４で説明したベイズ更新について、図１３で示した２入力１出力モデルに適用した例を示す図である。

なお、詳細は後述するが、２次元の確率分布で示される伝達率のデシベル値（図では「伝達率レベル」と表現している）Ａ_１、Ａ_２の存在確率分布の表現としてはマルコフチェーンモンテカルロ（Ｍａｒｋｏｖ Ｃｈａｉｎ Ｍｏｎｔｅ Ｃａｒｌｏ：以後ＭＣＭＣと略す）と呼ばれる推定パラメータの離散点群表現が有効である。このため以降の図では確率分布は、このような点群で表現する。図４同様、事前情報が無情報であるとすれば、伝達率レベルＡ_１と伝達率レベルＡ_２の事前確率分布１５１は一様であると仮定するが、上記のように確率密度分布を点群で表現する場合、一様密度の確率分布に従う点群の計算方法としては、ラテン・ハイパーキューブ・サンプリング（Ｌａｔｉｎ Ｈｙｐｅｒｃｕｂｅ Ｓａｍｐｌｉｎｇ：ＬＨＳ）などがある。

また図４に示したコイントスの問題同様、ここにデータベースＤＢからＴｉｍｅＩｎｄｅｘが０のときの入出力データ１４ａをデータセット０（Ｘ_ｍ１ ^（０）、Ｘ_ｍ２ ^（０）、Ｙ_ｍ ^（０））として切り出し、これに対応する尤度１６１を求め、それを先ほど得られた事前確率分布１５１に掛け合わせることでＴｉｍｅＩｎｄｅｘが０のデータを入手した後の事後確率分布１５２を得る。

次にデータベースＤＢからＴｉｍｅＩｎｄｅｘが１のときの入出力データ１４ａをデータセット１（Ｘ_ｍ１ ^（１）、Ｘ_ｍ２ ^（１）、Ｙ_ｍ ^（１））を切り出し、これに対応する尤度１６２を求め、それを先ほど得られた事前確率分布１５２に掛け合わせることでＴｉｍｅＩｎｄｅｘが１のデータを入手した後の事後確率分布１５３を得る。

これを繰り返すことで、データベースＤＢから入出力データ１４ａが得られる度に事後確率分布が逐次更新され、最終的に伝達率レベルＡ_１と伝達率レベルＡ_２の存在確率分布が次第に真値に近づいていくものと推定される。図５では、１５４が最終推論値を示すものとする。なお最終推論においては、真値が伝達レベルＡ_１と伝達率レベルＡ_２の確率分布の期待値として推定される。

ここで尤度Ｐ（Ｄ｜Ｈ）算出の考え方を簡単に説明しておく。まず尤度Ｐ（Ｄ｜Ｈ）とは、仮説Ｈが成立しているときにデータＤが得られる条件付き確率のことである。この尤度は例えば以下のようにして定めることができる。例えば図５においてＴｉｍｅＩｎｄｅｘが０のときのデータセット０（Ｘ_ｍ１ ^（０）、Ｘ_ｍ２ ^（０）、Ｙ_ｍ ^（０））について、入力データＸ_ｍ１ ^（０）、Ｘ_ｍ２ ^（０）と、それぞれに仮定した伝達率レベルＡ_１と伝達率レベルＡ_２から（１）式を実行し、出力データｙ^（０）を推定する。そのうえで、推定した出力データｙ^（０）と計測した出力データＹ_ｍ ^（０）の近似度合いを評価して尤度とする。尤度は伝達率レベルＡ_１と伝達率レベルＡ_２を可変設定しながら繰り返し実行され、複数の伝達率レベルの組み合わせの時の尤度をグラフ表示したものが、図５のグラフ１６１、１６２、１６３などである。

図６は従来手法である最小二乗法による線形重回帰分析と、本発明の伝達率デシベル値のベイズ推定について比較したものである。なお図６では、縦軸にデータが１個、２個、Ｎ個のケースを想定し、横軸に最小二乗法による線形重回帰分析と、伝達率パワー値の確定論的算出と、伝達率デシベル値（図６では伝達率ｄＢ値と表記）のベイズ推論による結果を対比して図示している。

図６において、最初に従来手法である最小二乗法による線形重回帰分析とその結果について説明する。まずデータが１個である図６上段左の表記について、従来の線形重回帰分析で一つのデータ１のみが与えられると、図６上段左のように入出力データｘ_１、ｘ_２、ｙで定義される入出力データ空間上に、このデータ１に対応する点Ｐ１（ｘ_ｍ１ ^（ｎ）、ｘ_ｍ１ ^（ｎ）、ｙ_ｍ ^（ｎ））が得られる。

この時、この点を含む面を形成する場合、伝達率ａ_１とａ_２の間には、ω_１＝ｘ_ｍ１ ^（ｎ）／ｙ_ｍ ^（ｎ）、ω_２＝ｘ_ｍ２ ^（ｎ）／ｙ_ｍ ^（ｎ）として、ａ_１ω_１＋ａ_２ω_２＝１という関係が成立する。すなわち、図６上段中のように伝達率空間［ａ_１、ａ_２］平面上で、一つの直線を描くこととなる。

ここにもう一点データ２が与えられると、点Ｐ１のほかに追加したデータ２に対応する点Ｐ２が得られる。このとき、図６中段左のように原点と点Ｐ１、Ｐ２の３点で形成する回帰平面の傾きａ_１、ａ_２が定まるが、この時、伝達率空間上では図６中段中央のように、ａ_１とａ_２の間の関係を示す直線がデータ１、データ２それぞれに対して１本ずつ、計２本描かれ、この２本の直線の交点が図６中段左の回帰平面の傾きａ_１、ａ_２を示す。

一方、さらにデータが沢山得られた場合、図６下段左のように重回帰分析ないし最小二乗法では無数の点群の間を平均的に通るように回帰平面ＰＬの傾きを定めるが、これは図６下段中央のように各点で定まるａ_１とａ_２の間の関係を示す複数の直線の交点の平均的な位置を求めていることになる。しかし、前述のように高負荷稼働状態で大きなばらつきが得られた場合、図６下段左のように重回帰分析ないし最小二乗法では、回帰平面ＰＬの傾きがこの高負荷稼働状態でのバラツキの影響を大きく受けることとなる。つまり、図６下段中においてこれまでの平均的な交点から大きく外れたところを通る直線が描かれ、これまでの傾向から大きく外れた直線への重みがかなり大きくなる。

これに対して、本発明の方法では、図６上段中のようにデータが一つ得られる度に伝達率パワー値空間［ａ_１、ａ_２］平面上で、一つの直線を描くのではなく、図６上段右のように伝達率デシベル値空間［Ａ_１、Ａ_２］平面上で、尤度に対応する確率分布を持った領域を描くこととなる。なお、図６では確率分布を領域で示しているが、先述のように点群で表してもよい。Ａ_１、Ａ_２の値はデシベル値であるので、この値が負値であったとしても、伝達率は非負のままとなる。

また、さらにデータが得られた場合、図６中段右のように、今得られたデータに対応する伝達率デシベル値空間［Ａ_１、Ａ_２］平面上の分布が得られることとなり、これらが重なった領域の伝達率が、真値の候補として確率が高いことを意味する。これを繰り返すことで図６中段右のようにそれぞれのデータの尤度に対応する伝達率デシベル値空間［Ａ_１、Ａ_２］平面上の領域がもっとも多く重なった部分の確率が高くなる。これにより、前述の伝達率空間［ａ_１、ａ_２］平面上では高稼働状態で大きなばらつきが得られた場合においてその影響を課題に受ける伝達率推論結果が、本発明の手法では単純に一回の尤度の掛け合わせに留まることで、その影響が大きくなることを回避できる。

このような本発明の手法は、図７ないし図８に示すようなシステムで構成することができる。

図７はシステム構成の全体像を示す図で、振動騒音センサー１１と、これらセンサー１１からの信号を周波数変換する周波数変換装置１２と、センサー信号データ１１ａないしはセンサー信号データ１１ａを周波数変換されたデータ１２ａを格納するデータ格納装置１３と、格納されたデータを演算するデータ演算装置１４と、演算された結果を元に、伝達経路推定結果表示する伝達経路推定結果表示装置１５からなる、機器や装置の振動騒音伝達経路を推定する振動騒音伝達経路推定システム９０を構成する。

図７のシステムにおいて、振動騒音センサー１１は、図１１の屋根、荷棚、窓、側、床など内装各部位に設置したセンサー、および図１１の車内騒音の部位に設置したセンサーなどがこれに相当している。またデータ格納装置１３は、図１３のデータベースＤＢを含んでいる。またデータ格納装置１３は、データ格納にあたり、（３）式で説明したモデル式１３ｂを実行する。１３ｂで示した（３）式は、モデルＭを表現するモデル式であり、データ格納装置１３は解析の対象とした２入力１出力モデルＭのモデル式と、パワー値をデシベル値に変換するときの変換式を実行する。

格納されたデータを演算するデータ演算装置１４は、少なくとも以下の処理を実行する。例えば図５の処理を行う場合の典型的な処理ステップの例は以下のようである。ただし、この前提としては、同一時刻に複数のセンサーで検知した騒音に関するデータがデシベル値として、サンプル入手時刻とともにデータベースＤＢに保存されているものとする。また事前確率分布１５１が準備されているものとする。

処理ステップＳ１：データベースＤＢから所定時刻に入手した入出力データ１４ａのデータセットを入手する。

処理ステップＳ２：入手した入出力データ１４ａのデータセットを用いて、このデータセットに対する尤度を算出する。尤度の算出は先述の手法が採用可能である。

処理ステップＳ３：事前確率分布（最初の状態では事前確率分布１５１）と尤度（最初の状態ではグラフ１６１）を用いて事後確率分布を得るとともに、求めた事後確率分布を次回の事前確率分布として扱う。

処理ステップＳ４：求めた事後確率分布が最終推論値に到達していることを確認し、到達していない場合には処理ステップＳ１に戻り、データベースＤＢにおける所定時刻を変更し、連続する時系列の次のデータセットについて、処理ステップＳ１から処理ステップＳ３の処理を繰り返し実行する。

処理ステップＳ５：処理ステップＳ４の処理により、求めた事後確率分布が最終推論値に到達していると判断されるとき、最終推論値として求められた伝達率レベルＡについて、この伝達率レベルＡを（１）式に代入してｙ^（ｎ）を求めることで伝達経路を求める。これにより、図１２右側に示した複数の振動源、音源と、伝達率レベルと、車内騒音ｙ^（ｎ）の関係を示す伝達経路が明らかになる。

処理ステップＳ６：作業員が視認しやすい形態で伝達経路推定結果を伝達経路推定結果表示装置１５に表示するための画像編集を実施する。

図８は伝達経路推定結果表示装置１５における振動騒音伝達経路推定結果１５ａの画面表示例を示したものである。たとえば、図では画面左側に示されたグラフィカルモデルで表現された伝達経路上のある要素を選択すると、その右の画面にその伝達率寄与の不確実性が周波数特性の確率分布として表現されるとともに、左下の画面では不確実性の影響で騒音の周波数特性やＯ．Ａ．などの全体への影響がどの程度かが視覚的に表示される。

これら図７ないし図８で示したシステム構成のうち、データ演算装置１４にてデータを演算して振動騒音伝達経路推定結果１５ａを算出する際に、図１２で示したデシベル値の尤度の確率推定を用いる。

これにより、モデル式が非線形なデシベル値で記述されていても、入出力データ、伝達率ともに物理量やパワーに還元することなく、デシベル値のままで推論および確率分布の更新が議論できる。これにより従来技術における３つの課題「伝達率が意図せず負値になる」「バラツキの取扱が高負荷稼動状態と低負荷稼動状態で等価でない」「バラツキをノイズとして切り捨ててしまうため不確実性ないしは結果の信頼性を議論できない」を解決する技術およびこれを搭載したシステムを提供することができる。

前述のようにベイズ更新の際は、事前確率分布に得られたデータに対応する尤度をかけることで事後確率分布を得るが、入出力データ１４ａの数が増えて、探査空間の次元が大きくなると、軸ごとに一様に探査する手法では価値のある情報を取得するのに計算時間がかかる。

これを解決する手段としてはマルコフチェーンモンテカルロ（Ｍａｒｋｏｖ Ｃｈａｉｎ Ｍｏｎｔｅ Ｃａｒｌｏ：以後ＭＣＭＣと略す）とよばれる推定パラメータの離散点群表現並びに探査方法が近年有効とされつつある。ＭＣＭＣでは、確率分布を点群の密度で示し、メトロポリス法（Ｍｅｔｒｏｐｏｌｉｓ法）のような効率化探査手法と組み合わせることで、反復計算回数を抑えながら効率よく確率分布を表現する点群を求めることが可能となる。

本発明の実施例２では、実施例１に記載の伝達率のデシベル値を取り扱うことを特徴とした事前・事後確率分布および尤度計算に対して、図５の１６１に例示するように、観測された入出力データ１４ａをもとに、尤度を多点サンプルとして表現することを特徴とした方法と、これを搭載したシステムを提供する。

実施例２に記載のＭＣＭＣ法では尤度の多点サンプルを得るために、反復演算を行わなければならないが、図１３に示したモデル式とデシベル値の取り扱いによる特殊性をうまく用いるとさらに演算を効率化できる。

これを説明したのが図９、図１０である。いま仮に、ＴｉｍｅＩｎｄｅｘがｎの時の入出力データ１４ａとしてＤ_ｎ（Ｘ_ｍ１ ^（ｎ）、Ｘ_ｍ２ ^（ｎ）、Ｙ_ｍ ^（ｎ））を得た時、この入出力データ１４ａに対する多点サンプル表現尤度として図９のような点群を得たとする。

また、次のＴｉｍｅＩｎｄｅｘであるｎ＋１の時の入出力データ１４ａとしてＤ_ｎ＋１（Ｘ_ｍ１ ^{（ｎ＋１）}、Ｘ_ｍ２ ^{（ｎ＋１）}、Ｙ_ｍ ^{（ｎ＋１）}）を得たとし、先ほどのＴｉｍｅＩｎｄｅｘがｎの時の入出力データ１４ａと、Ｘ_ｍ１ ^{（ｎ＋１）}＝Ｘ_ｍ１ ^（ｎ）＋ΔＸ１、Ｘ_ｍ２ ^{（ｎ＋１）}＝Ｘ_ｍ２ ^（ｎ）＋ΔＸ２、Ｙ_ｍ ^{（ｎ＋１）}＝Ｙ_ｍ ^（ｎ）＋ΔＹなる関係があったとする。つまり、ＴｉｍｅＩｎｄｅｘがｎの時の入出力データ１４ａに対してΔＸ１、ΔＸ２、ΔＹ、だけの変化があったとする。ここで図１３に示した入出力データ１４ａのデシベル値（Ｘ_１、Ｘ_２、Ｙ）と伝達率尤度のデシベル値（Ａ_１、Ａ_２）、それにモデル式１３ｂを考える。

（６）式はこの場合のモデル式１３ｂの例を示している。（３）式のモデル式にさらにデシベル値（Ａ_１、Ａ_２）のモデル式が加味されている。（６）式の関係から、さらに（７）式を導くことができる。

（８）式は、（７）式において、連続するＴｉｍｅＩｎｄｅｘであるｎとｎ＋１における各値Ｘ_１、Ｘ_２、Ｙ、Ａ_１、Ａ_２の関係を示している。また図１０は（８）式の関係を模式的に示した図である。

これによれば、ＴｉｍｅＩｎｄｅｘがｎの時の入出力データ１４ａに対する多点サンプル表現尤度（Ａ_１ ^（ｎ）、Ａ_２ ^（ｎ））に対して、ＴｉｍｅＩｎｄｅｘがｎ＋１の時の入出力データ１４ａに対する多点サンプル表現尤度（Ａ_１ ^{（ｎ＋１）}、Ａ_２ ^{（ｎ＋１）}）は、Ａ_１ ^（ｎ）を−ΔＸ１＋ΔＹ、Ａ_２ ^（ｎ）を−ΔＸ２＋ΔＹだけ移動させたものとなり、つまり（Ａ_１ ^{（ｎ＋１）}、Ａ_２ ^{（ｎ＋１）}）＝（Ａ_１ ^（ｎ）−ΔＸ１＋ΔＹ、Ａ_２ ^（ｎ）−ΔＸ２＋ΔＹ）となる。この特性を利用すれば、データを得られる毎に、尤度をＭＣＭＣ法を用いて反復計算をする必要はなく、計算が効率化できる。

１１：振動騒音センサー
１１ａ：センサー信号データ
１２：周波数変換装置
１２ａ：周波数変換信号データ
１３：データ格納装置
ＤＢ：データベース
１３ｂ：モデル式
１４：データ演算装置
１４ａ：入出力データ
１６１、１６２、１６３：尤度
１５１、１５２、１５３、１５４：確率分布
１５：伝達経路推定結果表示装置
１５ａ：振動騒音伝達経路推定結果
９０：振動騒音伝達経路推定装置

Claims (5)

1. 複数の振動騒音センサーからの信号を格納する格納装置と、前記格納装置に格納されたデータを用いて振動騒音を演算する演算装置と、前記演算装置の演算結果を元に、振動騒音の伝達経路推定結果を表示する表示装置からなる振動騒音伝達経路推定装置であって、
   前記格納装置は、複数の前記振動騒音センサーからの信号をデシベル値として測定時刻ごとに記憶し、
   前記演算装置は、同一測定時刻における複数の前記振動騒音センサーからの信号を入出力データとし、入出力データと伝達率を用いて振動騒音現象を表現するモデル式に適用し、同一測定時刻の入出力データに対応する伝達率のデシベル値に対する尤度を計算し、計算した尤度を事前に想定している伝達率の事前確率分布に乗算することで伝達率の推定値の事後確率分布を得、前記同一測定時刻を逐次変更して前記事後確率分布を繰り返し求めることで前記振動騒音の伝達経路を推定することを特徴とする振動騒音伝達経路推定装置。
2. 請求項１に記載の振動騒音伝達経路推定装置であって、
   前記入出力データをもとにした尤度の算出の際に、反復演算を基にした多点サンプルによる表現を用いることを特徴と振動騒音伝達経路推定装置。
3. 請求項２に記載の振動騒音伝達経路推定装置であって、
   入出力データをもとにした尤度算出の際に、一度、基準入出力データに対して反復演算を基にして得られた多点サンプル表現基準尤度多点サンプルを算出した後、これをデータとしてメモリ上にストアしておき、逐次別の入出力データが得られた際に、多点サンプル表現基準尤度格納用メモリ上にストアされた多点サンプル表現基準尤度を所定の演算方法にて加減算することで他の入出力データに対する多点サンプル表現尤度の計算を効率化させ計算速度を向上させたことを特徴とする振動騒音伝達経路推定装置。
4. 請求項１に記載の振動騒音伝達経路推定装置であって、
   入出力データと伝達率を用いて振動騒音現象を表現する前記モデル式において、出力データは入力データと伝達率の積和として定義され、仮設定された伝達率と計測した入力データの積和として求めた仮の出力データと計測した出力データの近似の度合いにより前記尤度を定めることを特徴とする振動騒音伝達経路推定装置。
5. 複数の振動騒音センサーからの信号をデシベル値として測定時刻ごとに記憶し、
   同一測定時刻における複数の前記振動騒音センサーからの信号を入出力データとし、入出力データと伝達率を用いて振動騒音現象を表現するモデル式に適用し、同一測定時刻の入出力データに対応する伝達率のデシベル値に対する尤度を計算し、計算した尤度を事前に想定している伝達率の事前確率分布に乗算することで伝達率の推定値の事後確率分布を得ることを特徴とする伝達率推定方法。

Images