JP2007033306A

JP2007033306A - 流体の流動計測システム及びその計測方法

Info

Publication number: JP2007033306A
Application number: JP2005218768A
Authority: JP
Inventors: Harutsugu Mori; 治嗣森; Hideaki Tezuka; 英昭手塚; Koichi Hishida; 公一菱田; Itsusei Tachibana; 逸成橘
Original assignee: Tokyo Electric Power Co Inc; Keio University
Current assignee: Keio University; Tokyo Electric Power Co Holdings Inc
Priority date: 2005-07-28
Filing date: 2005-07-28
Publication date: 2007-02-08

Abstract

【課題】粒子軌跡追跡法を用いた流動計測、特に３次元流動計測において、より信頼性が高くより精度が高く精密な測定が可能な流体の流動計測システム及びその計測方法。
【解決手段】レーザ発振装置１１から発振されたレーザ光を流体流動場１４内にシート状に投入させるレーザシート形成用光学系１３と、レーザシート１７で照明されて形成される２次元粒子画像と２次元粒子軌跡画像を選択的に撮像する画像撮像手段２９Ｌ、２９Ｒと、微小時間異なる２時刻における２次元粒子画像の輝度パターンを比較・解析して粒子群の移動方向及び移動量を計測する第１画像処理手段と、微小時間異なる２時刻における２次元粒子軌跡画像の輝度パターンを比較・解析して個々の粒子の移動方向及び移動量を計測する第２画像処理手段とを備え、第２画像処理手段で、第１画像処理手段で計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定める。
【選択図】図１

Description

本発明は、流体の流動計測システム及びその計測方法に関し、特に、閉空間内を流動する流体の流速及び流れ方向を高精度かつ精密に測定する流体の流動計測システム及びその計測方法に関するものである。

近年、圧力容器内のような外界と環境を異にする閉空間の流動場を把握することはプラント等の効率や安全性の向上等から非常に要望されている。しかし、参照窓の必要性や光学系等の設置の困難さから粒子画像流速測定法（以下、ＰＩＶ（Particle Imaging Velocimetry）と呼ぶ。）による測定は従来では不可能であった。非特許文献１では、一対のファイバスコープを使用することで、閉空間内等へ適用可能であるステレオＰＩＶシステムを構築し、２次元３方向成分の速度測定可能なシステムを構築した。しかし、この場合、ファイバスコープを通すことで画像が劣化し粒子の情報量が極めて減少することにより、測定精度の向上が必要とされていた。

このような状況下、特許文献１において、ファイバスコープを使用したＰＩＶシステムにおいて、ファイバスコープを通して撮像された粒子画像を各ピクセル毎に輝度パターン（輝度値）の加算平均値を求め、各々の粒子画像の輝度パターン（輝度値）から加算平均値を引くことで、光ファイバのファイバ配列の輝度パターン除去を行ってファイバスコープ像のノイズ除去を行い、流体の流速や流れ方向等の流動状態を精度良く正確に測定する方法が提案されている。

さらに、特許文献２、非特許文献２において、ＰＩＶを発展させた粒子追跡法（以下、ＰＴＶ（Particle Tracking Velocimetry ）と呼ぶ。）においてファイバスコープを使用し、かつ、粒子の流跡を記録する粒子軌跡追跡法（以下、ＰＴＶ−ＳＳ（PTV with Streak Searching ）と呼ぶ。）が提案されている。そのとき、特許文献１と同様のファイバスコープ像のノイズ除去を行っている。

さらに、非特許文献３において、ファイバスコープを使用したＰＴＶ−ＳＳを３次元に拡張したファイバスコープステレオＰＴＶ−ＳＳを提案している。

さらに、非特許文献４において、ＰＴＶにおいて、各粒子の流跡を検出するのに、各粒子の流跡の検出の閾値をその背景に応じて個々に定める動的２値化法を提案している。
特開２００２−２２７５９号公報特開２００３−８４００５号公報日本機械学会熱工学講演会論文集，Ｎｏ．０２−２２（２００２），ｐｐ．７〜８ＥｘｐｅｒｉｍｅｎｔｓｉｎＦｌｕｉｄｓ，３３（２００２），ｐｐ．７５２〜７５８日本機械学会熱コンファレンス２００３講演論文集［２００３−１１．１５−１６，金沢］ｐｐ．２７３〜２７４Ｍｅａｓ．Ｓｃｉ．Ｔｅｃｈｎｏｌ．，Ｖｏｌ．１１（２０００），ｐｐ．６０３〜６１６

本発明は従来技術のこのような状況に鑑みてなされたものであり、その目的は、粒子軌跡追跡法を用いた流動計測、特に３次元流動計測において、より信頼性が高くより精度が高く精密な測定が可能な流体の流動計測システム及びその計測方法を提供することである。

上記目的を達成する本発明の流体の流動計測システムは、レーザ光を発振させるレーザ発振装置と、発振されたレーザ光を流体の流動場内にシート状に投入させるレーザシート形成用光学系と、このレーザシート形成用光学系からのレーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を選択的に撮像する画像撮像手段と、前記画像撮像手段で撮像された微小時間異なる２時刻における２次元粒子画像の輝度パターンを比較・解析して粒子群の移動方向及び移動量を計測する第１画像処理手段と、前記画像撮像手段で撮像された微小時間異なる２時刻における２次元粒子軌跡画像の輝度パターンを比較・解析して個々の粒子の移動方向及び移動量を計測する第２画像処理手段とを備え、前記第２画像処理手段では、前記第１画像処理手段で計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めるように構成されていることを特徴とするものである。

この場合、前記レーザ発振装置は半導体レーザで構成され、前記半導体レーザの発光時間の差によって２次元粒子画像と２次元粒子軌跡画像を形成することが望ましい。

また、前記レーザシート形成用光学系からのレーザシートの両側から同時に２次元粒子画像と２次元粒子軌跡画像を選択的に撮像する２つの画像撮像手段を備え、各々の前記画像撮像手段に対してそれぞれ前記第１画像処理手段と、前記第２画像処理手段とを備え、各々の前記第２画像処理手段では、各々の前記第１画像処理手段で計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めるように構成することができる。

その場合、前記第２画像処理手段各々で計測された個々の粒子の移動方向及び移動量を相互に対応付けて個々の粒子の３次元的な移動方向及び移動量に再構成する３次元再構成手段を有し、個々の粒子の３次元的な移動方向及び移動量を計測することができる。

また、前記画像撮像手段は、前記レーザシート形成用光学系からのレーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を光伝送する画像伝送手段と、光伝送された２次元粒子画像と２次元粒子軌跡画像を撮像する撮像手段とを備えたものとして構成することができる。

その場合に、前記画像伝送手段は、多数本の光ファイバを束ねて一体化させ、両端面が平面加工された可撓性のイメージガイドで構成され、前記撮像手段は、ＣＣＤカメラで構成されており、前記イメージガイドは、対物レンズにより一方のファイバ端面に結像された画像を光ファイバで各画素に分解してＣＣＤカメラ側の他端面まで同一画像として伝送するように構成することができる。

また、以上において、前記レーザシート形成用光学系と前記画像撮像手段とは、一体的に組み合わせて構成されていることが望ましい。

また、前記第１画像処理手段、前記第２画像処理手段はそれぞれ２次元粒子画像と２次元粒子軌跡画像から前記イメージガイドに基づくファイバノイズを除去するファイバノイズ除去手段を備えていることが望ましい。

また、前記第２画像処理手段で計測された個々の粒子の移動方向及び移動量をサブピクセル精度化するためのサブピクセル処理手段を備えていることが望ましい。

本発明の流体の流動計測方法は、レーザ発振装置から発振されるレーザ光を流体の流動場にシート状に照射してレーザシートを形成し、レーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を選択的に画像撮像手段で撮像し、撮像された微小時間異なる２時刻における２次元粒子画像の輝度パターンを比較・解析して粒子群の移動方向及び移動量を計測し、また、撮像された微小時間異なる２時刻における２次元粒子軌跡画像の輝度パターンを比較・解析して個々の粒子の移動方向及び移動量を計測し、その個々の粒子の移動方向及び移動量の計測の際に、先に計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めることにより、流体流動場における流体の流動を計測することを特徴とする方法である。

本発明による流体の流動計測システム及びその計測方法によれば、流体の流動場にシート状に照射してレーザシートを形成し、レーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を選択的に撮像し、撮像された微小時間異なる２時刻における２次元粒子画像の輝度パターンを比較・解析して粒子群の移動方向及び移動量を計測し、また、撮像された微小時間異なる２時刻における２次元粒子軌跡画像の輝度パターンを比較・解析して個々の粒子の移動方向及び移動量を計測し、その個々の粒子の移動方向及び移動量の計測の際に、先に計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めることにより、複雑な流動場における個々の粒子の流速や流れ方向を測定して流体の流速や流れ方向等の流動状態を、より信頼性高く、より精度高く、より精密に測定することができる。

また、本発明によれば、複雑な流動場における流体の流動状態の測定を画像処理手段で粒子画像と粒子軌跡画像を処理することにより行うことができるので、遠隔地から流体の流動状態を測定でき、環境的に厳しい閉空間内の流体の流動状態を高信頼性、高精度、高精密に測定することができる。

以下に、本発明の流体の流動計測システム及びその計測方法を図面を参照して説明する。

図１は、本発明の流体の流動計測システムの１実施例の概要を示す図である。

本発明の流体の流動計測システム１０には、高出力のレーザ光を発振させるレーザ発振装置１１が光源として備えられる。レーザ発振装置１１には、例えば、小型で、発光のタイミングを制御することが容易で、高輝度、高繰り返し周波数の可能な半導体レーザを用いる。

レーザ発振装置１１から発振されたレーザ光は、送光用光ファイバ１２を通ってレーザシート形成用光学系１３に案内される。レーザシート形成用光学系１３はレーザ発振装置１１から発振されたレーザ光を流体流動場１４にシート状に投射させるようになっている。

流体流動場１４は、外界と環境を異にする閉空間内に形成される。また、場合に応じて、外界に開放された形態で流体流動場１４が形成されていてもよい。流動場１４は、気体、液体の内の少なくとも何れかの流体からなる。具体的には、例えば、原子炉圧力容器のダウンカマ部、炉心シュラウド、火力発電プラントの熱交換器や蒸気発生器内等の閉空間における流体の速度場を対象とする。また、船舶、航空機あるいは自動車等のような、水流や気流と接触する構造部位や機材周りにおける流体の流動場も対象とすることができる。

レーザシート形成用光学系１３は、流体流動場１４にレーザ発振装置１１からの発振レーザ光をシート状に照射して流体流動場１４にレーザシート１７を形成している。このレーザシート１７により流体流動場１４中に混入された微細なトレーサー粒子を照明して可視化している。

レーザシート形成用光学系１３は、図２（ａ）にレーザシート１７のシート面に沿う断面を、図２（ｂ）にそれに直交する断面を示すように、レーザシート１７の作成のためのレーザ光照射光学系である。レーザ発振装置１１として例えば上記のような半導体レーザを用いると、半導体レーザからの発振レーザ光は、送光用光ファイバ１２から大きな広がり角を持って出射される。出射されたレーザビームは軸対称の光となっておらず、非点較差があり、レーザビームの形状変換は複雑となる。

このため、レーザシート形成用光学系１３には、送光用光ファイバ１２から出射されたレーザビームからレーザシート１７を形成するためのレーザシート形成用光学手段としてのレンズ群１８が用いられる。このレンズ群１８は、送光用光ファイバ１２から出射したレーザビームを例えば３枚の平凸レンズ１９ａ、１９ｂ、１９ｃによってコリメートを行い、平凹シリンドリカルレンズ１９ｄによってレーザシート１７を作成する構成となっている。このレーザシート形成用光学系１３によって出力されるレーザーシート１７の厚さは例えば３．０ｍｍとなる。

レーザシート形成用光学系１３から流体流動場１４にレーザシート１７が照射され、流体の流動場１４がシート状に可視化される。そして、レーザシート１７の測定範囲１７ａに両側から対向して左右の画像伝送手段２２Ｌ、２２Ｒのそれぞれの対物レンズ２７ａが配置され、各対物レンズ２７ａを経て結像された左右の粒子画像は、左右の画像伝送手段２２Ｌ、２２Ｒを介して左右の画像撮像手段２９Ｌ、２９Ｒに送られ、そこで撮像されるようになっている。この実施例の流動計測システムにおいては、以下のように、この画像伝送手段２２Ｌ、２２Ｒとして、イメージガイド２６と対物レンズ２７ａから構成されるファイバスコープを使用している。

ここで、左右の画像伝送手段２２Ｌ、２２Ｒの先端のプローブ部は、図３（ａ）に斜視図、（ｂ）に断面図を示すように、レーザシート形成用光学系１３を間に挟んで対称に各画像伝送手段２２Ｌ、２２Ｒのイメージガイド２６の先端部が一体に配置され、各イメージガイド２６の先端に、光路を９０°曲げてレーザシート１７の測定範囲１７ａの両側に向くようにする４５°直角プリズム２７ｂと、その入射側に配置され測定範囲１７ａの粒子画像を撮像するための対物レンズ２７ａとが一体に配置されている。

そして、この左右の画像伝送手段２２Ｌ、２２Ｒで伝送されたレーザシート１７上の測定範囲１７ａの左右の粒子画像は、左右の画像撮像手段２９Ｌ、２９Ｒで撮像可能となっている。

左右の画像撮像手段２９Ｌ、２９Ｒは、それぞれ左右の画像伝送手段２２Ｌ、２２Ｒで伝送された粒子画像を撮像する撮像手段としてのカメラレンズ２８を備えたＣＣＤカメラ２４Ｌ、２４Ｒを有する。画像伝送手段２２（左右を区別するＬ、Ｒは省く。）は、図４に示すように、数千本から数万本の光ファイバ２５を束ねて一体化（溶融）させた可撓性のイメージガイド２６が用いられる。イメージガイド２６は、図４に示すように、束ねられた光ファイバ２５の両端面の各位置が正確に対応するように並設されて束ねられる。イメージガイド２６の両端面は平面に仕上げられ、対物レンズ２７ａにより一方の光ファイバ端面に結像された画像を各画素光ファイバとしての光ファイバ２５に分解し、同一画像をＣＣＤカメラ２４側の光ファイバ端面まで伝送している（図４では、４５°直角プリズム２７ｂは省いてある。）。光ファイバ端面からカメラレンズ２８を経てＣＣＤカメラ２４で撮像される。なお、画像撮像手段２９Ｌ、２９Ｒは、以上のようなイメージガイドを備えずに、レーザシート１７上の粒子画像を直接撮像できるようにしてもよい。

イメージガイド２４の画像の限界（解像度）は、単一の光ファイバ２５の直径と配列の仕方によって決まり、像の見え方は、イメージガイド２６に束ねられる光ファイバ２５の総本数と並べ方によって決まる。光ファイバ２５のコア径より小さな光量の光分布は、光ファイバ２５を伝わる間に平均化され、各コア断面内では明るさは一様となる。したがって、イメージガイド２６は光ファイバ２５のコア径が画素（ピクセル）としての空間分解能の最小単位となる。

実際のイメージガイド２６は、多数本の光ファイバ２５が俵積み（六方稠密配列）構造で整然とかつ稠密に配列されている。このイメージガイド２６で光ファイバ間隔と略等しい間隔の格子状物体を観察すると、物体の像の強度分布と無関係なモアレ干渉縞が生じる。

また、光ファイバ２５間にクラッド部分が融着されているイメージガイド２６では、クラッドの厚さが薄くなると、１つの光ファイバ２５のコアに入った光が隣のコアに漏れてしまう場合が生じる。このため、像の境界が滲み、像のコントラストが低下する。このため、イメージガイド２６は、六方稠密配列構造を採用したとき、クラッドの厚さが薄過ぎないことが条件となる。

一方、本発明の流体の流動計測システム１０には、シグナルジェネレータ３０を備えており、パソコン３１内のソフトウェアによりシグナルジェネレータ３０を制御することで、ＴＴＬ信号を任意に作成することができ、その信号をレーザ発振装置１１のレーザドライバ３２、ＣＣＤカメラ２４Ｌ、２４Ｒ両方に伝送し、ＣＣＤカメラ２４Ｌ、２４Ｒのフレームの切り替わり、レーザ発振装置１１のレーザの照射時間・照射間隔を任意に調節可能になっている。なお、レーザドライバ３２とレーザ発振装置１１の間には、レーザ発振装置１１の半導体レーザの冷却装置としてのペルチエドライバ３３が配置されており、レーザドライバ３２からこのペルチエドライバ３３を介してレーザ発振装置１１が駆動されるようになっている。

また、ＣＣＤカメラ２４Ｌ、２４Ｒはそれぞれ左右のインターフェース３４Ｌ、３４Ｒを介して撮像された左右の粒子画像をパソコン３１へ転送するように構成されている。そして、パソコン３１に取り込まれた粒子画像は即座に以下に説明する画像処理ソフトウェアによって、粒子画像流速測定法（ＰＩＶ）と粒子軌跡追跡法（ＰＴＶ−ＳＳ）とを利用して３次元速度ベクトルデータに変換することが可能となっている。

以上のような流動計測システムにおいて、図５に示すような画像処理ソフトウェアによって高信頼性、高精度、高精密に流体流動場１４の３次元速度ベクトル情報を取得する。以下に、図５の処理の流れを順に説明する。

同時刻に左右のＣＣＤカメラ２４Ｌ、２４Ｒで撮像され、インターフェース３４Ｌ、３４Ｒを通じてパソコン３１に送られた左右のカメラ画像４１Ｌ、４１Ｒそれぞれに対して、以下のような処理を並列して行う。図５では、ＰＩＶ処理４２Ｌ、ＰＴＶ−ＳＳ処理４３Ｌ、光ファイババンドルによる画像劣化のサブピクセル処理４４Ｌまでは、左カメラ画像４１Ｌについてのみ図示され、右カメラ画像４１Ｒについては図示を省いているが、右カメラ画像４１Ｒについても同様に行われる。

まず、ＰＩＶ処理４２Ｌは、特許文献１で提案された方法と同様に、ファイバスコープのイメージガイド２６の光ファイバ２５の配列パターンを除去した後に、通常のＰＩＶ処理を行う処理である。この処理で算出されるベクトルは、ＰＴＶ−ＳＳ処理４３Ｌにおける処理の際の補助的なベクトル情報となるものである。

ＰＴＶ−ＳＳ処理４３Ｌは、ＰＩＶ処理４２Ｌの後の処理で、従来のＰＴＶ−ＳＳ法に加えて前処理として行ったＰＩＶ処理４２Ｌで算出された流速ベクトル（ＰＩＶベクトル）の情報を利用してさらに精度良く粒子の対応付けを可能とするものである。

まず、ＰＩＶ処理４２Ｌについて説明する。ＰＩＶ処理は、得られた粒子画像を輝度の分布であると考え、粒子群の輝度パターンが微小時間内に移動する量を定量的に解析するものである。そのために、粒子移動量の解析には相互相関法が多く用いられる（相互相関ＰＩＶ処理４２２）。しかし、ＣＣＤカメラ２４Ｌで取得されるカメラ画像４１Ｌには、ファイバスコープのイメージガイド２６を通して撮像しているため、粒子群の輝度パターンだけでなく、光ファイバ２５の配列の輝度パターンも存在する。このため、光ファイバの配列が粒子群の輝度パターンの相関値のピークの検出の際に弊害となり、結果として粒子群の誤対応（過誤ベクトル）を招く。その誤対応の例を図６に示す。図６を見ると、明らかに粒子群ではなく、光ファイバの配列パターンを相関値のピークとして検出していることが分かる。

このように、カメラ画像４１Ｌ内に記録された光ファイバ２５の配列が粒子群の輝度パターンに対してノイズとなって過誤ベクトルの発生の要因となっている。そこで、ＰＩＶ処理４２Ｌにおいては、ファイバノイズの除去処理４２１を相互相関ＰＩＶ処理４２２の前に行う。このファイバノイズの除去処理４２１は次のようにして行う。

ファイバスコープのイメージガイド２６を通してＣＣＤカメラ２４Ｌで撮像されたカメラ画像４１Ｌを処理し、時系列に得られた粒子画像を各ピクセル毎に輝度パターン（輝度値）の加算平均値を求め、各々の粒子画像の輝度パターン（輝度値）から加算平均値を引くことで、光ファイバ２５の配列の輝度パターン除去を行う。その手法を説明すると、測定データ（粒子画像（カメラ画像）４１Ｌにおける輝度値）をＩとすると、
Ｉ＝ｓ＋ｎ・・・（１）
で表わされる。ただし、ｓは信号成分（光ファイバ配列の輝度パターン）、ｎは雑音成分（粒子の輝度パターン）を意味する。

すなわち、信号（光ファイバ配列の輝度パターン）に雑音（粒子の輝度パターン）が重なったものを測定しているものと考える。加算平均はＳＮ比が小さく、信号と雑音の構成周波数にも大差のない場合、また、同じ条件で何度も測定を繰り返すことが可能な場合に有効な手法である。式（１）を時系列として表すと、次のように表せる。

Ｉⁱ（ｋ）＝ｓⁱ（ｋ）＋ｎⁱ（ｋ）・・・（２）
ここで、ｉはｉ枚目の画像、ｋはｉ枚目の画像のｋ番目のピクセル輝度値を表す。

粒子画像の撮影において、取り込んだ粒子画像の枚数をＭ枚とすると、Ｍ枚の粒子画像に対する加算平均Ｉ（ｋ）は、次のように表せる。

_M
Ｉ（ｋ）＝（１／Ｍ）Σ Ｉⁱ（ｋ）
ⁱ⁼¹
_M _M
＝（１／Ｍ）Σ ｓⁱ（ｋ）＋Σ ｎⁱ（ｋ）・・・（３）
ⁱ⁼¹ ⁱ⁼¹
光ファイバ配列の輝度パターンｓは、同一の条件（送光系・受光系）であれば、連続する粒子画像においては同じパターンが出現するため、Ｍ回の加算によりＭ倍の輝度値となり、これを平均する（Ｍで除する）と、元の輝度値のままである。すなわち、加算平均操作により、粒子の存在しない光ファイバ配列のみ撮った画像を得ることができ、光ファイバ配列の存在する粒子画像４１Ｌから各画素毎に光ファイバ配列の輝度値（Ｉ（ｋ）の各画素毎の輝度値）を引けば、光ファイバ配列の輝度パターンの除去が可能となる。

図６と同じ例において、上記のファイバノイズの除去処理４２１を施した測定データに、次のＰＩＶ処理４２２で相互相関計算を行って得られた流速ベクトル（ＰＩＶベクトル）分布を図７に示す。図６と比較して、ファイバノイズの除去処理４２１を行うことで、過誤ベクトルがよく取り除かれていることがよく分かる。

なお、図８は、ファイバノイズの除去処理４２１によって光ファイバ配列の輝度パターンが除去され、粒子像のみが残る画像処理の流れを任意のラインの輝度値分布で示した図である。

ファイバノイズの除去処理４２１で光ファイバ配列の輝度パターンの除去を行った粒子画像は、次の相互相関ＰＩＶ処理４２２で、相互相関法に基づく粒子画像流速測定処理を行う。ＣＣＤカメラ２４Ｌで撮像された第１時刻における２次元粒子画像（処理４２１でファイバノイズの除去が行われている。）をデジタル処理をした粒子画像を「参照画像」、第１時刻とは微小時間異なる第２時刻における２次元粒子画像（同様に、処理４２１でファイバノイズの除去が行われている。）のデジタル処理粒子画像を「探索画像」、上記参照画像中のある限定された矩形領域を「参照窓画像」、探索画像中の限定された矩形領域を「探索窓画像」と定義する。デジタル化された粒子画像中のある１点の値を「輝度値」と定義し、この輝度値が矩形領域に分布しているものを「輝度パターン」と称する。

「参照窓画像」及び「探索窓画像」は輝度値の矩形状分布を表わしており、輝度パターンを示している。

粒子画像から流体中の粒子群の移動量を推定するアルゴリズムは、得られた粒子画像の輝度分布の変化であると判断し、レーザシート１７上の粒子群の輝度パターンが微少時間内に移動する量を定量的に解析するものである。

粒子群移動量の解析には、相互相関法が用いられ、この処理４２２においても相互相関法の解析手法が用いられる。

相互相関法による粒子群の移動量推定は、「探索窓画像」と「参照窓画像」との間で次式で表わされる相関値Ｒを持つことが知られており、この相関値Ｒをもって輝度パターンの類似度を求め、比較検討して２画像間の粒子群の移動量を求める手法である。

式（４）及び式（５）に示すように、相関値Ｒは、探索窓画像及び参照窓画像の輝度値の二乗平均で除する（正規化）ことによって、−１から１までの値をとる。相関値Ｒが１のときに２画像は完全に一致し、相関値Ｒの値が大きい程窓画像間の類似度が大きいことを意味する。

式（４）及び式（５）において、Ｉ₁及びＩ₂は参照窓画像と探索窓画像の各画素の輝度を表し、ξ、ηは探索窓画像と参照窓画像の相対的な位置で表わす。相関値Ｒが最大となる位置ξ、ηが画像内での粒子群の移動量ΔＸ′、ΔＹ′に相当する。参照窓画像の大きさ及び参照窓画像との相対的な位置は、予測される最小及び最大速度から決定される。参照窓画像の大きさｎ×ｍは参照窓画像の実座標系におけるその大きさが流れ場に存在する最小の渦よりも小さく、かつ、参照窓画像の中にトレーサ粒子を６〜８個以上含むことが望ましい。

以上のように、ファイバノイズの除去処理４２１と相互相関ＰＩＶ処理４２２を経ることで、レーザシート１７の測定範囲１７ａ内の図７に例示されるような流速ベクトル（ＰＩＶベクトル）分布情報が得られる。

そして、ＰＩＶ処理４２Ｌの最後の処理として、誤ベクトル除去及び空間平均処理４２３を行う。この処理では、上記のファイバノイズの除去処理４２１と相互相関ＰＩＶ処理４２２を経て算出されたベクトルに対して誤ベクトルと考えられるものを除去するため、正規分布に基づいたばらつき誤差の棄却を行う。相互相関ＰＩＶ処理によって算出されたベクトル間隔をＮピクセルとすると、算出された全てのベクトル（ｕ，ｖ）に対して、各々のベクトルの周囲２Ｎピクセルにあるベクトル成分が例えば平均値±３σ（σ：標準偏差）以内に収まらない場合は、誤ベクトルとして棄却する。もちろん、この基準値は±３σに限定されず、他の値を用いてもよい。そして、残った位置の異なるベクトルを、空間的に規則正しい交点座標位置近傍のものを抽出して平均化し、交点座標上で代表値にする空間平均処理を行う。

これらの処理を行った結果のベクトルをＰＩＶ処理によるベクトルと定義し、これを次のＰＴＶ−ＳＳ処理４３ＬでＰＴＶ−ＳＳ法を行う際の補助的なベクトル情報とする。

なお、以上のＰＩＶ処理４２Ｌによって流体流動場１４のベクトルの概略の情報が取得できる。しかし、粒子群の移動量を見積もるＰＩＶ法では、本来参照窓画像内に最低でも粒子は６〜８個あることが望ましいとされるため、参照窓画像を比較的大きくとらなくてはならなくなり、計算量が莫大に増加しアルゴリズムとしては現実的ではないものとなる（探索する領域が増えると、その４乗に比例して計算量は増加する。）。つまり、取得される画像パラメータ（粒子像径、数密度）が非常に制限を受け、様々な実流動場への対応がＰＩＶ処理のみでは困難となってしまう。そこで、上記のＰＩＶ処理４２Ｌでは粒子数を減少させてでも参照窓サイズを小さく限定し、精度を落としベクトルの概略を取得するというスタンスで解析を行う必要がある。

さて、本発明においては、ＰＩＶ処理４２Ｌの次段階の処理として、計測精度を高めるためのＰＴＶ−ＳＳ処理４３Ｌを適用するが、その際に生じる大きな問題の１つに、粒子の誤対応の問題がある。粒子の誤対応とは、図９に示すように、光ファイババンドルによる画像の低解像度化によって粒子情報量が減少することにより、画像間での粒子の対応付けの間違いが生じることである。図９は、粒子の誤対応を説明するための図であり、この図では、分かりやすいように、参照画像、探索画像を同一画像に表示している。正しい対応付け結果は、図９（ａ）のように、大きな矢印が指す位置であるが、上記の理由による情報量の不足により、図９（ｂ）の矢印のように、粒子の誤対応が起きる可能性が高い。

そこで、本発明においては、ＰＴＶ−ＳＳ処理４３Ｌの左右のカメラ画像４１Ｌ、４１Ｒとして、レーザシート１７により流体流動場１４を照明するレーザ発振装置１１に用いている半導体レーザの発光時間を比較的長くとって粒子の流跡を記録した画像を用い、粒子群ではなく個々に粒子の移動量を見積もる粒子軌跡追跡法（Particle streak tracking velocimetry：ＰＴＶ−ＳＳ）（特許文献２、非特許文献２）を適用することで、ファイバスコープのイメージガイド２６を通した画像からベクトル情報を算出している。

ＰＴＶ−ＳＳ処理４３Ｌは、粒子検出の際のファイバノイズの除去処理４３１と、その後の粒子検出処理４３２、ＰＩＶベクトルを用いたＰＴＶ−ＳＳ法４３３の３段階に分かれている。これらの処理について以下に説明する。

ファイバノイズの除去処理４３１について説明する。ＰＩＶ処理４２Ｌにおいては、ファイバノイズの除去処理４２１で、光ファイバ配列の存在する粒子画像４１Ｌから各画素毎に光ファイバ配列の輝度値である粒子画像の加算平均画像Ｉ（ｋ）を減算する手法を用いてノイズ除去を行った。このファイバノイズの除去処理４３１においても、同様の手法を採用してもよい。

しかしながら、この処理では加算平均を用いる都合上、どうしても必要成分である粒子の輝度変化の情報が多少なりとも損失してしまう。このために、粒子検出処理４３２により粒子検出数が減少してしまうことが考えられる。また、ＰＩＶ処理４２Ｌにおいては、ファイバノイズにより起因する最大の問題は、粒子検出においては光ファイバの配列パターンを相関値のピークとして検出することであり、ＰＩＶ処理における粒子検出の際は、この問題は発生しない。

そこで、本実施例においては、ＰＴＶ処理における粒子画像４１Ｌのファイバノイズの除去処理４３１として、輝度成分情報をできる限り落とさないようにするために、粒子画像４１Ｌに対して光ファイバ２５のパターンである画像の高周波数成分を取り除く高周波数フィルタとして空間平均処理を行う。この処理を行うに当たって、イメージガイド２６の空間周波数を求め、光ファイバ２５のコア径がＣＣＤカメラ２４ＬのＣＣＤアレイの何ピクセルに相当するかを算出する。図１０（ａ）、（ｂ）は実際に撮影したイメージガイド２６の画像の任意ラインの輝度値分布と周波数解析の例である。光ファイバ２５のコアとクラッドの明暗による輝度値の変化が支配的な空間周波数として検出されると考えられ、この例では、０．３２０pixels^-1がＣＣＤカメラ２４Ｌにおいて撮影されるイメージガイド２６の支配的な空間周波数となる。そこで、この例では光ファイバ２５のパターンの成分を取り除くために、３．１２４pixels以下の範囲において空間平均処理を行うことにより、０．３２０pixels^-1以上の空間周波数を取り除いて、ファイバノイズの除去処理を行っている。このファイバノイズの除去処理のためには、上記のような空間平均処理の代わりに、電気信号処理による高周波数の空間周波数を除去する空間周波数フィルタリング処理によってもよい。

以上のようにして光ファイバ配列の存在する粒子画像４１Ｌからファイバノイズ除去処理を行った画像から次の粒子検出処理４３２で粒子検出処理を行う。ＰＴＶ処理を用いてベクトル情報を算出する際に重要な問題の１つとして、粒子の特定を正確に行う点がある。撮像された粒子画像４１Ｌは輝度の濃淡によって表されるので、粒子の特定の方法としては、ある輝度の閾値を用いた単純２値化法や、マスク演算子を用いた粒子抽出法がある。しかしながら、単純２値化法は、画像のノイズ等の要因により粒子間の輝度変動が大きい場合には粒子の検出数が大幅に減少したり、一部分重なっている粒子を同一粒子と認識するといった問題点がある。また、マスク演算子を用いた粒子抽出法には粒子形状が推定パターンと一致しなければなないという問題点がある。

そこで、本発明においては、２値化法の閾値をそれぞれの粒子毎に決定することで、輝度変動が大きい場合への適用が可能な動的２値化法（非特許文献４）を適用する。動的２値化法では、個々の粒子に対して異なる閾値を設定することにより、画像全体における輝度分布のばらつきによって粒子検出精度が低下するという問題を解消している。この動的２値化法の概要を図１１に示す。初期閾値は、画像を分割し（図１２）、その領域内における輝度の平均値Ｉ_avgを利用する。単に画像の輝度平均値を初期の閾値（ｘ，ｙ）として２値化を行うと、隣接する領域の閾値が異なるとき、領域の境界で誤って粒子を抽出するという問題があるため、補間を行い滑らかに初期閾値を設定する。図１２に示すように、分割領域の中心座標（ｘ₁，ｙ₁）、（ｘ₁，ｙ₂）、（ｘ₂，ｙ₂）、（ｘ₂，ｙ₁）における初期閾値は、平均輝度Ｉ_avg1、Ｉ_avg2、Ｉ_avg3、Ｉ_avg4として、ｘ₁＜ｘ＜ｘ₂、ｙ₁＜ｙ＜ｙ₂の範囲において、座標（ｘ，ｙ）の初期閾値（ｘ，ｙ）は周り４点の平均輝度から補間を行う。２次元で、格子の目毎の双１次補間（bilinear interpolation）であり、その式は、
ｉ≡（ｘ−ｘ₁）／（ｘ₂−ｘ₁）
ｊ≡（ｙ−ｙ₁）／（ｙ₂−ｙ₁）
・・・（６）
となり、これらｉ、ｊで表されるパラメータを用いて、次式で表される。

Θ（ｘ，ｙ）＝（１−ｉ）（１−ｊ）Ｉ_avg1＋ｉ（１−ｊ）Ｉ_avg2
＋ｉｊＩ_avg3＋（１−ｉ）ｊＩ_avg4 ・・・（７）
双１次補間は、補間点が格子の目から目へ移動する際に、補間関数（ｘ，ｙ）の値は連続的に変化するため、粒子の検出の精度が向上する。その後、ラベリングを行い、２値化画像の粒子として認識された領域の面積Ｓを求め、設定した範囲（Ｓ₁ ＜Ｓ＜Ｓ₂）の値であれば粒子とみなし、閾値（ｘ，ｙ）は変更しない。Ｓ₁ より小さければ２値化された領域は粒子ではなくノイズであると見なして、その領域の閾値（ｘ，ｙ）を最大値に増加させる、逆に、Ｓ₂より大きければ、２値化された領域は複数の粒子を含むとして閾値（ｘ，ｙ）を増加させる。閾値設定、２値化、ラベリング、判定のサイクルを繰り返し、画面全領域で閾値の変化が行われなくなるまで繰り返す。

後記のように、本発明ではＰＴＶ−ＳＳ法を適用するため、粒子はストリークの形状で撮像されているので、個々の粒子面積は光ファイバ１本によって生じる輝度分布に比べ十分大きくなる。そのため、上記の動的２値化処理のパラメータである最小粒子面積Ｓ₁ を光ファイバ１本による輝度分布より大きく設定することで、光ファイバのパターンを粒子と認識しなくすることが可能である。

以上のファイバノイズの除去処理４３１と粒子検出処理４３２を経て、光ファイバのパターンであるファイバノイズが除去され、動的２値化法により個々の粒子に対して異なる閾値を設定して２値化された微小時間異なる２時刻の粒子画像（「参照画像」と「探索画像」）４１Ｌを用いて、ＰＩＶベクトルを用いたＰＴＶ−ＳＳ法４３３において、個々の粒子の移動量を求める。ただし、この際用いる粒子画像４１Ｌは、前記したように、半導体レーザの発光時間を比較的長くとった粒子軌跡画像を用いる。

以下、まず、粒子軌跡追跡法（ＰＴＶ−ＳＳ）を説明し、その際に、ＰＩＶ処理４２Ｌによって得られたＰＩＶベクトルを補助的なベクトル情報として利用してＰＴＶ−ＳＳ法で検出される２つの可能なベクトルの中の１つに限定する点を説明する。

ＰＴＶ−ＳＳ法は、前記のように、照明光源のレーザ発振装置１１として用いる半導体レーザの発光時間を比較的長くとり、粒子の流跡を記録する方法であり、粒子群ではなく個々に粒子の移動量を見積もる方法である。これは、流跡を記録することにより、必要以上に粒子像径を大きくとることなく、かつ、最小の空間分解能となる光ファイバ２５単繊維に対して十分な粒子像の情報を与えることができるためである。ＰＴＶ−ＳＳ法は、レーザ発振装置１１の発光時間を容易に制御することができる半導体レーザを用いることによって可能となる。

ここで、本発明で用いているＰＴＶ−ＳＳ法と従来のＰＴＶ法の比較を図１３に示す。図１３中、（ａ１）と（ａ２）はそれぞれの方法の参照画像と探索画像を示す図であり、（ｂ１）と（ｂ２）はそれぞれの方法の半導体レーザの発振タイミングを示す図である。また、ＰＴＶ−ＳＳ法の粒子の軌跡の部分のみに関する図を図１４に示す。

ＰＴＶ−ＳＳ法は、図１３（ａ１）に示すように、レーザ照射時間ｔ_dを粒子移動速度に対して長めに設定することで、参照画像中、探索画像中に図１４に示すように長軸と短軸を持つ移動軌跡を記録するようにする。これに対して、ＰＴＶ法では、図１３（ａ２）に示すように、レーザ照射時間を粒子移動速度に対して短いので、粒子像のみで移動軌跡は記録されない。

ここで、レーザ照射時間をｔ_d、レーザ照射間隔をΔｔ、移動軌跡の長軸、短軸の長さをそれぞれＬ_p1、Ｌ_pとすると、粒子の速度が次の撮影画像まで一定という条件の下に、次の式（８）で次の粒子位置を推定することができる。

Ｌ＝（Ｌ_p1−Ｌ_p）Δｔ／ｔ_d ・・・（８）
ここで、Ｌは、次の撮影画像までの推定移動距離である。この値が粒子の軌跡画像から推定できることによって、図１３（ａ２）、（ｂ２）のような従来の相互相関のみを用いたＰＴＶ法に比べて、効率的に探索窓画像を設定することができ、計算量の減少、つまり、ＰＴＶの計算の速度の向上が行われる。また、広い探索窓画像を開く図１３（ａ２）の場合に比べ、粒子の移動位置を推定することで、粒子の検出の誤対応が大幅に減少する。なお、このようにして設定された「探索窓画像」と「参照窓画像」との間で式（４）と式（５）に基づいて相互相関を計算することにより、個々の粒子の移動量が求まる。

このように、本発明においては、ＰＴＶ−ＳＳ法を適用することで、粒子の軌跡長さから次の時間の画像の粒子位置を推定することで粒子の誤対応を減少させることが可能となる。しかし、粒子の軌跡からのみの情報では、移動先が唯一に決まる訳ではなく、図１５（ａ）に示すように、相互に１８０°異なる２つのベクトルが移動候補となる。従来は、方向情報に関しては流動場によって人為的に設定値を与えることで唯一のベクトルとなるように処理していたが、様々な流動場への対応を行うにはこれは不十分である。

そこで、本発明においては、この方向情報を決定するために、ＰＩＶ処理４２Ｌによって得られたＰＩＶベクトルの方向を用いて、図１５（ｂ）のように粒子の軌跡からのみの情報では存在した２つのベクトルの可能性を、そのＰＩＶベクトルの方向に向いた１つのベクトルに限定するようにする。なお、図１５（ｂ）の背景の白い矢印群がＰＩＶベクトルであり、×印を付されたベクトルが可能性を否定されたベクトルであり、○印を付されたベクトルが可能性を肯定され、残るベクトルである。

図１６に、例示として、トレーサ粒子の混入したグリセリン溶液を回転円盤上に設置して形成した擬似的な流れを撮影した粒子軌跡画像４１Ｌから、上記のＰＩＶベクトルを用いたＰＴＶ−ＳＳ法４３３で算出した移動量検出結果を示すベクトルマップを示す。

このように、ＰＩＶベクトルを用いたＰＴＶ−ＳＳ法４３３において、ＰＩＶ処理４２Ｌによって得られたＰＩＶベクトルを補助的に用いて、人為的なパラメータの設定をしないで、様々な流動場において自動で粒子の移動方向と移動量を計測することができる。

ところで、ＰＴＶ−ＳＳ処理４３Ｌにおける相互相関法による移動量の計測精度は、ＣＣＤカメラ２４Ｌ、２４Ｒによって得られる画像が空間的に離散化されたものであることから、得られる移動量は整数値である。そのため、整数値以下の移動量は相互相関法では導き出すことができない。したがって、このような流速計においては、１ピクセル以下の粒子移動量（サブピクセル移動量）の評価が必要であり、この特性を十分に考慮する必要がある。一般的には、近似関数として正規分布関数を使用して計測結果をサブピクセル精度化することが多い。これは、相互相関に用いる粒子画像が正規分布関数と同じで、そのためそれらの相互相関値マップも同様にガウス分布関数の形をとるからである。しかしながら、本発明においては、ファイバスコープを通して粒子画像を撮像するため、粒子画像が劣化して粒子の輝度分布が正規分布関数の形状をなさない。図１８はファイバスコープを通した場合の粒子軌跡画像での相互相関マップの例を示す図である。ファイバスコープを通した粒子画像の場合は、光ファイバ２５のコアによる影響、クラッドによる影響を大きく受け、図１８のように光ファイバ２５の周波数成分が大幅に突出した相関値マップになっていることが分かる。このため、図１８の相関値マップに対してサブピクセル移動量算出を行う際に、正規分布関数を用いて近似を行うことは適正ではない。

そこで、本発明においては、ＰＴＶ−ＳＳ処理４３Ｌに続く光ファイババンドルによる画像劣化のサブピクセル処理４４Ｌにおいて、計測結果をサブピクセル精度化するための近似関数として参照画像の自己相関マップを使用する。参照画像の自己相関マップの例を図１７に示す。そして、算出された相互相関値マップ（図１９）に対して、２次曲線近似を行い、最終的なサブピクセル移動量を算出する。

ここで、１つの粒子のサブピクセル移動量算出手順をまとめると、以下のようになる。

１．参照画像の自己相関値をならべた自己相関値マップを作成する（図１７）。

２．参照画像と探索画像の相互相関値をならべた相互相関値マップを作成する（図１８）。

３．相互相関値マップ（図１８）においてサブピクセル精度で最大点を求めるため、自己相関値マップ（図１７）を近似関数として用いる。

４．自己相関値マップ（図１７）と相互相関値マップ（図１８）の相互相関値を並べた相互相関値マップ（図１９）を作成する。

５．算出された相互相関値マップ（図１９）に対して、２次曲線近似を行い、最終的なサブピクセル移動量を算出する。

図２０に、自己相関値マップと相互相関値マップの相互相関値から得られた相互相関値マップに２次曲線近似を行った一例を示す（ただし、この例は１次元方向の２次曲線近似である。）。

光ファイババンドルによる画像劣化のサブピクセル処理４４Ｌにおける以上の手順を踏むことで、粒子径や粒子の速度により粒子画像に違いが生じた場合においても、それぞれに対応した近似関数を設定することが可能になる。また、これによって光ファイバによる低解像度化とクラッドによる像の欠落がサブピクセル移動量算出に及ぼす影響を大幅に軽減することが可能となる。

以上のＰＩＶ処理４２Ｌ（ファイバノイズの除去処理４２１、相互相関ＰＩＶ処理４２２、誤ベクトル除去及び空間平均処理４２３）、ＰＴＶ−ＳＳ処理４３Ｌ（ファイバノイズの除去処理４３１、粒子検出処理４３２、ＰＩＶベクトルを用いたＰＴＶ−ＳＳ法４３３）、光ファイババンドルによる画像劣化のサブピクセル処理４４Ｌを経て、次の最終的なベクトルデータ処理４５Ｌにおいて、左のＣＣＤカメラ２４Ｌで撮像された左のカメラ画像４１Ｌから算出される２次元速度ベクトルデータ情報が得られる。

右のＣＣＤカメラ２４Ｒで撮像された右のカメラ画像４１Ｒについても、ＰＩＶ処理４２Ｌ（ファイバノイズの除去処理４２１、相互相関ＰＩＶ処理４２２、誤ベクトル除去及び空間平均処理４２３）、ＰＴＶ−ＳＳ処理４３Ｌ（ファイバノイズの除去処理４３１、粒子検出処理４３２、ＰＩＶベクトルを用いたＰＴＶ−ＳＳ法４３３）、光ファイババンドルによる画像劣化のサブピクセル処理４４Ｌと同様の処理を経て、次の最終的なベクトルデータ処理４５Ｒにおいて２次元速度ベクトルデータ情報が得られる。

これら左右のカメラ画像４１Ｌ、４１Ｒについて得られる２次元速度ベクトルデータ情報ぞれぞれを本発明に基づく流動計測結果としてもよいが、以上の実施例においては、左右視差を持った一対の２次元速度ベクトルデータ情報が得られているので、その一対の２次元速度ベクトルデータ情報から、次の３次元再構成処理４６により３次元速度ベクトルデータ情報を得るようにして、３次元性を有する流動場の３次元流動計測が可能になる。以下に、この３次元再構成処理４６について説明する。

ここで利用するのは、ステレオＰＴＶと呼ばれる方法であり、２台の撮像系（ＣＣＤカメラ２４Ｌ、２４Ｒ及びその光学系２７ａ、２７ｂ、２８等）と較正による３次元情報とを用いることによって、速度の第３方向成分の計測を可能にするものである。すなわち、レーザシート１７でシート状に照明された領域を貫く方向に存在している方向成分を計測し、３方向成分速度ベクトルの２次元的分布を同時に計測する手法である。そして、マッピング関数による実座標とカメラ内座標との対応付けを行うことにより、光学窓及び作動流体による屈折の影響を無視して計測することが可能となっており、したがって、湾曲した光学窓を有する流路等これまで計測が困難とされた系の計測を行うことも可能である。

ステレオＰＴＶは、図２１のように、較正によって実座標と２台のカメラのカメラ内座標（以下、ピクセル座標と称する。) との対応付けを行うことで、実座標系での速度の算出を行うものである。実座標とピクセル座標との対応付けは、レンズの収差、流体及び光学窓による屈折、遠近効果を伴うために、その関数を厳密に決定することは困難である。特に、ファイバスコープの対物レンズ２７ａはレンズ径が小さく、焦点距離も短い。焦点距離の短いレンズ系では、歪曲収差が無視できないものとなる。また、ＣＣＤカメラ２４Ｌ、２４Ｒの撮像面に対しレーザシート１７が傾きを持っているため、カメラ内座標と実座標との対応付けが重要になってくる。また、その較正によって得られたデータは、粒子の３成分位置の算出を行う上で非常に重要になってくる。しかし、実座標とピクセル座標との対応付けは、レンズの収差、流体及び光学窓による屈折、遠近効果を伴うために、その関数を厳密に決定することは困難である。そこで、本実施例では、マッピング関数と呼ばれる実座標とピクセル座標を対応付ける関数を用いる方法（非特許文献３）を採用する。

すなわち、マッピング関数は、実座標ｘ，ｙをピクセル座標値Ｘ，Ｙの多項式によって近似するものである。本実施例では、以下の式（９）に示されるような３次の多項式への近似を採用する。

式（９）において、Ｘ_i，Ｙ_iは較正によって得られるピクセル座標であり、ｘ_i，ｙ_iは実座標において対応する座標値（既知）である。また、ａ_i，ｂ_iは係数であり、較正によって算出する必要がある。ここで、ｎは較正に用いる点の数であり、ｎ＝１６のときに係数ａ_i，ｂ_iの数と一致し、逆行列の算出から係数を算出することが可能となる。しかし、実際に得られるＸ_i，Ｙ_iの位置はカメラの解像度に依存する不確かさを持っているので、これによって発生する対応づけ誤差を抑制するためにはより多くの較正点から最適な係数を決定することが必要となる。そこで、本実施例では、係数決定法として最小自乗法を用いる。最小自乗法による係数決定は、式（１０）に示される正規方程式から求める。

ｘ＝ＸＡに対して、
Ｘ^Tｘ＝Ｘ^TＸＡ・・・（１０）
ここで、Ｘは式（９）のピクセル座標値行列に、ｘが式（９）の実座標値行列に、Ａが式（９）の係数行列にそれぞれ対応する。正規方程式によって転置行列との積Ｘ^TＸは正方行列となり、その逆行列（Ｘ^TＸ）^-1は簡便に求めることが可能である。係数行列Ａを算出することにより、カメラ内座標を実座標へ変換することが式（９）を介して可能となる。

その逆に、実座標からカメラ内座標への変換は、係数行列の逆行列を算出することによって容易に行うことができる。また、２台のカメラ内座標の対応付け同様の方法によって行うことが可能であり、２台のカメラ内の座標と実座標とが相互に対応づけられる。

本実施例における具体的な較正手順について以下に説明する。

まず、図２２に示されるような１６点の格子点を持つキャリブレーションシート（白部分は光が透過）を印刷して作成する。本実施例においては、格子点間隔が２．５ｍｍのものを作成した。このキャリブレーションシートを、図２３に示すように３箇所（ｚ＝ａ，０，ｂ）それぞれの位置に移動させて、左右２台のカメラＡ、Ｂによりそのキャリブレーションシートの画像を撮像する。撮影されたピクセル画像を元に１６点の格子点の位置を検出する。検出の際には、粒子検出処理４３２で説明した粒子検出アルゴリズムを用いる。ピクセル座標を算出した後に、式（９）、式（１０）より左右のカメラにより３箇所ずつ撮像された計６枚のピクセル画像に対しての係数行列を求める。ここで、カメラＡにおけるｚ＝ａの位置での係数行列をＡ_a、ｚ＝ｂにおける係数行列をＡ_b、ｚ＝０における係数行列をＡ_c、カメラＢにおけるｚ＝ａの位置での係数行列をＢ_a、ｚ＝ｂにおける係数行列をＢ_b、ｚ＝０における係数行列をＢ_cとする。

このようにして算出された変換行列を元に、実際の３次元座標を算出していく。その手法を次に説明する。

上記の較正方法を利用して、粒子の３次元実座標を算出する方法を説明する。図２４は、粒子の実座標の算出の概念図であり、レーザシートによって照明された領域を上方から見た図であり、まず片方のカメラに関して、ｚ＝ａ，ｂの左右における２枚の較正行列に対して式（９）により実座標を算出する。これにより、２点の実空間座標（ｘ_a，ｙ_a，ａ），（ｘ_b，ｙ_b，ｂ）が求まる。これを結ぶ直線の式は、（ｘ−ｘ_a）／（ｘ_a−ｘ_b）＝（ｙ−ｙ_a）／（ｙ_a−ｙ_b）
＝（ｚ−ａ）／（ａ−ｂ）・・・（１１）
となる。同様にしてもう一方のカメラに関して直線の式を求め、この２直線が交わる位置が実際の粒子の位置となる。しかし、実際にはカメラの解像度が有限であることに起因する誤差、２次元ベクトルに内在する誤差、ノイズの影響等によって２直線は厳密には交差しない。本実施例では、２直線間の距離が最小となる２直線上の点の中点を以って交差点と見なして処理を行う。

しかし、この処理を行うに当たっては、まず左右のカメラにおける粒子の対応付けを行わなくてはならない。そこで、実際には、まず粒子の対応候補を探すために、全ての粒子がｚ＝０に存在すると仮定し、ｚ＝０における係数行列Ａ_c、Ｂ_cを用いて実座標に変換をする。変換された実座標値は粒子がｚ＝０と仮定した位置であるので、左右のカメラの粒子の実座標値は必ずしも一致しない。そこで、左カメラのある粒子の周囲に存在する右カメラの粒子の候補を求め、その候補全てに対して図２４の２直線の対応付け処理を行う。この処理を行い、最も直線間距離が短い場合を同一粒子とみなし、粒子の対応付け処理を行う。

３次元再構成処理４６においては、最終的なベクトルデータ処理４５Ｌ、４５Ｒにおいて、左右のＣＣＤカメラ２４Ｌ、２４Ｒで同期して撮像され、ＰＴＶ−ＳＳ処理４３３で対応づけられた個々の粒子の相互に微小時間異なる第１時刻と第２時刻における実座標値を上記の手法で求めることにより、３次元速度ベクトルデータ情報が得られる。

以上の本発明の流動計測方法の実流動場への適用性を確認するために、図２５に示すような配置で、軸対称噴流（出口流速２１４ｍｍ／ｓ，Ｒe ＝５３００）の高さ２Ｄ（Ｄ＝ノズル直径) における速度計測実験を行った。第３方向（ｚ方向）成分が検出されているかを検証するために、流動場に対して約１５°プローブを傾け計測した結果の空間平均ベクトルは図２６、図２７に示すようなものであり、平均１３．３°の傾きが検出された。

また、次の表１は、従来（例えば非特許文献１）と本発明の計測方法の手法を用いて５０組１００枚の同一画像を処理した結果、検出された３次元ベクトル数の関係である。従来の手法においては、ｘ軸を中心としてレーザーシートを回転させた場合、つまり、レーザーシートの面を貫く方向に流動場が発生している場合においてのベクトル算出数が約１／３に減少していたが、本発明の方法を用いた場合は、そのような減少は見られなかった。また、ベクトル検出数においても、平均して約２５倍程度の増加が見られた。

表１ ┌────┬──────┬──────┬───────┐ │傾斜角度│旧手法（本）│本手法（本）│検出割合（倍）│ ├────┼──────┼──────┼───────┤ │ ０°│ ２８９ │ ３５８３ │ １２．４０ │ ├────┼──────┼──────┼───────┤ │ １５°│ ９９ │ ３８４４ │ ３８．８３ │ ├────┼──────┼──────┼───────┤ │ ３０°│ １２７ │ ３９６７ │ ３１．２４ │ └────┴──────┴──────┴───────┘ 。

以上のように、本発明によると、より信頼性が高くより精度が高く精密な測定が可能な流体の流動計測システム及びその計測方法が達成できた。

本発明の流体の流動計測システムの１実施例の概要を示す図である。レーザシート形成用光学系のレーザシートのシート面に沿う断面とそれに直交する断面を示す図である。左右の画像伝送手段の先端のプローブ部の斜視図と断面図である。画像伝送手段の構成を説明するための図である。本発明の１実施例の流体の流動計測システムの画像処理ソフトウェアを説明するための図である。相互相関ＰＩＶ処理における粒子群の誤対応の例を示す図である。本発明のＰＩＶ処理で相互相関計算を行って得られた流速ベクトル（ＰＩＶベクトル）分布の例を示す図である。ファイバノイズの除去処理によって光ファイバ配列の輝度パターンが除去され、粒子像のみが残る画像処理の流れを任意のラインの輝度値分布で示した図である。光ファイババンドルによる画像の低解像度化による粒子の誤対応を説明するための図である。実際に撮影したイメージガイドの画像の任意ラインの輝度値分布と周波数解析の例を示す図である。動的２値化法の概要を示す図である。動的２値化法における画像の分割と補間を説明するための図である。本発明で用いているＰＴＶ−ＳＳ法と従来のＰＴＶ法の比較を示す図である。ＰＴＶ−ＳＳ法の粒子の軌跡の部分のみに関する図である。本発明のＰＩＶベクトルを用いたＰＴＶ−ＳＳ法でのベクトルの可能性とその選択を示す図である。擬似的な流れを撮影した粒子軌跡画像から本発明の方法で算出した移動量検出結果の例を示すベクトルマップである。参照画像の自己相関値をならべた自己相関値マップの例を示す図である。参照画像と探索画像の相互相関値をならべた相互相関値マップの例を示す図である。図１７の自己相関値マップと図１８の相互相関値マップの相互相関値を並べた相互相関値マップの例を示す図である。自己相関値マップと相互相関値マップの相互相関値から得られた相互相関値マップに２次曲線近似を行った一例を示す図である。ステレオＰＴＶにおいて較正によって実座標と２台のカメラのカメラ内座標との対応付けを行うことで実座標系での速度の算出を行う過程を模式的に示す図である。キャリブレーションシートの例を示す図である。図２２のキャリブレーションシートを位置に移動させて左右２台のカメラによりそのキャリブレーションシートの画像を撮像する様子を示す図である。図２３の較正方法を利用して粒子の３次元実座標を算出する方法を説明するための図である。本発明の流動計測方法の実流動場への適用性を確認するための配置例を示す図である。図２５の配置例で計測した結果の空間平均ベクトル分布を示す図である。図２５の配置例で計測した結果の空間平均ベクトル分布を示す別の図である。

符号の説明

１０…流動計測システム
１１…レーザ発振装置
１２…送光用光ファイバ
１３…レーザシート形成用光学系
１４…流体流動場
１７…レーザシート
１７ａ…測定範囲
１８…レンズ群
１９ａ、１９ｂ、１９ｃ…平凸レンズ
１９ｄ…平凹シリンドリカルレンズ
２２、２２Ｌ、２２Ｒ…画像伝送手段
２４、２４Ｌ、２４Ｒ…ＣＣＤカメラ
２５…光ファイバ
２６…イメージガイド
２７ａ…対物レンズ
２７ｂ…４５°直角プリズム
２８…カメラレンズ
２９Ｌ、２９Ｒ…画像撮像手段
３０…シグナルジェネレータ
３１…パソコン
３２…レーザドライバ
３３…ペルチエドライバ
３４Ｌ、３４Ｒ…インターフェース
４１Ｌ、４１Ｒ…カメラ画像
４２Ｌ…ＰＩＶ処理
４２１…ファイバノイズの除去処理
４２２…相互相関ＰＩＶ処理
４２３…誤ベクトル除去及び空間平均処理
４３Ｌ…ＰＴＶ−ＳＳ処理
４３１…ファイバノイズの除去処理
４３２…粒子検出処理
４３３…ＰＩＶベクトルを用いたＰＴＶ−ＳＳ法
４４Ｌ…光ファイババンドルによる画像劣化のサブピクセル処理
４５Ｌ、４５Ｒ…最終的なベクトルデータ処理
４６…３次元再構成処理

Claims

レーザ光を発振させるレーザ発振装置と、発振されたレーザ光を流体の流動場内にシート状に投入させるレーザシート形成用光学系と、このレーザシート形成用光学系からのレーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を選択的に撮像する画像撮像手段と、前記画像撮像手段で撮像された微小時間異なる２時刻における２次元粒子画像の輝度パターンを比較・解析して粒子群の移動方向及び移動量を計測する第１画像処理手段と、前記画像撮像手段で撮像された微小時間異なる２時刻における２次元粒子軌跡画像の輝度パターンを比較・解析して個々の粒子の移動方向及び移動量を計測する第２画像処理手段とを備え、前記第２画像処理手段では、前記第１画像処理手段で計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めるように構成されていることを特徴とする流体の流動計測システム。
前記レーザ発振装置は半導体レーザで構成され、前記半導体レーザの発光時間の差によって２次元粒子画像と２次元粒子軌跡画像を形成することを特徴とする請求項１記載の流体の流動計測システム。
前記レーザシート形成用光学系からのレーザシートの両側から同時に２次元粒子画像と２次元粒子軌跡画像を選択的に撮像する２つの画像撮像手段を備え、各々の前記画像撮像手段に対してそれぞれ前記第１画像処理手段と、前記第２画像処理手段とを備え、各々の前記第２画像処理手段では、各々の前記第１画像処理手段で計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めるように構成されていることを特徴とする請求項１又は２記載の流体の流動計測システム。
前記第２画像処理手段各々で計測された個々の粒子の移動方向及び移動量を相互に対応付けて個々の粒子の３次元的な移動方向及び移動量に再構成する３次元再構成手段を有し、個々の粒子の３次元的な移動方向及び移動量を計測することを特徴とする請求項３記載の流体の流動計測システム。
前記画像撮像手段は、前記レーザシート形成用光学系からのレーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を光伝送する画像伝送手段と、光伝送された２次元粒子画像と２次元粒子軌跡画像を撮像する撮像手段とを備えていることを特徴とする請求項１から４の何れか１項記載の流体の流動計測システム。
前記画像伝送手段は、多数本の光ファイバを束ねて一体化させ、両端面が平面加工された可撓性のイメージガイドで構成され、前記撮像手段は、ＣＣＤカメラで構成されており、前記イメージガイドは、対物レンズにより一方のファイバ端面に結像された画像を光ファイバで各画素に分解してＣＣＤカメラ側の他端面まで同一画像として伝送するように構成されていることを特徴とする請求項５記載の流体の流動計測システム。
前記レーザシート形成用光学系と前記画像撮像手段とは、一体的に組み合わせて構成されていることを特徴とする請求項１から６の何れか１項記載の流体の流動計測システム。
前記第１画像処理手段、前記第２画像処理手段はそれぞれ２次元粒子画像と２次元粒子軌跡画像から前記イメージガイドに基づくファイバノイズを除去するファイバノイズ除去手段を備えていることを特徴とする請求項６又は７記載の流体の流動計測システム。
前記第２画像処理手段で計測された個々の粒子の移動方向及び移動量をサブピクセル精度化するためのサブピクセル処理手段を備えていることを特徴とする請求項１から８の何れか１項記載の流体の流動計測システム。
レーザ発振装置から発振されるレーザ光を流体の流動場にシート状に照射してレーザシートを形成し、レーザシートで照明されて形成される２次元粒子画像と２次元粒子軌跡画像を選択的に画像撮像手段で撮像し、撮像された微小時間異なる２時刻における２次元粒子画像の輝度パターンを比較・解析して粒子群の移動方向及び移動量を計測し、また、撮像された微小時間異なる２時刻における２次元粒子軌跡画像の輝度パターンを比較・解析して個々の粒子の移動方向及び移動量を計測し、その個々の粒子の移動方向及び移動量の計測の際に、先に計測された粒子群の移動方向に基づいて個々の粒子の移動方向を定めることにより、流体流動場における流体の流動を計測することを特徴とする流体の流動計測方法。