周波数フィルタを用いた超解像計測の検証と水中気泡呼吸モードの界面変位測定への応用

全文

(1)可視化情報学会論文集 Vol.41, No.3 (2021 年 3 月）pp.1-10. 周波数フィルタを用いた超解像計測の検証と水中気泡呼吸モードの界面変位測定への応用＊奥隆夫 1), 平原裕行 2), 姜東赫 2). Validation of a super-resolution measurement using frequency filtering technique and optical displacement analysis of bubble surface in breathing mode oscillation Takao OKU, Hiroyuki HIRAHARA, Donghyuk KANG ABSTRACT A sound is generated when a bubble detaches from a nozzle installed in the water. The sound emission is known as a monopole sound source. This frequency is expressed by Minnaert's equation, and many studies have been conducted on bubble deformation, detachment, and coalescence. However, many physical processes still involve the unsettled phenomena regarding the mechanical behavior of the sound source. In this study, we report that the breathing mode was directly analyzed by an optical measurement, following the simulation of the breathing mode by CFD in the previous report. The present superresolution measurement gives a sufficient result to discuss the oscillatory motion of bubble surface where the displacement of the bubble is so small that the normal observation cannot capture. In the present research, the numerical simulation of the measurement method was conducted, and measurement of the vibration disc were performed as the preliminary experiments. Thus, the optimum measurement parameters were extracted, and they were applied to the displacement measurement of the breathing mode of the bubble. The results show a good agreement with the CFD ones. Keywords :. Bubble dynamics, Optical measurement, Super-resolution, Breathing mode. 議論されている．Minnaert バブルにおいては，周囲の水によ. １．緒論. ってバブル放射音圧が変化する問題や高次変形モードが呼. 水中に置かれたノズルから気泡が離脱する際には，気泡のもつ固有振動数による単極子振動が発生する．このような水中気泡の挙動は Rayleigh-Plesset 方程式 1) と呼ばれる運動方程式によって解析がなされており，これまでに基本的な運動のみならず複雑な変形を伴う気泡の運動および振動も含めて多くの研究がなされてきた．水中バブルの振動挙動に関してはバブル内の気体が，水蒸気である場合と空気の場合とに分けられる．前者はキャビテーションに関連して非常に多くの解析がなされており，バブルの振動に関する解析や実験が行われ多くの成果が挙げられているのは言うまでもない．一方，後者は Minnaert バブルと呼ばれており 2)，これを対象とした研究では，音の発生源としてのバブルの振動の力学が * 原稿受付 2020 年 12 月 9 日 1) 非会員埼玉大学理工学研究科 (〒338-8570 埼玉県さいたま市桜区下大久保 255, E-mail:[email protected]) 2) 正会員埼玉大学理工学研究科. 吸モード振幅に及ぼす影響などが議論されているが，いずれも不明な問題が残されている．著者らも種々のノズル口径に対して高速度カメラによる可視化実験を行い 3),4)，異なる気泡直径の気泡離脱時の放射音を計測した．カメラの画像取得周波数は 2000 fps であり，撮影対象の観測解像度は 40μm/pixel である．一例として Fig.1 に実験によって得られたバブルの初期の変形過程とバブルから 21 ㎜離れた位置に設置された水中マイクロフォンによって得られた圧力波形の関係を示す．得られた圧力波形をもとにして，Rayleigh-Plesset 方程式から計算される振動振幅値および VOF 法を用いた CFD シミュレーションを行い 4)，圧力変動と気泡体積変化を比較し，半径の変動振幅はほ. ぼ数 μm のオーダであることを確認した．しかしながら，このような可視化観察でバブル界面の変位を直接観察するには限界があり，呼吸モードにおける気泡界面の振動振幅を従来の手法で計測することは困難である．そこで，本研究ではこの気泡の呼吸モードの振幅を計測. ─1─.

(2) 奥隆夫，平原裕行，姜東赫するために，画像処理によりカメラの解像度以下であるサブピクセルオーダの変位を抽出する超解像解析による手法を. ２．画像処理手法の概要. 応用して呼吸モードの変位計測を行った．水中のノズルから. バブル離脱時の呼吸モードの例のように，気液二相流にお. 離脱した気泡振動は呼吸モード以外の低周波数、高次変形モ. いては気液界面が大きく変動するとともに種々の固有振動. ードが重なった界面運動を行うことが報告されている．. が随伴されることがよく見受けられる．このような場合には，. Fig.1(b)には実験により得られた圧力波形を FFT 解析した結. 可視化観察の解像度の制約から微細な振動の観察が困難と. 果を示した．最も大きなピークである 1016Hz がこの時の呼. なる．このような中において，撮影された画像から種々の情. 吸モードの共振周波数である．350Hz の低い周波数は実験容. 報を読み取ることを目的として解像度の低いデータ画像か. 器内で発生する共鳴およびノズル内の気柱共鳴との干渉に. ら高周波数成分を復元し、表現する手法として超解像 5)と呼. よって発生する変動成分４）であるが，現時点ではその騒音. ばれる手法が種々考案されてきた．この種の画像処理では，. 源は特定できていない．また，バブル変形の高次変形モード. 平行変位した画像情報を画素間で直接内挿して求める手法. は呼吸モードが支配的となる音の発生には関係していない．. 6)-8). 本研究では，呼吸モードの基本周波数を抽出し，その周波数. できる．また，画像強度の変化から撮影対象の周波数情報を. における振動変位を定量化した．実測に先立ち，計測手法の. 抽出し，解像度以下の変化に対しては特定周波数での振動を. 妥当性を検討するため，撮影周波数，画像解像度，輝度階調，. 選択し増幅することで，目視により確認することが不可能で. および振動する界面の画像界面幅をパラメータとして，数値. あった変動を画像で可視化する手法も提案されている. 実験と並行して動電型加振器ディスク境界面を対象として. 15). 校正実験を行った．これらの予備検討より，撮影解像度以下. flow の解析では非常に多くの画像解析手法が考案されてい. の振動振幅に対して，超解像解析が可能であることが示され. るのは言うまでもない 16),17)．Optical flow は対象物体の変位. た．これらの確認実験の条件を用いて，内径 4 および 6mm. による強度分布は，その関数形状を維持しつつ平行移動する. のノズルから離脱する気泡の呼吸モードの振動振幅を計測. という仮定に基づき変位量を算出する手法である．本研究で. したのでその結果について報告する．. は本手法に基づいて水中気泡の呼吸モードの撮影解像度以. と位相相関を用いて求める方法 9)-11)とに大別することが. 12)-. ．また一般に物体の変位抽出を目的としている Optical. 下の変位を抽出することを試みた．. 2.1 周波数フィルタを用いた Optical flow 変位計測初めに，平行移動する移動物体の画像に対する微小な 2 時刻間の Optical flow について 1 次元空間で考える．対象とする物体輪郭関数を I (x ,t ) とする．ここで，恣意的に輪郭関数と言う言葉を用いたのは I (x ,t ) に可視化で得られた実像のみを用いることを意味しないためであり，これについては後述する．時刻 t 0 において得られた画像関数を用いて，平行移動する移動物体の輪郭関数を f (x ) とすれば，次のように定義できる．. I (x ,t0 )  f (x ). (1). 時刻 t1 においては. (a) Bubble deformation and pressure signal.. I (x ,t f (x  x ) 1). (2). となる．ここに  x は t 0 ～ t1 間に物体が移動した距離である．  x が十分に小さければ，I (x ,t1) に対する 1 次の Taylor 展開. によって次の関係が成り立つ．. I (x ,t f (x  x )  f (x )  f (x )x 1). (3). これにより，変位量  x は強度分布における空間方向の１次. (b) FFT spectrum of pressure signal. Fig.1 Bubble detachment and following deformation process. Nozzle diameter φ=6mm.. 勾配を用いて，式(4)のように推定することができる． f (x  x )  f (x ) (4) x   f (x ) このように，強度の分布を維持し変位する場合，任意の座標 x における変位はその位置における強度の変化が既知であ. れば Taylor 展開による近似を考えることによりその変位を. ─2─.

(3) 周波数フィルタを用いた超解像計測の検証と水中気泡呼吸モードの界面変位測定への応用推定することができる．今回の計測では，特定の周波数に対する振動成分の変位計. 3.1 数値試行実験 3.1.1 数値試行実験の概要. 測を行うことが目的である．変位の推定には前述の Optical. 対象とする画像の原関数の形状を次式で定義する．. flow 解析を基本とするが，これに時間周波数フィルタを組 12).  j (x ). 1 2. (erf(cx )  1). (14). を参考として，本研究. ここに erf( ) は誤差関数である．今回の計測においては，ピ. では変位を定量化する手法を検討して解析に応用した．以下. クセル以下の解析精度が重要であるために，画像計測に則し. にその手順を概説する．. た輝度分布を適切に表現する必要がある．そこで，画像上に. み合わせて画像強調を行った手法.  (t ) と表すことにす今，変位  x を時間の関数として x  る．これにより式(2)は次のように表される．. I (x , t ) f (x   (t )). おける各画素の強度を 1 pixel 当たりの透過光の積分で計算して離散化することにする．i 番目の画素座標を xi ，画素間. (5). 隔を  x とするとこの点の輝度は各時刻に対して，次式で与. ここで時刻を差分的に表して tn  t0  nt とすれば，式(5)を. えられる．. 書き直して，.  J i (x i ,t ). I (x ,t f (x   (tn )) n). j (x   (t ))dx. (15). ここで，式(14)における定数 c は実像における焦点の度合，. 上式に対して Taylor 展開を適用すれば次式となる．. 撮影倍率などによって変化する界面の滑らかさを変化させ. (7). るための係数である．Fig.2 に初期状態における関数形状の.  (tn ) の 2 次以上の項を無視して次式を得る． I (x ,tn ) f (x   (tn ))  f (x )  f (x ) (tn ). . xi x / 2. (6). I (x , tn ) f (x   (tn ))  f (x )  f (x ) (tn )  O ( (tn ) 2 ). x i  x / 2. 例を示す． c =0.3, 0.7, および 50 における界面幅(以下  と. (8). 次に  (tn ) の時間変動について考える．フーリエ変換，逆フ. する)はそれぞれ 10, 4, および 2 pixels である．. ーリエ変換を一般的な形式で以下のように表すことにする．.  G ( ). g (t )e dt , g (t ) G ( )e  2 it. 2 it. d. (9). 今回の計測では，  (tn ) を周波数分解し，バンドパスフィルタを施して，特定周波数の変位量を求めることが目的である．ここで，特定の周波数 m に対するバンドパスフィルタリング操作の演算子を m として G ( ) に対して行うと，以下のようになる．. gˆm (t ) . 2 it 2 it  m  g (t )e dt e d. Fig.2 Original profile of model interface.. (10). このようにして撮影界面を模擬した後に，次式.  (t )  a sin(2Tt ). ここに“^”はバンドパスフィルタリング後の関数を表す．式 (8)にこの変換を施すと次式となる．  f (x )ˆm (tn ) Iˆm (x , tn )  fˆm (x )  f (x )ˆm (tn ) . よって次式により特定周波数の変位が得られる． Iˆ (x ,t ) ˆm (t )   m f (x ). (16). により測定対象周波数  T および振動振幅 a で界面を振動さ (11). せ，界面幅，サンプリング周波数  S および輝度階調について数値実験による評価を行った．解析の対象とする輪郭関数としては，. (12). I (x ,t )  J (x ,t ). (17). 今回の報告では界面変動を 1 次元方向に解析して考察する. として，原画像関数のラプラシアン（2 階微分値）を用いた．. が，この方法は多次元問題に容易に拡張可能である．また，. これは，一般的に知られているように 2 階微分値をとるこ. すべての周波数を含んだ実変位は  (t )  ˆ (t ). とによって画像輝度の最短空間周波数の変化を表現するこ. . m. とができるためである．. (13). 数値実験の手順は以下のとおりである．式(15)で表される. m. によって求めることができる．. 原画像関数 J (x ,t ) のラプラシアンを輪郭関数 I (x ,t ) とする．. ３．本画像処理法の精度検証. 界面変動を検出する位置は，初期状態に設定した関数. 本手法の妥当性を検証するために，ヘヴィサイド関数状の. f (x )  I (x ,t0 ) の界面解析領域において f  が最大となる点. 界面が光学観察によってやや滑らかとなった形状を誤差関. とした． f  は 1 次の中心差分によって求めている．この検. 数で表した数値実験と動電型加振器ディスクを使った実測. 出位置において式(12)による変位検出を行い，得られた変位. による実験の二つを行った．以下にそれぞれについて説明す. 時系列データの最大変位を算出した．界面幅 10pixels( c =0.3). る．. の例を Fig.3 に示す．. ─3─.

(4) 奥隆夫，平原裕行，姜東赫. Fig.3 Representation of the interface. c =0.3.. t ) j (x   (t )) (a) J (x ,. 3.1.2 数値実験の処理手順本手法の解析過程の一例を Fig.4 に示す．図に示した数値実験の条件は， c =0.3, 振幅 a =0.5 pixels, 測定対象周波数に対するサンプリング周波数の比は  S /  T = 100 である．図にはサンプリング開始から 3 時刻の波形を示した． c =0.3 は.  =10 pixels でやや広めの界面状態である．本解析では 0.01pixels 程度までの計測分解能を目指しているためにさまざまな振動変位条件で数値試行実験を行ったがその結果は Fig.5 において述べる．解析の手順は以下の通りである． (i)：原画像の取得（作成），（Fi.4(a)参照） (b) I (x ,t )  J (x ,t ). (ii)：原画像からラプラシアン画像へ変換し I (x ,t ) を決定（Fi.4(b)参照） (iii)： f  の最大となる点を決定 (iv): I (x ,t ) の時間軸方向データをサンプリングし FFT 分析後，抽出する周波数を特定 (v)：特性周波数のフィルタリング (vi)：逆 FFT により時間軸へ変換 (vii)：式(12)により変位を決定以上で変位を求めることができる．Fig.4 の各図は上のプロ. (c) Iˆm (x ,t ). セスの中で(i)⇒(a), (ii)⇒(b), (vi)の各時刻の空間分布波形⇒ (c), (vi)⇒(d)として対応した波形のグラフが描かれている． Fig.4(d)に示した結果について説明する．図は解析によって得られた各時刻の変位量の空間分布が示されている．図の点線とその右端に示された値が条件で与えた各時刻の変位量である．全体にわたって良好な変位の計測値が得られているが，x=±2.3 近傍では計測値の発散が見られる．これは式(12) における分母の f (x ) が 0 に漸近するためである．このような誤差は，本解析においては輪郭関数の分布が平坦な場所で見られる影響であるので適切なカットオフ値を用いて処理. (d) ˆm (t ). することで回避することができる．本解析ではその値を. Fig.4 Data processing and result of numerical inspection. c =0.3, a =0.5 pixels, and  T =10 Hz.. 1.0×10-5 として処理を行った．. 3.1.3 界面幅の影響誤差評価の結果を Fig.5 に示す．示された誤差は，界面に. 果を Fig.5(a)に示す．横軸は振幅 a であり  の値は 1,5,10 で. おける時系列変化の変位を FFT 解析し計測した振幅𝑎𝑎�と与. 16bit である．図より，界面幅が広いほど誤差が小さく，変. えた振動振幅 a と差で示した．初めに  の影響についての結. ある．サンプリング条件は  S /  T =100 であり，輝度階調は. ─4─.

(5) 周波数フィルタを用いた超解像計測の検証と水中気泡呼吸モードの界面変位測定への応用. 0.16. 0.1. 0.8. θ=5. Error amplitude [pixel]. Error amplitude [pixel]. θ=1. θ=10 0.6 0.4 0.2 0. 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 0.08 0.06. 0.14.  S /  T =160. 0.12.  S /  T =200. 0.1. 0.04. 0. 1. 8bit 12bit 16bit. 0.08 0.06 0.04. 0.02. 0.02 0. 0. 0.1. 0.2. 0.3. True amplitude [pixel]. (a) Influence of interface width.  S / T =100..  S /  T =100. Error [pixel]. 1. 0.4. 0.5. 0.6. 0.7. 0.8. True amplitude [pixel]. 0.9. 1. 0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.6. 0.7. 0.8. 0.9. 1. True amplitude [pixel]. (b) Influence of the sampling rate at (c) Influence of the pixel resolution at interface width  =10. interface width  =10 and  S / T =100. Fig.5 Error estimation by numerical inspection.. 動振幅が小さい方が推定の精度が高い結果となった．これは. ディスクの上下動を高速度カメラ(nac HX-7s)により LED 光.  S / T = 160, 200 についても同様な結果であった．これらの. 源(Sumita, LS-L150)を物体背面から照明し透過光を撮影し. 結果は式(12)で示した界面での 1 階微分値 f  は界面幅が広. た．加振器ディスクの変位は 80mm 上方に設置したレーザ. い方が，精度が良く，振幅は小さい場合では Taylor 展開による. 変位計(Keyence LK-080)によって計測した．レーザ変位計の. 1 次近似で十分であることを示唆している．なお，ここでは示. サンプリング周波数は 1kHz とした．変位計の分解能は 3μm. していないが，  >10 については測定精度がより向上するこ. であり,測定された加振器ディスクの変位は A/D コンバータ. とを確認しているが，本研究においては画像計測においては. (NI,USB-6251)により PC に取り込み処理した．加振器ディス. 界面を含めて可能な範囲で広範囲を観察することが望まし. クを 10Hz で上下動させ，高速度カメラのサンプリング周波. いとの観点から  =10 程度が妥当な値であろうと判断し，以下  =10 を基準として解析を進めることとする．. 数を 1kHz，撮影倍率を 147μm/pixel，撮影の輝度階調を 12bit とした．加振振幅は 0.1～1pixels の間で 0.1pixels ごとに変化. 3.1.4 サンプリング周波数の影響. させた．また、撮影では加振器ディスクの界面幅を変化させ.  =10 の時，10Hz で界面を振動させたとき，  S / T. るために，焦点位置での撮影条件を Case-1 とし，後方に 1mm. =100,160,200 と変化させた場合に，与えた振動振幅 a と測定. および 2mm の位置に焦点をずらしたボケ像の撮影条件をそ. 振幅𝑎𝑎�の誤差を評価した結果を Fig.5(b)に示す．この  S / T. れぞれ Case-2 および Case-3 とした．良好なボケ像を得るた. の範囲では誤差は振幅に線形的に比例して増大し a =1pixel. めに今回の実験では Leica のマイクロスコープ Z16 を用い. の時に誤差は 0.06pixels であった．別途行った  =1 および 5. ている．以上の条件により撮影を行い，画像処理から得られ. についても結果は同様であり，変動振幅に対して誤差は 1/10. た変位をレーザ変位計で計測された変位と比較することに. 程度となる結果が得られた．. より，画像処理の精度を検証した．実験は単一周波数および. 3.1.5 輝度階調の影響. 複数周波数の合成波形の 2 種類について実施した．. 次に  =10， S / T = 100 とし，輝度階調を 8，12，16bit と変化させた場合に，振動振幅 a と測定振幅 a との誤差を評価した結果を Fig.5(c)に示す．この結果から，8bit の階調では誤差の変動がみられるが，12 および 16bit では，誤差は同程度であると言える．ここで，それぞれの階調での分解能は 3.9×10-3，2.4×10-4，および 1.5×10-5 となる．このことと Fig.5(c) で検証した結果から，本研究が目的としているサブピクセルでの解析を目的とした場合には，  を 10 ( c  0.3)程度とし，.  S / T >100，輝度階調を 12bit 以上とすれば界面の変位に対. Fig.6 Image acquisition of oscillating head for calibration on the present image analysis.. して誤差は 1/10 程度の計測が十分に行えると言える．. 3.2 加振ディスクの変位計測による精度検証 3.2.1 精度検証実験の概要. 3.2.2 単一周波数による検証実験の解析および結果高速度カメラにより撮影した加振器ディスクの画像を撮. この節では，これまでに説明した手法を用いて実際の加振体を計測して精度検証を行ったのでそれについて記述する．. 影条件 Case-1, 2, 3 について Fig.7 に示す．ディスクを含む上. Fig.6 に動電型加振器を使った検証実験装置の概略図を示す．. 下を 240×64 pixels で撮影し，取得画像から界面近傍の. ファンクションジェネレータ(Gwintek GAG 810)により発生. 60×20pixels を示した．なお，実験でのディスク振動方向は. された電圧信号をアンプ(Emic, 371-A)で増幅し，防振台上に. 鉛直であるが，図には振動方向を水平方向として示している．. 設置した加振器(Emic Shaker, 511-8)を駆動した．この加振器. 画像左側の黒色部分が加振器ディスクである．本画像における加振器ディスクの振動方向の輝度分布を 12bit で抽出し，. ─5─.

(6) 奥隆夫，平原裕行，姜東赫式(15)に対して最小二乗法による曲線近似を行い，係数 c を. れたが，0.1～0.5 pixel 程度の振動振幅であれば，推定された. 算出して，この値より Case-1～3 における界面幅を評価した．. 振幅のばらつきも小さく，再現性の精度も高い．しかしなが. Fig.8 に各撮影条件における界面近傍の画像強度の分布およ. ら，数値実験同様に，変動振幅が大きくなると推定精度が悪. び係数 c の値，界面幅の値を Table 1 に示す．これらの結果. 化し，バラツキも大きくなることも含め，数値実験で確認さ. から，光学観察では界面幅は 7～11 pixels となり，3.1 節で. れた本手法による変位計測の特徴を確認することができた．. 得られた画像処理の精度が良い  =10 近傍であることを確. 誤差の主たる原因の一つは f  の評価が考えられる．これに. 認した．. は，前述のように中心差分を用いた．今回は界面強度の平滑化などの処理は一切用いていない．本来，実像がサブピクセ. (a) Case-1. ル程度となる撮影倍率のものを，界面幅が 10pixels 程度にぼ. (In-focus). かして捉えられているので f  は精度よく評価できるように. (b) Case-2. 配慮したものの，画像ノイズはそのまま結果に反映されてい. (Out-of- focus). るものと考えられる．. (c) Case-3 (Out-of- focus) Fig.7 In- and out-of-focus images of oscillating disc interface. The image size is 60×20 pixels.. Fig.9 Comparison between laser measurement and image processing measurement.. 3.2.3 複数周波数による検証実験の解析および結果ここではいくつかの周波数変動を含む場合についてバンドパス解析をそれぞれに対して施すことによって変位が適切に回復されることを示す．そのために，加振器ディスクに. Fig.8 Intensity distribution of the oscillating diaphragm.. Case. Table 1 Captured interface properties Corresponding Approximated c-value in Eq.(14) interface width  (pixels). Case-1. 7. 0.486. Case-2. 9. 0.325. Case-3. 11. 0.247. 以上の Case-1～3 に対し以下の手順で検証実験を行った． Case-1～3 のそれぞれの撮影条件において，ディスクの振幅. を 0～150μm の範囲において 10 条件で変化させ，それぞれ. (a) Disc displacement. 5 秒間(5,000 フレーム)撮影した．振幅 150μm はカメラの解. 像度換算でほぼ 1pixel 相当する．得られた画像の中央付近の 5 行を選び、得られた原画像関数をもとにして 3.1.2 節で述. べた処理(i)～(vii)を実行して変位を求めた．画像処理から得られた時系列変位の時間平均 opt とレーザ変位計により実測された変位の時間平均値 las を比較した結果を Fig.9 に示す．いずれの平均値も計測 1 秒間における 10 周期の peak to peak の 1/2 で評価した．この図より，. las を 0.1～1 pixel と変化させたときに対して，本手法による画像処理によって良好な振動振幅の計測ができていることが確認された．las が大きい時にはやや大きな誤差が見ら. ─6─. (b) Displacement for each oscillation mode. Fig.10 Displacement analysis for pulse excitation..

(7) 周波数フィルタを用いた超解像計測の検証と水中気泡呼吸モードの界面変位測定への応用正負の矩形波を入力して振動させ、矩形波に対する応答を計. 像度は過去に行った解析において界面振幅が数 mm と予測. 測した．前述の手順により高速度カメラを用いて撮影し，画. されることから Table 2 に示した空間解像度で撮影を行った．. 像を本解析手法で処理した．応答によって発生したディスク. また，前述の予備検証をもとにして，輝度階調を 12bit，界. の変位はそれぞれの離散周波数でバンドパス処理して，その. 面幅  を 10pixel となるように調整した．計測は各ノズルサ. 後，これらの変動を合成して，レーザ変位計で計測した波形. イズの条件について 6 回行った．. と比較した．これらの結果を Fig.10 に示す．矩形波の周波数は 10Hz である．したがって応答周波数はその奇数倍の周波数が出現する．よって画像解析では基本周波数を含めて三つの周波数，10,30,および 50Hz に対して処理を施した．フィルタの帯域はそれぞれ 10±1Hz，30±2Hz，および 50±1Hz である．Fig.10(a)の実線はレーザ変位計で得られた変位を，計測結果との比較のために画像計測条件での pixel に変換して示した．Fig.10(b)は画像計測から求められた. 各周波数での変位の時間変動を示した．得られたこれら三つの信号波形を線形的に合成したものが Fig.10(a)に破線で示した曲線である．実線と破線は良い精度で一致していることが確認でき，本解析手法によって固有の周波数ごとの解析が Fig.11 Experimental apparatus for bubble release. 有効であることが示された．. 4．．離脱気泡の呼吸モードの可視化計測. Table 2 Image acquisition conditions of bubble detachment Acquisition rate Image resolution Nozzle diameter  mm  S Hz mm/pixel. 前章で記述した検討結果をもとにして，水中のノズルから離脱するバブルが呼吸モードで振動する際のバブル界面の. 4. 振幅を計測した結果を以下に示す．Fig.11 がバブルの呼吸モ. 6. 100k. 5.08 9.85. ード計測のための実験装置である．水槽は横 400mm 奥行き 200mm 高さ 600mm で厚さ 2mm のアクリル製である．深さ. 全体画像，界面部拡大画像(320×88 pixels)の画像の例と圧. 500ｍｍの水中の水深 250mm の位置にガラス製ノズルが設. 力波形，界面部分の強度分布，および計測によって得られた. 置されている．ノズルの内径  は 4 および 6mm，長さは. 呼吸モード界面振動振幅の時系列変位を Fig.12 に示す．時. 20mm であり，外部のシリンジポンプから空気を供給した．. 刻原点は，気泡離脱直後の圧力波形の最初の最小負圧となる. 空気流量は 0.186ml/s とし，1 秒間に 1～2 個の気泡が発生す. 時刻を選んだ．ここでは詳細は記載しないが，既報 4)と同様. るように流量を設定した．離脱の挙動は 2 台の高速度カメ. に，本実験においても，水中マイクの圧力波形を FFT 解析. ラで撮影した．気泡界面の局部拡大画像撮影では 1 台は気. し，第 1 主ピークの周波数が，測定された気泡の等価直径か. 泡全体を，2 台目は気泡界面を拡大して同時に撮影している．. ら Minnaert の式によって計算される呼吸モード周波数であ. 全体像は CAM-2 (IDT , XS-5)に中焦点レンズ(Nikon, f=85mm). ることを確認している．各サイズのノズルで測定された圧力. で行い，界面の撮影は CAM-1 (nac HX-7s)にマクロレンズ. 波形には低周波数のピークも含まれるが，これらの低周波数. (Leica, Z16APO) を装着し，対向する位置から LED 光源. モードが画像計測に与える影響は，本解析においてはバンド. （Sumita, LS-L150)を用いて透過光によりバブルのシャドウ. パスフィルタによって除去されている．呼吸モードの共振周. を撮影した．同時に，ノズル出口から水平方向に 15mm，垂. 波数は， =4 および 6mm に対してそれぞれ，1163±52Hz お. 直方向に 15mm 離れた位置に水中マイクロフォン (B&K,. よび 1044±45Hz であった．このばらつきは，離脱時の液滴. 8103, 感度 26μV/Pa)を設置して離脱時の圧力波形をカメラ. の等価直径のばらつきによるものである．次に可視化画像計測による変位計測の結果について述べ. 撮影開始と同期して記録した．水中マイクのサンプリングレ. る．前章で説明した手続きにより，Fig.12(a)で示されている. ートは 40kHz であり，気泡全体像の撮影は 2kHz である．. 界面局部拡大画像の中心近傍の 30 行に対して各行毎に 3.1.2 節の手順(i)～(vii)を実行して  (t ) を求めた．この時，バンド. 4.1 気泡界面を拡大した画像解析および結果気泡界面拡大像の撮影条件を Table 2 に示す．各ノズル径. パスの周波数は呼吸モード固有周波数±10Hz 程度とした．. の条件において，発生した気泡サイズはノズル内径 4 およ. 以上によって求めた結果を Fig.12(d)に示す．図には画像解析. び 6mm のいずれの場合においてもおよそ直径 6~7mm であ. から得られた気泡界面の振動波形を橙色、マイク位置での圧. った．撮影周波数は発生した気泡の等価直径から求められる. 力波形を青色で重ねてプロットしている．変位は圧力波形と. 呼吸モード周波数の 100 倍程度である 100kHz とし、空間解. ほぼ同じ周波数で変化していることがわかる．変位計測にお. ─7─.

(8) 奥隆夫，平原裕行，姜東赫いてはバンドパスを行っているにもかかわらず変位曲線に. 致するはずの圧力の最大値と界面振動波形の最小値の位置. は高い周波数変動成分が見られる．これはバンドパスにおい. は 1/4 周期弱ずれていることが判明した．一方，バブル界面. て不適切な処理が行われているのではなく，処理手順(vii)に. と圧力計測位置の距離がずれていることによって生じる音. おける式 (12) に用いる輪郭関数の勾配が変動しているため. 波の伝播による時間遅れはおよそ 15ms であるので，周期に. に生じたものである．既報 4). において行われた数値解析の. 換算すれば 1/100 周期程度のずれしかない．加えて，既報 4). 結果に見られるように，圧力波形と振幅変位波形とは位相は. で行った数値計算の結果をもとにして考えてみると，数値計. 1/2 周期ずれることを示しており，この関係は Rayleigh-. 算においては当然ながら粘性を考慮しているので，数値計算. Plesset 方程式においても同様である．このことを考慮したう. の結果に顕在化した位相遅れが見られないことから，水の粘. えで圧力波形と界面の振動波形を比べると以下のことが言. 性による位相遅れも無視できると判断できる．変位計測にお. える．初めに，界面の振動波形は t=0 において圧力が最小値. いて，呼吸モードの振動振幅に対応した変位が明確に観測さ. をとる位置で最大値をとるはずであるがそうはならず，変位. れた意義は大きいが，このような圧力と変位との位相関係が. は 0 よりわずかに小さな値をとった．すなわち，バブルは離. 得られたことに関してはその原因は解明できていない．他方，. 脱直後にわずかに収縮した後に膨脹に転ずる．また，本来一. バブル離脱の現象においては，圧力の最低ピーク時刻は離脱. φ=4mm. φ=6mm. (a) Overall and local images of bubble. 125. 125. Pa. 250. ｐ Pa. 250. 0. ｐ. 0. -125. -125. -250 -0.001. 0. 0.001. 0.002 0.003 t sec. 0.004. 0.005. -250 -0.001. 0.006. 0. 0.001. 0.002. 0.003. 0.004. 0.005. 0.006. sec. t. 1. 1. 0.8. 0.8. Intensity (-). Intensity (-). (b) Pressure waveform. 0.6. 0.4. 0.6. 0.4. 0.2. 0.2. 0. -40. -20. 0. 20. x (pixel). 0. 40. -20. 0. 20. 40. x (pixel). 60. (c) Intensity distribution across the interface (t=0) 250. p. 0.15. δ. 125. 250. δ. 125. 0.05. 0. 0. δ pixel. 0.05. -0.05 -125. -0.1. 0. 0.001. 0.002. 0.003. 0.004. 0.005. 0. 0 -0.05. -125. -0.1. -0.15 -0.2 -0.001. p. 0.1. p Pa. δ pixel. 0.1. 0.15. -0.15 -0.001. -250 0.006. 0. 0.001. 0.002. 0.003. t sec. p Pa. (d) Displacement of bubble-water interface Fig.12 Experimental results of breathing mode of detached bubble.. ─8─. 0.004. 0.005. -250 0.006. p Pa. 0.2.

(9) 周波数フィルタを用いた超解像計測の検証と水中気泡呼吸モードの界面変位測定への応用時刻より遅れるとの報告もある 18)．これは，離脱による呼吸. おいては周囲や残留気泡との共鳴を含めた解析が重要とな. モードの振動が離脱による表面張力の不連続的な変化によ. るであろう．. る単純減衰振動ではないことを暗に示している．このことは既報 4)において行った CFD の結果においても確認すること. 5．．結論. ができ，離脱直後の振動の周期はバブル形状が大きく変形す. Optical flow 解析法と時間周波数フィルタを組み合わせた. るために固有振動数は呼吸モードの基本周波数からごく僅. 画像解析法を用いて，水中ノズルからバブルが離脱する際に. かにずれている．本解析では周波数フィルタリングにより呼. 発生する呼吸モードの振動現象の界面変位計測を行った．本. 吸モードの固有振動数を抽出して解析を行っているために. 報告では，数値実験と加振ディスクの振動変位計測による検. 変位計測結果にこのようなずれが発生したものと考えられ. 証実験も合わせて行った．得られた結果をまとめると以下のとおりである． (1)サブピクセルでの画像計測を対象とした本手法によって. サブピクセルでの変位計測において変位の 1/10 の精度での計測が可能であることが数値実験，および加振ディスク変位計測によって実証された (2)本手法において界面移動計測を行う場合には，計測界面. 幅はある程度の大きさの幅が必要であり，サブピクセルの計測においては 10pixel 程度が良好である．また輝度階調は 12bit 以上とすることが望ましい． (3)本手法をバブル離脱時の呼吸モードの界面変位計測に適. 用して変位を可視化計測することが可能であることを示し. Fig.13 Averaged amplitude in breathing mode.. た．計測結果は数値計算によってシミュレーションされた既. る．バブル離脱時の詳細な変動に関しては，今後，更なる非. 報 4)の結果とほぼ一致し，計測の妥当性が示された．. 定常解析が必要である．次に，界面変位の大きさについて整理した結果について述. 謝辞. べる．界面振動は減衰振動であるので，振幅値を整理するにあたって，前章の加振器ディスクの場合とは異なり，振幅の peak to peak の 1/2 の値  をとって整理した．著者らの一連. 本研究の実験遂行において埼玉大学理工学研究科土田大騎君にご協力いただいた．ここに記して謝意を表する．. の研究においてバブル離脱時の初期変動が重要であるとみなしているために，ピークは Fig.12(d)に矢印で示したように変位の最初の極小値のピークと 2 番目の極大値のピークを選んだ．Fig.13 にその結果を示す．横軸は，バブルの等価半径の値をとって示した．ここで縦軸  は長さの次元で示してあるので注意されたい． は，ノズル径が大きくなるほど大きく，今回の計測で得られた振幅値はノズル内径 4mm で. 参 1) 2) 3). 0.8±0.22μm，6mm で 0.9±0.1μm であった．気泡の離脱の状態. は，ノズルへ供給する空気流量によって影響を受ける．ノズ. 4). ル直径が大きい場合には，離脱過程はやや不安定で，ノズル内の残存気体の影響も大きく受けると言える．空気流量を極めて少なくして準静的に気泡を膨張させる場合には周囲の. 5). 水がノズル内に侵入する影響を受けやすいという影響が出てくる．実験では 1 秒間に 1~2 個程度のバブルが生成され. 6). るように調整している．離脱バブルの周囲条件は綿密に考慮しなければならないが，今回の報告ではこれらについての議論はしない．著者らが行った既報の文献４）で報告している. 7). 数値計算によるシミュレーション結果と比較すると，今回の計測結果は良い一致を示していると言える．既報 4)でも考察. 8). しているように，気泡離脱後の運動は，気泡離脱時にノズル内に残留する気泡の影響も受けるうえにタンクなどの境界条件にも影響を受けると考えられる．従って，今後の研究に. 9). ─9─. 考. 文. 献. Plesset, M.S., Prosperetti, A., :Bubble Dynamics and Cavitation, Ann. Rev. Fluid Mech., Vol.9(1977), pp.145-185. Martin, D., Thierry, H., Jean-Claude, B., Valentin, L. :The Minnaert bubble: a new approach, HAL Id : hal-00145867,pp.1-13. Oku, T., Hirahara, H., Akimoto, T. :Visualization of deformation and sound emission from bubble in water using VOF method. 18th Int. Symp. Flow Vis. (2018), pp.1–9. Oku, T., Hirahara, H., Akimoto,T., Tsuchida,D. :Numerical simulation of breathing mode oscillation on bubble detachment, Fluids, Vol.5, Issue 2 (2020), pp.1-17. 田中正行、奥富正敏:画素数の壁を打ち破る複数画像からの超解像技術, 映像情報メディア学会誌, Vo.62, No.3(2008), pp. 337-342. Qi, T., Michael, N.H. :Algorithms for Subpixel Registration, Computer Vision, Graphics and Image Processing, Vol.35, Issue.2(1985), pp.220-233. Gustavo, K., Akram, A., Dennis, M.H. :Interpolation Artifacts in Sub-Pixel Image Registration, IEEE Transactions on Image Processing ,Vol.18, Issue.2(2009), pp.333-345. Istvan, H., Sergiu, N. :Design of Interpolation Function for SubpixelAccuracy Stereo-Vision System, IEEE Transactions on Image Processing,Vol.21, Issue.2(2012), pp.889-898. Hagiwara, M., Abe,M., Kawamata，M. : Estimation Method of Frame Displacement for Old Films Using Phase-Only Correlation,.

(10) 奥. 10). 11). 12). 13). 14). 15). 16). 17) 18). 隆夫，平原. 裕行，姜東赫. Journal of Signal Processing, Vol. 8, No. 5(2004), pp.421-429. Manuel, G., Samuel T., James R.: Efficient subpixel image registration algorithms, Optics Letters, Vol. 33, Issue.2(2008), pp.156-158. Amr, Y., Jiang, L., Mohammad, K. :High-speed Image Registration Algorithm with Subpixel Accuracy, IEEE Signal Processing Letters, Vol.22, Issue.10(2015), pp.1796-1800. Hao-Yu, W., Michael, R., Eugene, S., John, G., Fredo, D., William, F. :Eulerian Video Magnification for Revealing Subtle Changes in the World, ACM Transactions on Graphics, Vol.31, No.4(2012), pp.1-8. Liu. L., Lu, L., Luo, J., Zhang, J., Chen, X., :Enhanced Eulerian Video Magnification, IEEE 7th International Congress on Image and Signal Processing(2014),PP.50-54. Neal, W., Michael, R., Frédo, D., William, T. F., :Phase-Based Video Motion Processing, ACM Transactions on Graphics, Vol.32, No.4(2013), pp.1-9. Neal, W., Hao-Yu, W., Abe, D., Michael, R., Eugene, S., Gautham, M., Justin, C., Oral, B., John, V., William F., Frédo, D.: Eulerian Video Magnification and Analysis, Communications of the ACM, Vol.60, No.1(2017), pp.87-95. Bruce, D.L., Takeo, K. :An Iterative Image Registration Technique with an Application to Stereo Vision, In Proceeding of IJCAI, Vol.2(1981), pp.674–679. Horn, B.K.P., Schumck, B.G.: Determining Optical Flow, Artificial Intelligence, Vol.17 (1981), pp.185-203. Manasseh, R., Riboux, G., Bui, A., Risso, F. :Sound Emission on Bubble Coalescence: Imaging, Acoustic and Numerical Experiments, Proceedings of the 16th Australasian Fluid Mechanics Conference, 16AFMC(2001), pp.857-860.. ─ 10 ─.

(11)