(57)【要約】【課題】流体の流れ解析を実環境下で容易に行う。【解決手段】流体を撮影するカメラ１１と、カメラ１１で所定期間撮影した流体の温度画像データから一定期間の温度データを、一定の解析ステップずつ時間方向にずらして切り出して部分温度データを生成する（ステップＳ１６）。そして、切り出した部分温度データ毎に解析を行い、一定期間毎の温度変動分布を検出する（ステップＳ１７、Ｓ１８）。【選択図】図４
2020052036.tif 000002
【特許請求の範囲】
【請求項１】
流体を撮影するカメラと、
当該カメラで所定期間撮影した前記流体の温度画像データから一定期間の温度データを、一定の解析ステップずつ時間方向にずらして切り出して部分温度データを生成する切り出し処理部と、
前記切り出し処理部で切り出した部分温度データ毎に解析を行い、前記一定期間毎の温度変動分布を検出する解析部と、
を備えることを特徴とする流体の流れ計測装置。
【請求項２】
前記解析部は、前記部分温度データ毎にフーリエ変換処理を行うことを特徴とする請求項１に記載の流体の流れ計測装置。
【請求項３】
前記フーリエ変換処理は、コサイン変換及びサイン変換のうちの少なくともいずれか一方であることを特徴とする請求項２に記載の流体の流れ計測装置。
【請求項４】
前記フーリエ変換処理の解析周波数を設定する解析周波数設定部を備え、
前記解析部は、前記解析周波数設定部で設定された解析周波数の基底関数を用いて前記フーリエ変換処理を行うことを特徴とする請求項２または請求項３に記載の流体の流れ計測装置。
【請求項５】
前記解析周波数設定部で設定された解析周波数設定値を、前記流体と対応付けて記憶する記憶部を備え、
前記解析部は、前記記憶部に記憶されている流体について前記フーリエ変換処理を行うときには、前記記憶部に記憶されている前記解析周波数設定値を用いて前記フーリエ変換処理を行うことを特徴とする請求項４に記載の流体の流れ計測装置。
【請求項６】
前記一定期間の期間幅と、前記解析ステップのステップ幅の少なくとも一方を設定するパラメータ設定部を備え、
前記切り出し処理部は、前記パラメータ設定部で設定された設定値を用いて前記部分温度データを生成することを特徴とする請求項１から請求項５のいずれか一項に記載の流体の流れ計測装置。
【請求項７】
前記パラメータ設定部で設定された設定値を、前記流体と対応付けて記憶する記憶部を備え、
前記切り出し処理部は、前記記憶部に記憶されている流体について前記温度データから前記部分温度データを切り出すときには、前記記憶部に記憶されている前記設定値を用いて切り出すことを特徴とする請求項６に記載の流体の流れ計測装置。
【請求項８】
表示装置と、
前記解析部で検出した前記温度変動分布を前記表示装置に表示する表示処理部と、
を備え、
前記表示装置と、前記カメラと、前記切り出し処理部と、前記解析部と、前記表示処理部とは、一体に形成されている請求項１から請求項７のいずれか一項に記載の流体の流れ計測装置。
【請求項９】
前記カメラで撮影した前記流体の温度画像データに対してデジタルフィルタ処理を行うフィルタ処理部と、
当該フィルタ処理部での処理後の温度画像データで表される画像を、温度変動分布を表す画像として表示装置に表示する表示処理部と、
前記解析部により前記温度変動分布を検出する第一の検出処理と、前記フィルタ処理部で前記デジタルフィルタ処理を行うことにより前記温度変動分布を検出する第二の検出処理と、のうちのいずれかを選択する選択部と、
を備え、
前記選択部で選択された検出処理により前記温度変動分布を検出することを特徴とする請求項１から請求項７のいずれか一項に記載の流体の流れ計測装置。
【請求項１０】
表示装置を備え、
前記表示装置と、前記カメラと、前記切り出し処理部と、前記解析部と、前記表示処理部と、前記フィルタ処理部と、前記選択部とは、一体に形成されている請求項９に記載の流体の流れ計測装置。
【請求項１１】
流体を撮影するカメラと、
当該カメラで撮影した前記流体の温度画像データに対してデジタルフィルタ処理を行うフィルタ処理部と、
当該フィルタ処理部で処理後の温度画像データで表される画像を、温度変動分布を表す画像として表示装置に表示する表示処理部と、
を備えることを特徴とする流体の流れ計測装置。
【請求項１２】
前記フィルタ処理部は、前記カメラから前記温度画像データが入力される毎に前記デジタルフィルタ処理を行い、
前記表示処理部は、前記フィルタ処理部で前記デジタルフィルタ処理が行われる毎に、処理後の温度画像データで表される画像を表示することを特徴とする請求項１０又は請求項１１に記載の流体の流れ計測装置。
【請求項１３】
前記フィルタ処理部は、ＩＩＲフィルタを用いて前記デジタルフィルタ処理を行うことを特徴とする請求項１２に記載の流体の流れ計測装置。
【請求項１４】
表示装置と、
前記カメラと、
前記フィルタ処理部と、
前記表示処理部とは、一体に形成されている請求項１１から請求項１３のいずれか一項に記載の流体の流れ計測装置。
【請求項１５】
前記デジタルフィルタ処理はバンドパスフィルタ処理であり、中心周波数は選択可能であることを特徴する請求項９から請求項１４のいずれか一項に記載の流体の流れ計測装置。
【請求項１６】
前記表示処理部は、
異なる時刻における前記温度変動分布から、互いに対応する箇所を、温度変動の対応点として検出し、
前記異なる時刻間における前記対応点の移動量に基づき、前記対応点の移動速度をベクトルで表示することを特徴とする請求項８から請求項１５のいずれか一項に記載の流体の流れ計測装置。
【請求項１７】
前記流体は気体であり、前記カメラは赤外線カメラであって、
前記赤外線カメラの入力部に、前記気体に対応した通過帯域を有するバンドパスフィルタ処理を行う赤外線バンドパスフィルタを備える請求項１から請求項１６のいずれか一項に記載の流体の流れ計測装置。
【請求項１８】
前記カメラは、視野が重複するように間隔を空けて設けられた２台のカメラを含み、
当該カメラ毎に、その温度画像データから前記温度変動分布を取得するようになっており、
前記カメラ毎に取得した前記温度変動分布それぞれから、同一箇所に対応する温度変動分布を抽出し、抽出した温度変動分布から、ステレオ立体視の原理を用いて前記抽出した温度変動分布に対応する部分の、実際の奥行方向の位置を演算する奥行情報算出部を、さらに備える請求項１から請求項１７のいずれか一項に記載の流体の流れ計測装置。
【請求項１９】
請求項１から請求項１８のいずれか一項に記載の流体の流れ計測装置を有するガス検知装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、流体の流れ計測装置及びガス検知装置に関する。
【背景技術】
【０００２】
従来、流体の流れの解析には、粒子画像流速計測法（ＰＩＶ  Particle Image Velocimetry）が広く用いられている（例えば、特許文献１参照。）。ＰＩＶ法は、トレーサ粒子と呼ばれる粒子を用いるものであり、解析の対象の流体内にトレーサ粒子を供給し、このトレーサ粒子をレーザ光で照射し、レーザ光の発振と同期して撮影することで、流体の流れを解析するようにしたものである。
また、流体自体を異なる二つの色で着色し、着色された流体を交互に流すことで、流体の流れを観察する方法等も提案されている（例えば、特許文献２参照。）。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１８−４０５７１号公報
【特許文献２】特開２００６−７１３８２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ところで、ＰＩＶ法においては、レーザ照射装置等の大がかりな装置を必要とし、且つ人の立ち入りが制限された環境下で行う必要がある。また、トレーサ粒子を供給する必要があることから、煩雑であり、実験室内での測定に限られていた。さらに、トレーサ粒子を測定する必要があることから、１回に測定できる時間も限られていた。
そのため、広い範囲を対象とする建物内や屋外、工場内、人が大勢いる駅等の公共の設備等、実際の構造物や環境に対しては、流体の流れの解析を有効に行うことは困難であった。
【０００５】
そこで、この発明は上記従来の未解決の問題に着目してなされたものであり、簡易であり且つ実環境下での流れ解析を容易に行うことの可能な流体の流れ計測装置及びガス検知装置を提供することを目的としている。
【課題を解決するための手段】
【０００６】
本発明の一態様によれば、流体を撮影するカメラと、カメラで所定期間撮影した流体の温度画像データから一定期間の温度データを、一定の解析ステップずつ時間方向にずらして切り出して部分温度データを生成する切り出し処理部と、切り出し処理部で切り出した部分温度データ毎に解析を行い、一定期間毎の温度変動分布を検出する解析部と、を備える、流体の流れ計測装置が提供される。
【０００７】
また、本発明の他の態様によれば、流体を撮影するカメラと、カメラで撮影した流体の温度画像データに対してデジタルフィルタ処理を行うフィルタ処理部と、フィルタ処理部で処理後の温度画像データで表される画像を、温度変動分布を表す画像として表示装置に表示する表示処理部と、を備える、流体の流れ計測装置が提供される。
【０００８】
さらに、本発明の他の態様によれば、上記いずれかの態様に記載の流体の流れ計測装置を有するガス検知装置が提供される。
【発明の効果】
【０００９】
本発明の一態様によれば、流体の流れ計測を簡易な構成で、且つ実環境下で容易に行うことができる。
【図面の簡単な説明】
【００１０】
【図１】本発明の第一実施形態に係る流体の流れ計測装置の一例を示す構成図である。
【図２】切り出し窓幅及び解析ステップを説明するための説明図である。
【図３】撮影時の処理手順の一例を示すフローチャートである。
【図４】温度変動分布解析時の処理手順の一例を示すフローチャートである。
【図５】切り出し窓幅の最適化を説明するための説明図である。
【図６】本発明の第二実施形態にかかるバンドパスフィルタ法による温度変動分布解析時の処理手順の一例を示すフローチャートである。
【図７】本発明の第三実施形態に係る流体の流れ計測装置の一例を示す構成図である。
【図８】表示装置を備えた赤外線カメラの一例を示す概略構成図である。
【図９】計測対象の気流が存在する領域の可視画像の一例である。
【図１０】計測対象の気流が存在する領域の温度画像の一例である。
【図１１】時系列に並べた温度変動分布画像の一例である。
【図１２】窓幅（時間分解能）と周波数分解能との関係を説明するための説明図である。
【図１３】エアダスターが噴射する気流が存在する領域の温度画像の一例である。
【図１４】時系列に並べた温度変動分布画像の一例である。
【図１５】フィルタホイールを備えたカメラの一例を示す概略構成図である。
【図１６】流れの速度ベクトルの演算方法を説明するための図である。
【図１７】流れの速度ベクトルの表示例である。
【図１８】本発明のその他の実施例を示す概略構成図である。
【図１９】本発明のその他の実施例を示す概略構成図である。
【図２０】本発明のその他の実施例を説明するための速度ベクトルの表示例である。
【発明を実施するための形態】
【００１１】
以下、図面を参照して本発明の実施形態について説明する。
なお、以下の詳細な説明では、本発明の実施形態の完全な理解を提供するように多くの特定の具体的な構成について記載されている。しかしながら、このような特定の具体的な構成に限定されることなく他の実施態様が実施できることは明らかである。また、以下の実施形態は、特許請求の範囲に係る発明を限定するものではない。また、実施形態の中で説明されている特徴の組み合わせの全てが発明の解決手段に必須であるとは限らない。
【００１２】
まず、本発明の第一実施形態を説明する。
図１は本発明の第一実施形態に係る流体の流れ計測装置の一例を示す構成図であって、流体の流れ計測装置１は、図1に示すように、カメラ１１と、入力装置及び表示装置を有する画像処理装置１２と、カメラ１１からの温度画像データを格納するためのハードディスク等からなる記憶装置（特許請求の範囲に記載の記憶部に対応）１３とを備える。流体の流れ計測装置１では、予め設定した所定期間の間、予め設定した周期でカメラ１１によって計測対象の気流が存在する領域を撮影し、温度画像データを取得する。カメラ１１の温度画像データは画像処理装置１２を介して時系列に記憶装置１３に格納される。画像処理装置１２は、記憶装置１３に格納された時系列の温度画像データ及び、予め記憶装置１３に格納されている温度変動分布解析に必要なデータベースをもとに、流体の流れ解析を行う。
カメラ１１としては、一般には赤外線カメラが好適に用いられるが、計測対象の流体の温度が高い場合（例えば、７００℃以上）には、可視カメラも使用することができる。
【００１３】
画像処理装置１２は、図２(ａ)に示すように、時系列温度データの一部を、一定期間の切り出し窓幅Δｔ１（特許請求の範囲に記載の一定期間に対応）の単位で切り出し、切り出した温度データを部分温度データとする。このとき、時系列温度データから部分温度データを切り出す時刻は、図２（ｂ）に示すように一定の解析ステップずつずらすものとする。これにより、切り出し窓幅相当分の期間を有する温度データからなる部分温度データが複数切り出される。
なお、解析される時系列温度データの温度変動に欠落が生じることを防ぐため、連続する部分温度データどうしは、温度データの一部が重なるように切り出されることが望ましい。画像処理装置１２は、切り出した部分温度データそれぞれに対してフーリエ変換を行う。
【００１４】
切り出し窓幅とは、図２（ｂ）に示すように、部分温度データとして切り出す期間を規定するものであり、例えばΔｔ１期間分の部分温度データを抽出する場合には、Δｔ１が切り出し窓幅となる。また、解析ステップとは、図２（ｂ）に示すように、部分温度データを抽出する際の切り出し窓のずれ量を表し、ずれ量がΔｔ２である場合には、Δｔ２ずつずれた時点における、Δｔ１期間分の部分温度データが抽出される。
画像処理装置１２は、切り出した部分温度データ毎に周波数解析法により、温度変動分布解析を行い、解析結果を表示装置に時系列に表示する。なお、本実施形態においては、時系列の温度データに対し、短時間の窓関数を順次ずらしながら掛けて切り出し、それに対してフーリエ変換を行っており、この方法を短時間フーリエ変換法という。短時間フーリエ変換法を用いることにより、時間変化に伴う温度変動を得ることができる。
【００１５】
一般に、流体の温度は一様ではなくばらつきがある。特に気体の流れを知りたいような対象、例えば、エアーコンディショナーから送出される空気の流れ、車両の排気管から排出される排気ガスの流れ、ドライヤから送出される熱風の流れ等、温度が大きく異なる気体が検知対象として挙げられることが多い。
ここで、計測された温度範囲に比べて温度変動量がわずかな場合、温度画像データそのものからは流れを認識することはできない。例えば、計測する温度範囲を１５℃から２５℃とし、温度変動を０．１℃とすると、温度画像データを表示する場合、ダイナミックレンジに対する変化量は１／１００となるため、温度変動成分は、背景に埋もれてしまい、流れとして表現することはできない。つまり、温度画像データを時系列に表示したとしても、温度画像データから温度変動成分を認識することは困難である。
【００１６】
そのため、一定期間内での解析により温度変動成分を算出し、これを指標とすることで、表示の温度範囲を狭めることができ、その結果、温度変動成分を、背景温度に埋もれることなく流体の流れとして表現することができる。特に背景の温度範囲が広い場合、つまり、ダイナミックレンジが広い場合に、効果的である。
記憶装置１３に格納されているデータベースは、例えば表１に示すように、計測対象の流体を特定する情報と、計測対象毎の切り出し窓幅、解析ステップ及び解析周波数等を含む。
【００１７】
【表１】
2020052036.tif 000003
【００１８】
このデータベースは、計測対象の流体に対して温度変動分布解析を行ったときに、新規登録または更新設定される。
画像処理装置１２では、計測対象の流体がデータベースに登録されている場合には、登録された切り出し窓幅、解析ステップ、解析周波数を用いて温度変動分布解析を行う。また、計測対象の流体がデータベースに登録されていないとき、また、ユーザの操作により入力装置からこれらパラメータの値が入力設定されたときには、設定された各種パラメータ値を用いて温度変動分布解析を行う。
【００１９】
図３は、計測対象の流体を撮影する時の画像処理装置１２で実行される処理の処理手順の一例を示すフローチャートである。
画像処理装置１２では、カメラ１１により、予め設定されたフレームレートで計測対象の流体を撮影し（ステップＳ１）、画像処理装置１２では、カメラ１１で撮影した温度画像データを記憶装置１３に時系列に格納する（ステップＳ２）。
【００２０】
図４は、温度変動分布解析を行うときの、画像処理装置１２で実行される処理の処理手順の一例を示すフローチャートである。
画像処理装置１２は、記憶装置１３に格納された温度画像データを読み込み（ステップＳ１１）、続いて入力装置により計測対象の流体を特定する情報が設定されると、画像処理装置１２は、設定された計測対象の流体が記憶装置１３のデータベースに登録されていれば（ステップＳ１２、Ｓ１３）、対応する、切り出し窓幅、解析ステップ及び解析周波数を読み出し、これらを温度変動分布解析用のパラメータ値として設定する（ステップＳ１４）。
【００２１】
一方、計測対象の流体がデータベースに登録されていなければ（ステップＳ１３）、画像処理装置１２は、入力装置により設定された切り出し窓幅、解析ステップ及び解析周波数を、温度変動分布解析用のパラメータ値として設定する（ステップＳ１５  特許請求の範囲のパラメータ設定部、解析周波数設定部に対応）。
次いで、画像処理装置１２は、ステップＳ１４またはステップＳ１５で設定された切り出し窓幅と解析ステップで、時系列温度データから部分温度データを切り出す（ステップＳ１６  特許請求の範囲に記載の切り出し処理部に対応）。
そして、切り出した部分温度データそれぞれについて短時間フーリエ変換法を用いて温度変動分布解析を行う（ステップＳ１７  特許請求の範囲に記載の解析部に対応）。具体的には、温度変動分布解析用のパラメータ値として設定された解析周波数の正弦波信号と余弦波信号を用いて、切り出した部分温度データ信号に対し、サイン変換とコサイン変換からなるフーリエ変換を行う。
【００２２】
このとき、入力装置からパラメータ値が再度設定されたならばステップＳ１５からステップＳ１６の処理を再度行い、例えば確定操作が行われたとき、このときのパラメータ値を温度変動分布解析用のパラメータ値としてデータベースに記憶する。さらに、部分温度データ毎の短時間フーリエ変換による温度変動分布解析の結果について、振幅を温度変動値とし、位相から温度変動の正負を判断する。
なお、ここでは、正弦波信号と余弦波信号とを用いてフーリエ変換を行っているが、正弦波信号及び余弦波信号のうちのいずれか一方のみを用いた変換処理を行ってもよい。
【００２３】
そして、記憶した温度変動分布解析用のパラメータ値を用いて、一つの画素の温度変動分布解析を行ったならば、次の画素について同様に処理を行う。そして、全ての画素について温度変動分布解析を行ったならば、次の部分温度データについて同様に温度変動分布解析を行い、全ての部分温度データについて温度変動分布解析を行ったならば、ステップＳ１８に移行し、解析結果の表示処理を行う。具体的には、部分温度データの単位で、温度変動分布解析の結果を時系列に表示する。
【００２４】
ステップＳ１８での、解析結果の表示処理は、例えば、後述の図１１に示すように、カメラ１１による温度画像における各画素について同一時点における部分温度データから得た、各画素の温度変動値からなる温度変動分布画像を、温度変動値の大きさに応じて表示形態を変える等して、表示装置に表示する。図１１に示すように、各部分温度データから得た温度変動分布画像を時系列に表示することによって、温度変動分布の変化状況を表示することができる。なお、図１１は、気流を計測対象としたものであり、天井における気流の流れを可視化したものである。
ここで、気流があるということは、ある一つの地点では、温度変動が現れることになる。そのため、温度変動の分布を計測し、温度変動の分布を時系列に表示することによって、計測対象の気流の流れを可視化することができる。つまり、流れがある場合、単位時間に温度変動が凸になる場合と温度変動が凹になる場合とが、交互に繰り返されるが、これを時間的に連続して描画することで、気流の流れを可視化することができる。
【００２５】
次に、第一実施形態の動作の一例を説明する。
まず、計測対象の流体の温度画像データを取得する。そのためユーザは、計測対象の流体が視野内に含まれるように、カメラ１１を固定する。そして、カメラ１１により予め設定したフレームレートで所定期間継続して撮影を行い、所定期間分の温度画像データを得る。この温度画像データは画像処理装置１２を介して記憶装置１３に時系列に格納される。
次に、記憶装置１３の時系列の温度画像データを画像処理装置１２に取り込む。
【００２６】
次に、ユーザは、計測対象の流体を特定する情報を入力し、初めて解析を行う流体である場合には、切り出し窓幅、解析ステップ、解析周波数を入力設定する。
画像処理装置１２は、記憶装置１３から温度画像データを読み込むと、設定された解析ステップだけ切り出し領域をずらしながら、設定された切り出し窓幅で切り出して、部分温度データを得る。
そして、設定された解析周波数を用いて、フーリエ変換処理を行う。ここで、解析周波数や、切り出し窓幅、解析ステップを調整する。
【００２７】
つまり、瞬間的に発生する温度変動分布解析を行うためには、図５（ａ）に示すように、発生する温度変動の時間周期と同等程度の切り出し窓幅で部分温度データを切り出すことが好ましい。すなわち、図５（ｂ）に示すように、温度変動の時間周期に比較して切り出し窓幅が広すぎると、解析の時間分解能が低下し、流れの時間変化を十分にとらえきれない可能性が生じる。逆に図５（ｃ）に示すように、温度変動の時間周期に対して切り出し窓幅が狭すぎると解析のためのデータ数が少ないため、ノイズが重畳すること、また、発生している温度変動の一部のみを解析することになることにより温度変動値が低く見積もられる可能性がある。そのため、発生する温度変動の時間周期と同等程度の切り出し窓幅で部分温度データを切り出すことができるように、各種パラメータを探索する。
【００２８】
そして、良好な部分温度データを切り出すことができたとみなすことができるときの各種パラメータを用いて得た温度変動値からなる温度変動分布画像を、解析結果として表示する。また、良好な部分温度データを切り出すことができたときの解析周波数、切り出し窓幅、及び解析ステップの各種パラメータ値を、計測対象の流体、例えば、空気、ガス、排気ガス等といった計測対象の流体と対応付けてデータベースに新規に登録または更新記憶する。
【００２９】
この計測対象の流体に対して、別の時点で取得した温度画像データに対して温度変動分布解析を行う場合には、データベースに記憶されている計測対象の流体のパラメータを読み出し、読み出した切り出し窓幅、解析ステップで温度画像データから抽出した時系列温度データから部分温度データを切り出し、部分温度データに対して読み出した解析周波数を用いて温度変動分布解析を行い、得られた解析画像を表示する。
ここで、同一の計測対象の流体における温度変動の周波数は、同等程度となる傾向にある。したがって、データベースに格納されたパラメータを利用することによって、良好な部分温度データを切り出すことのできる各種パラメータの調整に要する処理時間を短縮することができる。
【００３０】
このように、上記第１実施形態に係る流体の流れ計測装置１においては、時系列温度データを切り出し窓単位で切り出し、その切り出した窓を時間方向にずらしながら周波数解析を行うようにしているため、不規則に生じる温度変動を精度よく解析することができる。また、切り出し窓幅及び解析ステップ、また解析周波数を、調整可能に構成しているため、良好な解析結果、すなわち、温度変動が得られるような切り出し窓幅及び解析ステップ、また解析周波数を探索し、探索した値を用いて温度変動分布解析を行うため、より精度よく、温度変動分布解析を行うことができる。
そして、計測対象の気流が存在すると予測される領域を撮影し、解析するだけで温度変動分布解析を行うことができるため、簡易な構成で容易に気流の流れを計測することができる。
【００３１】
次に、本発明の第二実施形態を説明する。
第二実施形態に係る流体の流れ計測装置は、第一実施形態において、温度変動分布解析時の処理手順が異なる。流体の流れ計測装置の構成は、第一実施形態における流体の流れ計測装置の構成と同様である。
第一実施形態に係る流体の流れ計測装置１では、カメラ１１により温度画像データを取得して記憶装置１３に時系列に格納し、記憶装置１３に格納した時系列の温度画像データをもとに、温度変動分布解析を行っている。第二実施形態に係る流体の流れ計測装置１では、温度変動分布解析を、第一実施形態で説明した短時間フーリエ変換法と、リアルタイムで温度変動分布を取得するバンドパスフィルタ法とのいずれかの方法で計測できるようになっている。
【００３２】
第一実施形態で説明したように、短時間フーリエ変換法では、解析の時間幅を狭め、温度画像データを少しずつ時間方向に切り出しながら解析するようにしている。そのため、演算量が膨大となり処理時間を要する。第一実施形態では、例えば、ハードディスク（ＨＤＤ）等の記憶装置１３に、カメラ１１により取得した温度画像データを格納し、その後改めて温度画像データを読み出して解析している。つまり、解析結果が得られるまでに時間がかかり、温度変動に与える様々な要因を特定するのに時間を要する。
そこで、バンドパスフィルタ法では、フーリエ変換処理を行う代りに、デジタルフィルタとしてバンドパスフィルタ処理を行うことで、処理時間の短縮を図り、リアルタイムでの温度変動分布解析を可能にしている。バンドパスフィルタとしては、ＩＩＲ(Infinite Impulse Response)フィルタを用いる。
【００３３】
例えば、温度変動分布の計測を行う場合、短時間フーリエ変換法では、赤外線カメラで１秒間に１００枚程度のレートで計測する必要がある。つまり、温度画像データを１秒間に０．０１秒毎に取得する必要がある。例えば温度変動の周期が３秒程度であるとすると、短時間フーリエ変換法を用いて解析を行うためには、同じく３秒程度で切り出しながら解析する必要がある。つまり、３００枚程度の連続する画像を部分温度データとして解析する必要がある。例えば、６４０×５１２画素の赤外線カメラを用いた場合、演算時間に約０．７秒を要することになり、実時間での処理は困難である。また、一旦温度画像データを記憶装置に記憶した後に、記憶した温度画像データを用いて演算するにしても、６０秒間の温度画像データとした場合、６０×１００×０．７＝４２００秒、つまり、約７０分かかることになる。
【００３４】
一方、バンドパスフィルタとしてＩＩＲフィルタの一つである双二次フィルタを用いた場合には、５枚の温度画像データを用いれば、簡易的ではあるが温度変動分布画像を得ることができる。つまり、演算時間は、単純には１／６０以下となる（実際には０．００３秒）。
したがって、温度画像データの取得間隔を０．０１秒とすると、温度画像データの取得間隔のタイミングで温度画像データから温度変動分布画像を十分得ることができる。そのため、温度画像データとともに温度変動分布画像を同時にリアルタイムで表示することが可能になる。
【００３５】
バンドパスフィルタとして双二次フィルタを用いた場合、フィルタ係数は次式（１）に示すように表すことができる。
ａ０＝１＋α
ａ１＝−２×ｃｏｓ（ω）
ａ２＝１−α
ｂ０＝α
ｂ１＝０
ｂ２＝−α
なお、ω＝２．０×π×中心周波数／サンプリングレート
α
＝sin(ω)×sinh{(log(２)／２)×帯域幅×ω／sin(ω)}          ……（１）
【００３６】
上記式（１）に示すフィルタ係数を用いたバンドパスフィルタ処理は、次式（２）に示す通りとなる。つまり、温度変動分布を表す画像を得るためには、５枚の温度画像で実現することができる。
出力画像＝（ｂ０／ａ０）×入力画像
＋（ｂ１／ａ０）×１つ前の入力画像
＋（ｂ２／ａ０）×２つ前の入力画像
−（ａ１／ａ０）×１つ前の出力画像
−（ａ２／ａ０）×２つ前の出力画像
……（２）
【００３７】
なお、デジタルフィルタには、大きく分けてＩＩＲフィルタとＦＩＲ（Finite Impulse Response）フィルタとがあるが、ＩＩＲフィルタの方がより少ないタップ（画像数）でフィルタを構成することができるため、演算時間も短くてすみ、実時間での処理に有利である。
このように、バンドパスフィルタ法を用いることで、簡易的な温度変動分布画像を得ることができる。
そこで、第二実施形態に係る流体の流れ計測装置１では、温度変動分布解析を、短時間フーリエ変換法と、バンドパスフィルタ法とのいずれの方法で行うかを選択可能に構成し、簡易的ではあるがリアルタイムでの温度変動分布解析と、リアルタイムでの演算処理は困難ではあるがより詳細な温度変動分布画像を得ることのできる短時間フーリエ変換法を用いた温度変動分布解析とを、切り替え可能としている。
【００３８】
図６は、第二実施形態に係る流体の流れ計測装置１における処理手順の一例を示すフローチャートである。
流体の流れ計測装置１では、画像処理装置１２の入力装置により、周波数解析法として短時間フーリエ変換法とバンドパスフィルタ法とのいずれかが設定されると（ステップＳ２１  特許請求の範囲に記載の選択部に対応）、短時間フーリエ変換法が設定されたときには、ステップＳ２２に移行し、図３に示す第一実施形態における、撮影時の処理と同様の処理を行い、撮影した温度画像データを記憶装置１３に格納した後、記憶した温度画像データをもとに、第一実施形態における短時間フーリエ変換法による温度変動分布解析を行う。
【００３９】
一方、周波数解析法としてバンドパスフィルタ法が設定されたときには、ステップＳ３１に移行する。そして、画像処理装置１２の入力装置により、計測対象の流体を特定する情報が入力されると、入力された計測対象の流体に対応する中心周波数を設定する（ステップＳ３２）。この中心周波数はバンドパスフィルタの中心周波数である。中心周波数は、計測対象の気流の温度変化の周期に対応させて設定している。
【００４０】
ステップＳ３２でバンドパスフィルタの中心周波数を設定したならばステップＳ３３に移行し、流体の流れ計測装置１では、予め設定されたフレームレートでカメラ１１により計測対象の流体を撮影する。画像処理装置１２は、カメラ１１で撮影した温度画像データを読み込み、ステップＳ３２で設定された中心周波数を有し、中心周波数を中心として予め設定した帯域幅の通過帯域を有する、ＩＩＲフィルタを用いてフィルタ処理を行う（ステップＳ３４  特許請求の範囲に記載のフィルタ処理部に対応）。そして、ステップＳ３５に移行し、フィルタ処理後の温度画像データを、温度変動分布画像として、画像処理装置１２の表示装置に表示する（特許請求の範囲に記載の表示処理部に対応）。なお、画像処理装置１２は、カメラ１１からの温度画像データを読み込むと共に、この温度画像データを記憶装置１３に時系列に格納するようにしてもよい。
【００４１】
なお、バンドパスフィルタの中心周波数や、帯域幅は、計測対象の流体に対応する温度変動の周波数に応じて設定すればよい。また、サンプリングレートは、温度変動の周波数、カメラ１１の仕様や画像処理装置１２の処理性能に応じて設定すればよい。
また、中心周波数は複数段階に設定可能に構成してもよく、また、帯域幅等といった各種パラメータは固定であってもよく、複数段階に設定可能に構成してもよい。また、中心周波数を、自動又は手動によってスイープさせるようにしてもよい。このようにスイープさせることで、中心周波数等といった各種パラメータを、その場で調整し最適な値に設定することができる。
【００４２】
このように、上記第二実施形態では、周波数解析法として、短時間フーリエ変換法とバンドパスフィルタ法とをユーザが選択できるようになっている。そのため、例えば、リアルタイムで流体の流れを確認したいときには、バンドパスフィルタ法を選択することによって、温度画像データを取得しつつその場で、温度変動分布画像を表示することができ、すなわち流体の流れを可視化することができる。逆に、リアルタイムでなくともよく、流体の流れを詳細に可視化したい場合には、短時間フーリエ変換法を用いることにより、流体の流れをより詳細に可視化することができる。
【００４３】
次に、本発明の第三実施形態を説明する。
第三実施形態に係る流体の流れ計測装置１ａは、第二実施形態に係る流体の流れ計測装置１において、周波数解析法としてバンドパスフィルタ法のみを備えたものである。この場合には、カメラ１１から入力される温度画像データを記憶しなくともよいため、図７に示すように、流体の流れ計測装置１ａは、温度画像データを記憶するための記憶装置を備える必要はなく、カメラ１１と、カメラ１１からの温度画像データの取得及びフィルタ処理を行う画像処理装置１２と、を備えていればよい。
【００４４】
なお、カメラ１１と温度画像データの取得及びフィルタ処理を行う画像処理装置１２とを、一体型としてもよい。例えば、図８に示すように、表示装置１１０ａを備えた携帯可能な赤外線カメラ１１０からなるカメラ１１に画像処理装置１２の機能を組み込んでもよい。この場合、例えば、画像処理装置１２の機能をＦＰＧＡ（field-programmable gate array）等の書き換え可能な回路に組み込み、この回路を、赤外線カメラ１１０に搭載すればよい。
また、画像処理装置１２として、携帯型のパーソナルコンピュータ等、簡易なパーソナルコンピュータを適用した場合であっても、リアルタイムで温度変動分布画像を表示することができ、使い勝手を向上させることができる。
【００４５】
なお、上記第二及び第三実施形態において、バンドパスフィルタ法では、リアルタイムで温度変動分布画像を取得するようにしているが、所定時間分の連続した温度画像データを一旦記憶装置に格納し、記憶装置に格納した連続した温度画像データに対してバンドパスフィルタ処理を行うことによって、簡易的な温度変動分布画像を取得するようにしてもよい。
また、第一実施形態及び第二実施形態に係る流体の流れ計測装置において、カメラ１１と画像処理装置１２と記憶装置１３とを一体に形成してもよい。
【実施例】
【００４６】
次に、本発明の実施例について説明する。
［実施例１］
本発明の第一実施形態に係る流体の流れ計測装置１を用いて、気流の流れを計測した。ここでは、図９に示すように、室内の天井付近の気流を計測対象とした。具体的には、エアーコンディショナーによって生じる気流の流れを計測した。なお、図９は、計測対象の気流が存在する領域の可視画像である。
カメラ１１として、中赤外帯域の冷却型（素子ＩｎＳｂ）の赤外線カメラを用いた。カメラ１１の温度画像の画素数は、６４０×５１２画素である。
【００４７】
撮影フレームレートは９７Ｈｚとした。露光時間は２０００μｓｅｃ．、撮影時間は６０ｓｅｃとした。
そして、計測対象の気流が存在する領域の温度画像データを記憶装置１３に時系列に記憶し、次に、画像処理装置１２において、短時間フーリエ変換により、温度変動分布解析を行った。
短時間フーリエ変換の切り出しの窓幅は２．４ｓｅｃ．とし、解析の周波数は０．３７Ｈｚとし、窓関数は矩形窓とし、窓関数で切り出す切り出し窓の解析ステップは０．２ｓｅｃ．とした。
【００４８】
図１０に、計測対象の気流が存在する領域の温度画像を示す。図１０に示すように、ある時点の温度画像からは、温度変動の分布を認識することはできない。
図１１に、計測対象の気流が存在する領域の短時間フーリエ変換による温度変動分布解析の結果を示す。温度変動が比較的大きい部分が時間の経過と共に、天井に上昇し、天井に沿って移動していることがわかる。つまり、時間の経過に伴う温度変動を表す濃淡の分布の変化から、気流の流れが良好に可視化できることがわかる。
【００４９】
本実施形態においては、連続する温度画像を短時間の切り出し窓で切り出して解析するため、時間分解能は高くなる反面、周波数分解能は低くなる。したがって、解析結果に及ぼす周波数の影響は小さい。そのため、解析の周波数の選定にそれほど注意する必要はない。解析の周期（周波数の逆数）が、切り出しの窓幅と同程度か、それよりもやや大きいぐらいが適当である。周波数特性に広がりを持つため、解析の周波数を厳密に設定しなくとも、一つの周波数で解析を行うことができる。
すなわち、窓幅（時間分解能）と周波数分解能とは相反する関係にあり、一方を向上させようとすると、他方が低下する。つまり、時間分解能を向上させると（窓幅を狭めると）周波数分解能が低下し、時間分解能を低下させると（窓幅を広げると）、周波数分解能が向上する。
【００５０】
この様子を図１２に示す。図１２は、サンプリング周波数を３１ｍｓｅｃとし、１Ｈｚの信号に対して切り出しの窓の幅を変化させて周波数解析した例である。図１２において、横軸は周波数Ｈｚ、縦軸は振幅である。図１２において、特性線Ｌ１は窓幅を４０９６点とした場合、特性線Ｌ２は窓幅を５１２点とした場合、特性線Ｌ３は窓幅を６４点とした場合を示す。なお、窓幅が例えば６４点とは、３１ｍｓｅｃのサンプリング周期で、６４枚の連続する温度画像を切り出して、周波数解析の対象とすることを意味する。
【００５１】
特性線Ｌ１に示すように窓幅が広い４０９６点では１Ｈｚ近傍の周波数の信号にしか感度を持たないが、特性線Ｌ３に示すように窓幅が狭い６４点では０.５〜１.５Ｈｚの範囲に感度を持ち、窓幅が５１２点では特性線Ｌ２に示すように窓幅が広い４０９６点での感度よりも周波数の範囲が広い１Ｈｚ前後に感度を持つことがわかる。したがって、窓幅が広い４０９６点では、少しの周波数の誤差であっても解析が行えないが、窓幅が狭い６４点ではある程度の周波数の誤差も許容される。
【００５２】
［実施例２］
次に、本発明の第三実施形態に係る流体の流れ計測装置１ａを用いて、気流の流れを計測した。
ここでは、エアダスター噴射による気流の流れを計測した。
カメラ１１として、中赤外帯域の冷却型（素子ＩｎＳｂ）の赤外線カメラを用いた。カメラ１１の温度画像の画素数は、６４０×５１２画素である。
撮影フレームレートは９７Ｈｚとし、露光時間は２０００μｓｅｃ．撮影時間は６０ｓｅｃ．とした。
バンドパスフィルタとしては、ＩＩＲフィルタの一種である双二次フィルタを用いた。フィルタ係数は（１）式に示す通りであり、バンドパスフィルタのフィルタ処理は（２）式に示す通りである。なお、中心周波数（カットオフ周波数）は０．３Ｈｚ、サンプリングレートは９７Ｈｚ、帯域幅は中心周波数を中心として２Ｈｚとした。
【００５３】
図１３に、エアダスターが噴射する気流が存在する領域の温度画像を示す。図１３に示すように、温度画像からは、気流を視認することはできない。
図１４に、バンドパスフィルタ法により解析処理を行うことにより得られる温度変動分布画像の一例を示す。図１４は、０．１秒毎の解析結果である。
図１４に示すように、エアダスターを噴射させたときに生じる温度差が、気流の流れとして可視化できることが確認された。
なお、中心周波数が０．３Ｈｚであり、帯域幅が２Ｈｚとした場合、バンドパスフィルタの通過帯域ｆは、０Ｈｚ≦ｆ≦１．３となる。すなわちローパスフィルタと同等となるため、バンドパスフィルタに代えてローパスフィルタを用いてもよい。
【００５４】
また、上記各実施形態においては、エアーコンディショナーによる気流の流れや、エアダスター噴射による気流の流れを計測した場合について説明したがこれに限るものではない。
上記各実施形態に係る流体の流れ計測装置を、例えば、石油化学プラントや、パイプライン等からのガスの漏出を検知するガス検知装置として適用することも可能である。
ここで、石油化学プラントやパイプライン等から漏出するガスの中には、人体に有害であるガスや、災害や環境汚染の引き金にもなり得るガス等を含む可能性があることから、迅速に検知する必要がある。
【００５５】
石油化学プラントやパイプライン等からのガスの漏出を検知する場合、単純に漏出したガスと周囲との温度差によってガスの漏出を検出することも可能であるが、本発明の各実施形態に係る流体の流れ計測装置を用いて温度変動分布解析を行い、ガスの流れを計測することによって、ガス漏出の検出精度をより向上させることができる。そのため、ガスの漏出の見逃しやガスの微量な漏出も検出することができ、また、ガスの漏出源を容易に特定することができる。
【００５６】
なお、漏出ガスの検知に限るものではないが、上記各実施形態に係る気体の流れ計測装置において、特定のガスの検知を対象とする場合には、そのガスの吸収帯の波長のみを取り込むようにしてもよい。例えば炭化水素ガスを対象とする場合には、吸収帯である３．２μｍ以上３．４μｍ以下の赤外線のみをカメラ１１に取り込むことを目的に、赤外線バンドパスフィルタを通過させてもよい。各実施形態に係る流体の流れ計測装置において、赤外線バンドパスフィルタを通過させた後のデータに対して温度変動分布解析を行うことによって、より微量な漏出ガスや、温度変化の小さな漏出ガスの検知を行うことができる。
【００５７】
つまり、赤外線バンドパスフィルタを設けたとき、特定のガスとして炭化水素ガスを適用した場合、炭化水素ガスからの赤外線放射と背景からの赤外線放射の両者を合成した成分が測定されるが、炭化水素ガスが存在すれば、炭化水素ガスの吸収帯（３．２μｍ以上３．４μｍ以下）の波長の背景からの赤外線放射は吸収されるため小さくなり、炭化水素ガスが存在しなければ、背景からの赤外線放射は吸収されないため小さくならない。そして、赤外線バンドパスフィルタを通過した炭化水素ガスの吸収帯の波長に狭められた赤外線放射に対して、温度変動分布解析を行うことになるため、背景からの放射の影響を軽減することができ、より精度よく炭化水素ガスの温度変動分布を求めることができ、結果的に漏出ガスの検知精度を向上させることができる。あるいは、逆に、吸収が大きいということは放射率が高いということなので、炭化水素ガス自体による温度変動分布が存在する場合には、漏出ガスの検知精度が向上するという効果もある。
【００５８】
一方で、赤外線バンドパスフィルタを通過させない場合にもあらゆる種類のガスの漏出の検知を行うことができるため、例えば図１５に示すように、カメラ１１に赤外線バンドパスフィルタの有無、或いは通過帯域を選択することの可能な、通過帯域の異なる複数のフィルタ１２０ａを有するフィルタホイール１２０を設け、検知対象のガスに応じた通過帯域を有する赤外線バンドフィルタを選択するようにしてもよい。
赤外線バンドパスフィルタは、カメラ１１のレンズの前に設けてもよいし、レンズと素子の間に設けてもよい。
【００５９】
また、ガス漏出の検知以外にも、車両の排気管から排出される排気ガスの流れ、ドライヤから送出される熱風の流れ、配管からのエアの漏れ、建物内の空調の流れ、乾燥工程でのエアの流れ等を計測することも可能である。また、気体の流れに限るものではなく、液体の流れ等を計測することも可能である。例えば、海洋・河川などの液体の流れ、射出成型時の樹脂の流れ、溶湯の流れ、塗装時の溶剤の流れ、などの可視化にも適用可能である。
また、可視化した温度変動分布画像から、流れの速度ベクトルを求めることも可能である。
具体的には、まず、時刻ｔにおける温度変動分布画像をＩ（ｘ，ｙ，ｔ）とする。ここで、ｘは横方向の画素番号（ｘ＝１，２，…，Ｎ）、ｙは縦方向の画素番号（ｙ＝１，２，…，Ｍ）、Ｎは温度変動分布画像の横方向の画素数、Ｍは縦方向の画素数である。
【００６０】
ここで、図１６に示すように、時刻ｔ１における温度変動分布画像上の点（ｘ１，ｙ１）が、時刻ｔ２において対応する点（ｘ２，ｙ２）に移動していたとすると、点（ｘ１，ｙ１）の移動速度Ｖ（ｘ１，ｙ１，ｔ１）は、次式（３）で表すことができる。
Ｖ(ｘ１，ｙ１，ｔ１)
＝Ｋ・{(ｘ２−ｘ１)２＋(ｙ２−ｙ１)２}１／２／(ｔ２−ｔ１)    ……（３）
なお、（３）式中のＫは、温度変動分布画像上の１画素が実際の寸法ではいくつかを表す画素分解能である。
【００６１】
また、点（ｘ１，ｙ１）が点（ｘ２，ｙ２）に移動したときの速度ベクトルの方位角φ（ｔ１）は、次式（４）で表すことができる。
φ(ｔ１)＝tan−１{(ｙ２−ｙ１}／(ｘ２−ｘ１)}      ……（４）
ただし、φは、右横向きを基準とし、反時計回りを正に定義している。また、ｙ２−ｙ１＞０のときは０＜φ＜１８０°、ｙ２−ｙ１＜０のときは１８０＜φ＜３６０°となるように定義する。
【００６２】
ここで、時刻ｔ１における温度変動分布画像上の点に対応する時刻ｔ２における温度変動分布画像上の点を抽出する手段について説明する。時刻ｔ1における温度変動分布画像上の点（ｘ1，ｙ1）と対応する時刻ｔ２における温度変動分布画像上の点（ｘ２，ｙ２）は、次式（５）で表される残差Ｓ（ｘ，ｙ）を最小にする（ｘ，ｙ）として求めることができる。なお、（５）式中のｉｓ，ｉｅ，ｊｓ，ｊｅは残差を計算する範囲を表し、流れの特徴が表されるような範囲に設定する。つまり、時刻ｔ1における温度変動分布は、時刻ｔ２においても同じように温度変動分布が生じるとみなし、時刻ｔ1における温度変動分布画像において流れの特徴が表されるような温度変動分布が生じている領域を特徴領域とし、時刻ｔ２における温度変動分布画像において、温度変動の分布が特徴領域と略相似な形状を有する領域を、時刻ｔ２における、特徴領域に対応する対応領域とみなし、時刻ｔ1における特徴領域と、時刻ｔ２における対応領域との移動状況から、特徴領域の速度ベクトルを演算する。
【００６３】
【数１】
2020052036.tif 000004
【００６４】
ここで、「流れの特徴が表されるような範囲」というのは、流れに対応する温度変化の空間的な分布が略相似な形状を有する領域を１対１で見つけることができるのに適当な大きさの範囲ということになる。範囲が小さすぎると類似の箇所が多く見つかってしまい、どれが本当に対応する箇所かがわからなくなってしまうし、逆に大きすぎると温度変動の分布が変わってしまって対応点が見つからなくなってしまう心配がある。
【００６５】
時刻ｔ1における温度変動分布画像上の点（ｘ1，ｙ1）と対応する時刻ｔ２における温度変動分布画像上の点（ｘ２，ｙ２）は、残差の代わりに、次式（６）で表される相関係数を最大にする点として求めることもできる。
【００６６】
【数２】
2020052036.tif 000005
【００６７】
ここで、（６）式中のＩave（ｘ，ｙ，ｔ）は、時刻ｔの（ｘ，ｙ）におけるＩの平均であり、次式（７）で表される。
【００６８】
【数３】
2020052036.tif 000006
【００６９】
また、（６）式中のσ（ｘ，ｙ，ｔ）は時刻ｔ、座標（ｘ，ｙ）におけるＩ（ｘ，ｙ，ｔ）の標準偏差であり、次式（８）で表される。
【００７０】
【数４】
2020052036.tif 000007
【００７１】
図１４に示す温度変動分布画像に対して、速度ベクトルを求めた結果の一例を、図１７に示す。
図１７において、温度変動分布画像上に重畳した矢印の大きさで速度の大きさを表し、矢印の向きで方向を表している。また、温度変動分布画像上の全ての点（画素）で対応付けを行うと誤りも多くなることから、標準偏差σ(ｘ，ｙ，ｔ1）が所定の値以上の特徴的な点のみに注目し、対応点を探すようにしている。
【００７２】
ここで、速度の大きさを絶対値で計算する際に使用した画素分解能Ｋは、流れの奥行き方向の位置が必要である。また、画像内に複数の奥行き方向の位置の流れが混在する場合などには、可視化された気流の流れが、奥行方向のどの位置に存在するかを知りたい場合もある。
その場合は、図１８に示すように、２台の赤外線カメラ２１Ｌ及び２１Ｒを設ける。このとき、赤外線カメラ２１Ｌ及び２１Ｒの光軸が平行となり、且つ互いの視野が重複するように配置する。
【００７３】
そして、赤外線カメラ２１Ｌ及び２１Ｒで撮影した温度画像データそれぞれに対し、上記各実施形態で説明した手順で温度変動分布画像を取得する。つまり、気流の流れを計測する。そして、各赤外線カメラ２１Ｌ及び２１Ｒの温度画像データに基づく温度変動分布画像をもとに、奥行情報算出部において、気流の流れの奥行方向の位置を算出する。例えば、この奥行情報算出部は画像処理装置１２に設け、画像処理装置１２において、上記各実施形態で説明した手順で温度変動分布画像を取得すると共に、奥行の演算も行う。
具体的には、各赤外線カメラ２１Ｌ及び２１Ｒからの温度画像データそれぞれに基づく温度変動分布の画像から、対応する同一の流れ部分を抽出し、通常の可視カメラではよく実施されるステレオ立体視の原理で、同一の流れ部分の奥行位置を演算する。
【００７４】
まず、同一の流れ部分を抽出する方法について説明する。
同一の流れ部分を抽出する方法は、上述した異なる時刻の温度変動分布画像間の対応を取る方法を、同時刻の赤外線カメラ２１Ｌ及び２１Ｒの温度画像データに基づく温度変動分布画像に適用したものである。すなわち、赤外線カメラ２１Ｌの温度画像データに基づく温度変動分布画像をＩＬ（ｘ，ｙ，ｔ）、赤外線カメラ２１Ｒの温度画像データに基づく温度変動分布画像をＩＲ（ｘ，ｙ，ｔ）（ｘ＝１，２，…，Ｎ，ｙ＝１，２，…，Ｍ）とする。
このとき、時刻ｔ０における赤外線カメラ２１Ｌに基づく温度変動分布画像上の点（ｘ０，ｙ０）と対応する赤外線カメラ２１Ｒに基づく温度変動分布画像上の点は、次式（９）で表される残差Ｓ（ｘ，ｙ）を最小にする点（ｘ，ｙ）として求めることができる。ｉｓ，ｉｅ，ｊｓ，ｊｅは残差を計算する範囲を表し、流れの特徴が表されるような範囲に設定する。
【００７５】
【数５】
2020052036.tif 000008
【００７６】
あるいは、残差の代わりに、次式（１０）で表される相関係数を最大にする点（ｘ，ｙ）として求めることもできる。
【００７７】
【数６】
2020052036.tif 000009
【００７８】
ここで、（１０）式中のＩＬave（ｘ０，ｙ０，ｔ０）は、時刻ｔ０，（ｘ０，ｙ０）におけるＩＬの平均、ＩＲave（ｘ，ｙ，ｔ０）は、時刻ｔ０，（ｘ，ｙ）におけるＩＲの平均であって、次式（１１）で表される。
【００７９】
【数７】
2020052036.tif 000010
【００８０】
また、（１０）式中のσＬ（ｘ０，ｙ０，ｔ０）は時刻ｔ０，（ｘ０，ｙ０）におけるＩＬの標準偏差、σＲ（ｘ，ｙ，ｔ０）は時刻ｔ０，（ｘ，ｙ）におけるＩＲの標準偏差であって、次式（１２）で表される。
【００８１】
【数８】
2020052036.tif 000011
【００８２】
奥行情報は、図１８に示すように、２台の赤外線カメラ２１Ｌ、２１Ｒが、間隔２Ｄを空けて平行に並んでいるとし、赤外線カメラ２１Ｌ、２１Ｒを含む平面内で考える。
ここで、赤外線カメラ２１Ｌ、２１Ｒを起点として奥行距離Ｚの位置にある点Ｐは、赤外線カメラ２１Ｌの撮像素子２２上の座標ｔＬの点で撮像され、同様に、赤外線カメラ２１Ｒの撮像素子２２上の座標ｔＲの点で撮像されるとする。なお、座標ｔＬおよび座標ｔＲは、いずれも光軸中心を原点とし、右向きを正とする。赤外線カメラ２１Ｌ、２１Ｒのレンズ２３の焦点距離をｆとし、Ｚが十分遠いとすると、座標ｔＬ及びｔＲは、近似的に次式（１３）及び（１４）で表すことができる。
【００８３】
−ｔＬ／ｆ＝ＸＬ／Ｚ          ……（１３）
−ｔＲ／ｆ＝ＸＲ／Ｚ          ……（１４）
ここで、ＸＬ−ＸＲは、次式（１５）で表すことができる。
ＸＬ−ＸＲ＝２Ｄ            ……（１５）
（１３）式〜（１５）式から次式（１６）が得られる。したがって、奥行距離Ｚは（１７）式から得ることができる。
（ｔＲ−ｔＬ）／ｆ＝２Ｄ／Ｚ          ……（１６）
Ｚ＝２Ｄｆ／（ｔＲ−ｔＬ）            ……（１７）
【００８４】
また、図１８では、２台の赤外線カメラ２１Ｌ、２１Ｒが平行に配置されている場合を考えたが、角度をつけて配置されている場合、つまり、赤外線カメラ２１Ｌ、２１Ｒの光軸が交差するように配置されている場合も、奥行を算出することができる。
図１９に、赤外線カメラ２１Ｌが角度θＬ０、赤外線カメラ２１Ｒが角度θＲ０だけ傾けられ、光軸が交差するように配置されている場合の図を示す。ただし、角度は、赤外線カメラ２台を平行に並べた時の光軸を基準とし、時計回りを正と定義する。
このとき、次式（１３ａ）、（１４ａ）、（１８）が成り立つ。
【００８５】
−ｔＬ／ｆ＝ｔａｎθＬ          ……（１３ａ）
−ｔＲ／ｆ＝ｔａｎθＲ          ……（１４ａ）
ＸＬ／Ｚ＝ｔａｎ（−θＬ０＋θＬ）……（１８）
したがって、（１３ａ）、（１４ａ）、（１８）式から次式（１３ｂ）、（１４ｂ）が成り立つ。
ＸＬ／Ｚ＝-tan(θＬ０＋tan−１(ｔＬ／ｆ))      ……（１３ｂ）
ＸＲ／Ｚ＝-tan(θＲ０＋tan−１(ｔＲ／ｆ))      ……（１４ｂ）
【００８６】
これら（１３ｂ）式及び（１４ｂ）式を、（５）式に代入すると、次式（７ａ）が得られる。
Ｚ・{tan(θＲ０＋tan−１(ｔＲ／ｆ))−tan(θＬ０＋tan−１(ｔＬ／ｆ))｝
＝２Ｄ
Ｚ
＝２Ｄ／{tan(θＲ０＋tan−１(ｔＲ／ｆ))−tan(θＬ０＋tan−１(ｔＬ／ｆ))}
……（７ａ）
θＲ０＝θＲ０＝０の場合は、（７）式と（７ａ）式とは一致する。
【００８７】
図２０は、奥行き方向の異なる位置に置いた２つのエアダスターからの噴射を可視化した例である。上述のステレオ立体視の方法で奥行き方向の位置を算出し、さらに、対応付けを行った点の時間方向の変化から速度ベクトルを算出した。このとき速度算出に用いている画素分解能は奥行き方向の位置を考慮して決定している。その結果、速度ベクトルを示す矢印の色、矢印を表す線種等といった表示形態を奥行き方向の位置により変えて表示するようにしている。
左右の赤外線カメラ２１Ｌ、２１Ｒそれぞれで取得した温度変動分布画像間で、同一位置を表す画素同士の対応を多くの点で取ることができる場合には、速度ベクトルを示す矢印は表示せずに、画像そのものに色等を付けて奥行き情報を表現することもできる。また、表示する奥行き方向の範囲を指定し、その範囲に含まれる流れのみ表示することも可能である。
【００８８】
なお、本発明の範囲は、図示され記載された例示的な実施形態に限定されるものではなく、本発明が目的とするものと均等な効果をもたらす全ての実施形態をも含む。さらに、本発明の範囲は、全ての開示されたそれぞれの特徴のうち特定の特徴のあらゆる所望する組み合わせによって画され得る。
【符号の説明】
【００８９】
１、１ａ  流体の流れ計測装置
１１  カメラ
１２  画像処理装置
１３  記憶装置
２１Ｌ、２１Ｒ  赤外線カメラ
２２  撮像素子
１１０  赤外線カメラ
１１０ａ  表示装置
１２０  フィルタホイール
【図１】
2020052036.tif 000012
【図２】
2020052036.tif 000013
【図３】
2020052036.tif 000014
【図４】
2020052036.tif 000015
【図５】
2020052036.tif 000016
【図６】
2020052036.tif 000017
【図７】
2020052036.tif 000018
【図８】
2020052036.tif 000019
【図９】
2020052036.tif 000020
【図１０】
2020052036.tif 000021
【図１１】
2020052036.tif 000022
【図１２】
2020052036.tif 000023
【図１３】
2020052036.tif 000024
【図１４】
2020052036.tif 000025
【図１５】
2020052036.tif 000026
【図１６】
2020052036.tif 000027
【図１７】
2020052036.tif 000028
【図１８】
2020052036.tif 000029
【図１９】
2020052036.tif 000030
【図２０】
2020052036.tif 000031