- 機械工学・建築工学・土木工学・航空宇宙工学・原子力工学・電気工学 -
筑波大学大学院 システム情報工学研究科 構造エネルギー工学専攻:機械・土木・建築・航空・宇宙・原子力・電気・環境・エネルギー:修士課程・博士課程・追加募集
EME_logo:機械工学科・建築学科・土木工学科・航空宇宙工学科・原子力工学科・電気工学科の融合体
フリーソフトウェア
PIV

PIV (Particle Image Velocimetry, 粒子画像流速測定法)

開発担当者:榊原 潤 准教授

粒子画像流速測定法(PIV、particle image velocimetry)は、流体に追従する粒子を流動場に混入させ、時間的連続撮影された可視化画像から微少時間dtにおける粒子の変位ベクトルdxを求め、速度ベクトルdx/dtを推定する方法である。

粒子の変位ベクトルを求める方法の一つである画像相関法は、光シートで照明された2次元的な粒子群をビデオカメラ等で撮影し、1時刻目t=t0および2時刻目t=t0+dtにおける粒子画像を得た上で(図1a)、1時刻目の画像における微小な領域(検査領域)内の輝度値分布と2時刻目の画像における領域(探査領域)内の輝度値分布との相互相関関数を求め、その最大値となる位置を検査領域内の粒子群の平均的な移動位置として推定し、変位ベクトルdxを求める方法である(図1b)。
  相互相関関数は粒子画像と同様に空間的に離散化されているため、求められる変位ベクトルは±0.5画素の誤差を伴う。そこで、離散化された相関関数に二次元正規分布を内挿して連続関数とした上で変位ベクトルを求めることで、誤差を±0.1画素程度に減少させる手法(サブピクセル解析)がとられる。ただし、粒子像の大きさが2画素を下回るときには真の変位量と求められる変位量の関係が線形にならない(ピークロッキング)ので、粒子像の大きさに注意する必要がある。

1C-PIV

図1a 1時刻目t=0および2時刻目t=dtにおける粒子画像

1C-PIV

図1b 画像相関法のために撮影された粒子画像、相関関数および速度ベクトル分布

一方、カメラ2台を用いてトレーサ粒子をステレオ視することで、二次元面内における速度三成分を計測するステレオPIVは、従来のPIVシステムを発展させることで比較的容易にシステムを構築できるため、次世代PIVとして注目されている。 ステレオPIVではトレーサ粒子を懸濁した流動場にパルスシート光を照射し、トレーサ粒子からの散乱光を2台のCCDカメラで撮影する。パルスシート光は短い時間間隔dtで2回照射し、それぞれの散乱光をCCDカメラの時間的に連続する2フレームに撮影する。図2は可視化された粒子画像の例である。

 

図2 粒子画像の例;時刻t=0(左)、時刻t=8ms(右)

 図3はシート光面を横切る断面で移動する粒子とCCD素子面に結像される粒子像の幾何学的関係を示している。ここで、シート面中央部を原点とする物理座標P、CCD(L)面およびCCD(R)面における画像座標XLXRをそれぞれ定義する。物理座標の位置Pにある粒子がdt後にP'へ移動したとき、CCD(L)面に結像される粒子像は第1フレームの画像座標XLから第2フレームのXL+dXLへ、CCD(R)面に結像される粒子像は画像座標XRからXR+dXRへそれぞれ移動する。ここで、点XL+dXLおよびXR+dXRには物理座標における直線AC およびBD上の全ての点がそれぞれ投影される。よって、画像座標における任意の点とそれに対応する物理座標における直線の関係が後述する位置校正により既知であれば、点XL+dXLおよびXR+dXRに対応する直線ACおよびBDがそれぞれ求まり、さらに、この2直線の交点を求めることで点P'が定まる。ここで、2直線は種々の誤差要因により完全に交わることはあり得ないので、その代わりに2直線の最短距離を結ぶ直線、すなわち2直線に互いに垂直に交わる直線を求め、その2交点の中間点をP'と定義する。最後に、点Pにおける速度ベクトルV

(1)
として近似する。  CCD面上の粒子移動ベクトルdXL, dXRは2CPIVに基づいて各々の画像から算出する。例えば、CCD(L)により得られた第1フレーム画像のXLにおける粒子パタンと第2フレームにおける近傍の粒子パタンとの相互相関関数を求め、相関係数の最大となるXL+dXLを探索する。  以上示した原理の他に、物理座標から画像座標への写像関数の微係数テンソルを求めておく方法も提案されている 。

図3 パルスシート光内の粒子位置とCCD素子面に結像される粒子像との関係



画像座標の任意の点とそれに対応する物理座標における直線の関係は、AB面およびCD面に図3に示すような校正板を挿入し、その画像から画像座標と物理座標を対応づける写像関数を作成する。校正板は既知の間隔の格子線を描いたものであり 、アームを介して微動送り台に固定されている。校正板には基準となる点(基準点)を示すためのマークを書き込んでおき(図4参照)、このマーク近傍(例えば右下)の格子点の位置を定めておく。


校正板

図4 校正板画像


画像座標の任意の点とそれに対応する物理座標における直線の関係は、AB面およびCD面に図3に示すような校正板を挿入し、その画像から画像座標と物理座標を対応づける写像関数を作成する。校正板は既知の間隔の格子線を描いたものであり 、アームを介して微動送り台に固定されている。校正板には基準となる点(基準点)を示すためのマークを書き込んでおき(図4参照)、このマーク近傍(例えば右下)の格子点の位置を定めておく。

校正手順は次のとおりである。
(1)校正板をシート光面中央(AB面, z=z0)に位置するよう設置する。
(2)校正板をCCDカメラで撮影し画像を記憶する(校正画像)。ここで、基準点が画像中にあるようにする。
(3)校正画像の格子点の画像座標を求めるために、あらかじめソフトウェア的に用意しておいた格子パタンと校正画像全体の相互相関関数を求め、相関係数のピーク位置を格子点とする。この段階で、各格子点の画像座標のリスト(格子点リスト)が出来上がる。

(4)次に、各格子点の物理座標を求める。物理座標が既知である基準点を対象点とし、対象点の4近傍にある格子点を格子点リストから探しだし、それら4点の物理座標を対象点の物理座標と格子間隔から算出する(図4(a))。これにより、対象点に加えてその近傍4点の物理座標が既知となるので、今度はこれら4点をそれぞれ対象点として同様の操作を再帰的に繰り返せば、全ての格子点座標を算出することができる(図4(b)(c))。こうして、各格子点の画像座標と物理座標の対応リストを作る。

(5)画像座標と物理座標の対応リストから最小二乗法を用いて画像座標から物理座標への写像関数fx(X,Y)fy(X,Y)およびその逆関数FX(x,y)FY(x,y)のそれぞれの係数ai, bi, ci, diを求める。ここでx,yは物理座標、X,Yは画像座標であり、写像関数は3次多項式とした。

(2a) (2b) (3a) (3b)

以上により、z=z0における写像関数が求まる。
(6)つぎに、微動送り台を用いて校正板をその面と垂直な方向に移動する(CD面, z=z1)。移動量は画像中で格子点が10〜20画素移動するのに相当する程度であり、後述の例では2mmである。
(7)上記(2)〜(5)と同様にしてz=z1における写像関数を求める。 以上の操作により、画像における任意の点に対応するz=z0およびz=z1面における物理的な位置が得られることになる。

piv.zip(1.3MB):PIV、ステレオPIV、校正プログラム一式(実行ファイルのみ) 

利用規約に同意の上、ダウンロードしてください。

同意の上ダウンロード

Copyright 2010 All Right Reserved. University of Tsukuba, Graduate School of Systems and Information Engineering
Tennodai 1-1-1, Tsukuba, Ibaraki, 305-8573, Japan