JP3554822B2 - 最適化排他的二値化相関計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラム - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は，最適化排他的二値化相関計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラムに関するものである。特に，移動する計測対象の信号である計測信号と参照信号の相関を計測するものである。高速に移動する物体の移動速度，移動方向等を観測することを可能にするものであり，地球大気のゆらぎの観測等にも利用できるものである。
【０００２】
星のまたたきは，星からの光が地球大気圏の屈折率ゆらぎによって波面（位相）を乱されることによって生じる。この星の瞬きが地表に生成するランダムパターンを観測することによりゆらぎ層における風向，風速，ゆらぎ層の高度，ゆらぎの強度等を計測することができる。このランダムパターンは高速に移動するので，高速な画像処理を必要とする。
【０００３】
【従来の技術】
このような，ランダムパターンの観測には，ランダムパターンの二次元画像の一部の画像信号を参照信号として高速に取得し，ランダムパターンの二次元画像を計測信号として，両者の相関をとることによりランダムパターンを計測する。
【０００４】
図１２（ａ）は計測信号を表し，図１２（ｂ）は参照信号を表す。
Ａは計測対象の計測信号，Ｂは参照信号である。
【０００５】
計測信号Ａと参照信号Ｂの相関をとり，両者が同一の信号であるかどうかを判別する。計測信号Ａと参照信号Ｂは，Ａ＝＜Ａ＞＋α，Ｂ＝＜Ｂ＞＋βで表すことができる。
【０００６】
【数１】

【０００７】
【数２】

【０００８】
であり，＜Ａ＞，＜Ｂ＞は，それぞれの信号強度の平均値を表す。αは計測信号Ａの平均値からの変動量，βは参照信号Ｂの平均値からの変動を表す。
【０００９】
計測信号Ａと参照信号Ｂとの相関Ｃは次の式で計算する。
【００１０】

ここに，＜＞はＡとＢの平均を表す。ＡおよびＢがランダムに変動するものであると，＜α＞＝０，＜β＞＝０である。
【００１１】
また，計測信号Ａと参照信号Ｂが同一でないとき時（Ａ≠Ｂ）
＜αβ＞＝０
計測信号Ａと参照信号Ｂが同一である時（Ａ＝Ｂの時）
＜αβ＞＝σ^２ となる。σ^２ は計測信号Ａと参照信号の分散である。
【００１２】
従って，Ｃ＝＜Ａ＞^２ （Ａ≠Ｂの時）
Ｃ＝＜Ａ＞^２ ＋σ^２ （Ａ＝Ｂの時）
で表される。
【００１３】
例えば，図１（ａ）の装置でランダムに変動するパターンが，パターン形状を変えないで移動する場合の移動速度を計測すると考える（図１（ａ）の詳細は後述する）。
【００１４】
アレイ素子光検出器５と単素子光検出器６は，同期して複数（ｉ＝１，２，・・・・，Ｎ）のデータ対を取得する。
【００１５】
計測信号Ａはランダムパターン１がアレイ素子光検出器５の任意の一つの画素で観測される信号Ａ_ｉ の集合である。参照信号Ｂはピンホール７を介して単素子光検出器６で観測される信号で，対となっているアレイ素子光検出器５のデータ取得時刻から遅延時間τだけ遅れて取得された信号Ｂ_ｉ の集合である。
【００１６】
図１（ａ）の装置で，参照信号と計測信号の一般的な相関Ｃは上記の式のようになる。
【００１７】
その結果からランダムパターン１が形状を変えずに移動している場合には，この相関Ｃを計算することにより，ランダムパターン１の移動する方向と速度を
求めることができる。
【００１８】
この一般的な相関Ｃの計算に対して，相関を計算する前に計測信号Ａのｉ番目の要素と参照信号Ｂのｉ番目の要素を比較して，それらが異なる場合は両方の要素の値をゼロにしてから排他的相関Ｃ’を
Ｃ’＝＜Ａ×Ｂ＞
を計算することにより
Ｃ’＝０ （Ａ≠Ｂ）
Ｃ’＝＜Ａ＞^２ ＋σ^２ （Ａ＝Ｂの時） （２０） のように相関の無い部分の値を減らすことで，ピークをより鮮明にすることができる。その結果として少ない信号要素数で相関を計算することができる。
【００１９】
【発明が解決しようとする課題】
上記の方法は，後述するように，計測信号か参照信号，あるいはその両方に計測雑音が含まれる場合には，Ａ＝Ｂであった信号が，Ａ＝Ｂ＋ε（εは計測雑音）と変化してしまうためそのまま適用すると上記式（２０）で示されたピークの高さが大幅に減少してしまう。そのため，Ａ＝Ｂ±δの範囲にある場合は，Ａ＝Ｂと見なす許容誤差幅δを導入する必要がある。この許容誤差幅δの値をどのように設定するかは，以下に説明するように非常に重要である。
【００２０】
このδのを適切に設定しないと，例えば，δが小さすぎる場合は，本来はＡ＝Ｂであった信号がＡ≠Ｂと判定される確率が高くなり，相関のピークが減少する。また，δが大きすぎる場合は，本来はＡ≠Ｂであった信号がＡ＝Ｂと判定される確率が高くなり，相関の無い部分の値が増加し，ピークを明瞭に判定できなくなる。そこで，本発明では信号変動幅と計測雑音を考慮して最適な許容誤差幅δを求め，排他的二値化を最適化して相関計測を行なう最適化排他的二値化相関計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラムを提供することを目的とする。
【００２１】
【課題を解決するための手段】
本発明は，参照信号と移動する計測対象の信号である計測信号の相関を計測する最適化排他的二値化計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラムにおいて，参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をし，参照信号と計測信号の相関を計測するようにした。
【００２２】
そして，許容誤差幅は，計測環境に応じて定めるものである。許容誤差幅は計測環境中の計測雑音に応じて定まるものであり，参照信号と計測信号の差と最適化された該許容誤差幅の大小関係に応じて排他的二値化を行なうようにした。
【００２３】
さらに，相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように該許容誤差幅を最適化するものである。
【００２４】
ここでＳＮ比とは，相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量で，ＳＮ比が大きいほど相関ピークが鮮明であることを示すものである。
【００２５】
本発明は，相関をとる二次元画像の計測信号と参照信号の信号強度が同じであるかないかを判断を，ＳＮ比ができるかぎり良くなるように求めた許容誤差幅に基づいて排他的に二値化して計測信号と参照信号の相関をとるようにしたので，ＳＮ比が良い相関が得られる。また，信号強度の異なる部分は画素値を，例えば，０にし，信号強度が同じ部分の画素値を，例えば，１にするように二値化してから，相関をとるようにしたので，保存データ量も少なくてすみ，簡単な計算で高速に相関処理を行なうことができる。
【００２６】
【発明の実施の形態】
図１は本発明の実施の形態１を示す図である。図１（ａ）は観測装置を示し，図１（ｂ）は観測方法の説明図である。
【００２７】
図１（ａ）において，１は計測対象の光パターンであり，ランダムパターンである。ランダムパターン１は，例えば，高速で移動する物体から放射される光，あるいは星が地表面につくるランダムパターン等である。３は集光レンズであり，ランダムパターン１から放射される光を集光するものである。集光レンズ３は望遠鏡でも良い。４はビームスプリッタであり，ランダムパターン１から放射される光を計測光８と参照光９に分離するものである。５はアレイ素子光検出器であり，例えば，ＣＣＤのような二次元の光検出器である。６は単素子光検出器であり，ピンホール７を通過する狭い視野の光を検出するものである。単素子光検出器６はアレイ素子光検出器５に比べて高速にサンプリングできるものである（例えば，光電子増倍管等を使用する）。７はピンホールである。８は計測光であり，アレイ素子光検出器５に入射する光である。９は参照光であって，単素子光検出器６に入射する光である。
【００２８】
図１（ａ）の構成において，ランダムパターン１からの放射光は集光レンズ３で集光され，アレイ素子光検出器５に入射される。また，その放射光の一部はビームスプリッタ４で分離され，ピンホール７で視野を狭められて単素子光検出器６に入射される。
【００２９】
アレイ素子光検出器５はシャッターを有しており，ランダムパターン１の移動速度に対して十分な短時間露光が可能である。ピンホール７の位置は，ランダムパターン１のうち，アレイ素子光検出器５の中心のアレイ素子で観測される部分が単素子光検出器でも観測されるように調整されている。アレイ素子（アレイ素子光検出器５のアレイ素子）と単素子（単素子光検出器６の単素子）が同期してデータを取得する。一回の露光で，時間毎のアレイ素子のフレーム画像と遅延時間τの異なる単素子の一組のデータ対ができる。
【００３０】
図１（ｂ）を参照して，図１（ａ）の観測装置によりランダムパターン１を観測する方法を説明する。図１（ｂ）はアレイ素子と単素子が同期して取得した一組のデータ対を表している。図１（ｂ）において，１５はアレイ素子光検出器５の一回の露光で得られたフレーム画像で，１６はフレーム画像１５と同期して取得された単素子光検出器のデータ列である。
【００３１】
図１（ｂ）のフレーム画像１５に対して，次のフレーム画像を取得するまでの間に実際に二次元アレイに投影される画像はたえず時間とともに変化している。単素子光検出器６はその間のフレーム画像の中心部の光強度変化を記録し続け，一連のデータ列を取得する。図１（ｂ）はこのフレーム画像１５と同期して取得された単素子光検出器のデータ列との一組のデータ対の関係を示している。
【００３２】
２５は計測光の像Ａ_０ であり，ランダムパターン１の一部Ｒ_０ が，アレイ素子光検出器５上に投影されたものである。
【００３３】
Ｒ_０ は，遅延時間τ経過した後には移動して，ピンホール７を介して単素子光検出器６によって参照光Ｂ_０ として検出される。
【００３４】
ｘ_ｖｐはアレイ素子光検出器５上の中心アレイから測ったアレイ素子の位置を示すベクトルである（以下添字のｖはその参照記号がベクトル量を表していることを示す）。
【００３５】
アレイ素子光検出器のフレーム画像１５は，遅延時間τ＝０の画像である。遅延時間τにおける単素子光検出器６のデータＢ_０ （光強度）と遅延時間０におけるアレイ素子光検出器５の各アレイの各画素との相関をとると，データＢ_０ の光強度と像Ａ_０ がもっとも相関が高いとする。このことは，計測光の像Ａ_０ （座標位置ｘ_ｖｐ）はτ秒後にはピンホールの位置，即ちアレイ素子光検出器５の中心のアレイ素子に移動することを表している。従って，ランダムパターン１が形状を変えずに移動するものであれば，フレーム画像１５中の計測光の像Ａ_０ と，データ列１６の遅延時間τにおける参照光のデータＢ_０ とは同一の値が記録されている。そこで，このようなデータ対を多数取得して相関計算をすると，フレーム画像１５における計測光の像Ａ_０ の位置に相関ピークが現れる。
【００３６】
このように，パターンの一部Ｒ_０ は，アレイ素子光検出器５の露光時には，座標位置ｘ_ｖｐにあり，遅延時間τ経過した後は，座標位置０_ｖ に移動することを観測できる。従って，図１（ｂ）の観測方法により，アレイ素子光検出器５に投影されるランダムパターンの移動速度と移動方向を観測することができる。
【００３７】
以下において，ＳＤ（Ｓｉｎｇｌｅ Ｄｅｔｅｃｔｏｒ）は単素子光検出器６，ＴＡＤ（Ｔｗｏ Ｄｉｍｅｎｓｉｏｎａｌ Ａｒｒａｙ Ｄｅｔｅｃｔｏｒ）はアレイ素子光検出器５を表す。ランダムパターン１は形を変えずに速度ｖ_ｖ （速さｖのベクトル）で視線に垂直な面を移動するものとする。ランダムパターンの強度分布はガウス分布（Ｎ（ｍ_Ｉ ，σ_Ｉ ^２ ））であるとする。ここでｍ_Ｉ は平均，σ_Ｉ は標準偏差である。パターン（前記のランダムパターン，以下同様）は，ＴＡＤ上に焦点が合わされ，隣接するピクセルと統計的に独立しているとする。入射光の半分はビームスプリッタ４で分離されてＳＤに投射される。ＳＤに入射される光の視野は，ＴＡＤの中央画素の領域に照射されるのと同じになるように制限される。
【００３８】
ＴＡＤはＭ画素×Ｍ画素であり，画素の幅はΔｄである。測定可能な最大変移はｄ_ｍａｘ ＝（Ｍ／２）Δｄである。画素間の感度の差はあらかじめ補正しておく。アレイ素子光検出器５および単素子光検出器６の露光時間Δｔは十分に短く，その間のランダムパターンの移動距離は１画素の幅以下であるようにする。フレーム間隔ｔ_ｍａｘ は十分に長く，アレイ素子光検出器５はその移動を検出できるものとする。ＳＤはＴＡＤの露光の開始に同期してデータサンプリングを開始し，次のフレームの開始までサンプリングを継続する。ＳＤもＴＡＤもパターンの強度変化を検出するのに必要な分解能を備えている。ここで，最小測定速度および最高測定速度はそれぞれ，ｖ_ｖｍｉｎ＝Δｄ／ｔ_ｍａｘ ，およびｖ_ｖｍａｘ＝ｄ_ｍａｘ ／で表される。速度分解能Δｖ_ｖ ／ｖ_ｖ はΔｔ／ｔ_ｍａｘ （ｖ_ｖ ＜Δｄ／Δｔの時）およびΔｄ／ｄ_ｍａｘ （ｖ_ｖ ＞Δｄ／Δｔの時）である。
【００３９】
ＴＡＤの画像フレーム（以下フレームと称する）の画素とＳＤのデータ列のデータ対の関係は，図１（ｂ）で前述した通りである。参照光（参照信号Ａに同じ）と計測光（計測信号Ｂに同じ）の平均をとるのに必要な十分な数の参照光と計測光のデータ対が測定される。画面の位置ｘ_ｖ の光強度をＩ（ｘ_ｖ ，ｔ）であらわす（ｘ_ｖ はベクトルである）。ｘ_ｖ はＴＡＤの画素の座標に対応する。座標ｘ_ｖ の原点はＴＡＤの中心とする。ＳＤはＴＡＤの画素のアレイの原点の光強度を観測する。ＴＡＤのｉ番目のフレームの座標ｘ_ｖ ，時刻ｔ_ｉ の画素値をＩ^ｉ （ｘ_ｖ ，ｔ_ｉ ）で表す。ｉ番目のフレームが開始されてからτ秒後のＳＤで記録されるデータの値をＩ^ｉ _ＳＤ（τ）で表す。ｉ番目のデータ対（計測光の各アレイ素子と参照光とのデータ対のことであって，以下同じ）は，Ｉ^ｉ _ＳＤ（τ）＝Ｉ（０，ｔ_ｉ ＋τ）であり，ここに，ｔ_ｉ は，ｉ番目のフレームの開始時刻である。τはｉフレームが開始されてからの遅延時間である。測定を実行した結果として，ＴＡＤの各画素の強度変化はＩ_ＴＡＤ （ｘ_ｖ ）＝｛Ｉ^ｉ _ＴＡＤ ｝（ｉ＝１，２，・・・のＩ_ＴＡＤ の集合），およびＳＤにおいてＩ_ＳＤ＝｛Ｉ^ｉ _ＳＤ（τ）｝（ｉ＝１，２，・・・）の各フレームにおける遅延時間τのＩ_ＳＤの集合）が得られる。位置ｘ_ｖｐに対して，Ｉ_ＴＡＤ （ｘ_ｖｐ）＝Ｉ_ＳＤ（τ）があるτに対して求まると，速度ｖ_ｖ は，ｖ_ｖ ＝−ｘ_ｖ ／τで求まる。ｘ_ｖｐを特定するためにＩ_ＴＡＤ （ｘ_ｖ ）とＩ_ＳＤ（τ）の相関を計算する。実際の測定では，必要なデータに計測雑音が含まれる。簡単化のために，検出器ＳＤとＴＡＤの各画素が独立に雑音を生成し，その分布は平均０，分散σ_Ｅ ^２ のガウス分布Ｎ（０，σ_Ｅ ^２ ）であるとする。
【００４０】
画素ｘ_ｖｐの光強度Ｉ_０ ＝Ｉ_ＴＡＤ （ｘ_ｖｐ）とし，他の画素の光強度Ｉ_Ｘ ＝Ｉ_ＴＡＤ （ｘ_ｖ ≠ｘ_ｖｐ）とする。パターンは形状を変えないで移動しているとする。装置に雑音がなければ，Ｉ_０ は対になったＳＤのデータ列の中で遅延時間τの光強度Ｉ_ＳＤ（τ）に等しい。しかし，計測雑音のために，ＴＡＤによるＩ_０ とＳＤの値とに違いを生じる。相関のある画素ｘ_ｖｐの強度をＩ_Ｐ とする。ＴＡＤのそれ以外の画素での強度をＩ_Ｂ とする（バックグラウンドレベル）。遅延時間τのＳＤのデータをＩ_Ｔ とする。それぞれの計測雑音をＥ_ｐ ，Ｅ_Ｂ およびＥ_Ｔ とする。この時，これらは次のように表される。
【００４１】
Ｉ_Ｐ ＝Ｉ_０ ＋Ｅ_ｐ
Ｉ_Ｂ ＝Ｉ_Ｘ ＋Ｅ_Ｂ （１）
Ｉ_Ｔ ＝Ｉ_０ ＋Ｅ_Ｔ
式（１）において，各変数のｉ番目のデータは，Ａ^ｉ のように，上付き添字で表す。
【００４２】
上記において，相関のある座標での強度（ピークレベル）とバックグラウンドでの強度では差がある。ピークレベルはＩ_Ｐ とＩ_Ｔ の相関であり，バックグラウンドレベルはＩ_Ｂ とＩ_Ｔ の相関である。ピークレベルとバックグラウンドの大きさの雑音による変動は，相関強度の決定にエラーを生じる。
【００４３】
本発明の計算方法について，二次元モデルを例として相関の求め方について説明する。
【００４４】
（１）一般的方法（従来の方法）
二次元の相関は次の式のようになる。
【００４５】

相関強度のピークは次の式で表される。
【００４６】
ここでは，あとの計算の単純化のために，前述した従来技術の説明におけるＣの計算をする前に平均値をひいてある。
【００４７】
Ｃ_ｐ ＝＜（Ｉ^ｉ _Ｐ −＜Ｉ^ｉ _Ｐ ＞）×（Ｉ^ｉ _Ｔ −＜Ｉ^ｉ _Ｔ ＞）＞ （３）
Ｅ_ｐ およびＥ_Ｔ は分散σ_Ｅ ^２ のガウス分布であるとする。また，計測対象のランダムパターンは分散σ_Ｉ ^２ のガウス分布であるとすると，式（１）および雑音がガウス分布であることから，Ｉ_Ｐ とＩ_Ｔ の分布はガウス分布Ｎ（ｍ_Ｉ ，σ_Ｉ ^２ ＋σ_Ｅ ^２ ）で表される。ここにｍ_Ｉ は平均であり，σ_Ｉ ^２ ＋σ_Ｅ ^２ は分散である。二つの信号Ｉ_Ｐ とＩ_Ｔ は式（１）の共通項Ｉ_０ によって統計的に依存している。Ｃ_ｐ の期待値Ｃ_ｐＥと非バイアス標準偏差σ_ＣＰはＮデータペアに対して
Ｃ_ｐＥ＝σ_Ｉ ^２ （４）
【００４８】
【数３】

【００４９】
である。
【００５０】
一方，バックグラウンドレベルに対しては，
Ｃ_Ｂ ＝＜（Ｉ^ｉ _Ｂ −＜Ｉ^ｉ _Ｂ ＞）×（Ｉ^ｉ _Ｔ −＜Ｉ^ｉ _Ｔ ＞）＞ （６）
Ｉ_Ｂ とＩ_Ｔ は独立であるから，Ｃ_Ｂ の期待値Ｃ_ＢＥと非バイアス標準偏差σ_ＣＢはＮ個のデータペアに対して次のようになる。
【００５１】
Ｃ_ＢＥ＝０ （７）
【００５２】
【数４】

【００５３】
結果的に，一つのデータ対当たりに規格化されたピークのＳＮ比は，次のようになる。
【００５４】
【数５】

【００５５】
ここに，^＊ は規格化されたことを表す。
【００５６】
σ_Ｅ ＜σ_Ｉ の時，ＳＮＲ_Ｇ は，光強度の変動が大きくなるにつれ，１／３^１／２ で一定になる（後述する）。ＳＮ比を向上させるにはＮを大きくする必要がある。その時，ＳＮＲ_Ｇ は，ほぼＳＮＲ_Ｇ Ｎ^１／２ に比例するように修正される。
（２） 排他的二値化法
本発明の最適化排他的二値化法（ＥＢＭ）による相関計算方法に従えば，上記のような一般的方法でのランダムパターンのＳＮＲ_Ｇ の限界を越えることが可能になる。ＥＢＭは，ＴＡＤの画素のデータのうち，ペアになるＳＤの値と等しくないもの（相関のない画素）のデータを除外する（画素値を０にする）。各ＴＡＤのフレームでのＩ^ｉ _ＴＡＤ （ｘ_ｖ ）は，ペアのＳＤの記録データであるＩ^ｉ _ＳＤ（τ）と比較され，次のように二値化される。即ち，両者が同じ値の時，画素値を１し，両者の値が異なる時，画素値を０にする。実際は雑音の影響があるので，両者が全く同じであることはないので，本発明では，両者の大きさがある許容誤差範囲δの範囲にある時に両者は等しいとする。そのため，本発明の排他的二値化画像は次のように表される。
【００５７】
Ｉ^ｉ _ＳＤ（τ）−δ≦Ｉ^ｉ _ＴＡＤ （ｘ_ｖ ）≦Ｉ^ｉ _ＳＤ（τ）＋δの時
Ｂ^ｉ （ｘ_ｖ ）＝１ （１０）
ここで，δ＝０として，Ｂ^ｉ （ｘ_ｖ ）＝１の代わりに，Ｂ^ｉ （ｘ_ｖ ）＝Ａ^２ とすると，前述した従来の技術の説明における排他的相関Ｃ’（式２０）と同じになる。
【００５８】
それ以外の時
Ｂ^ｉ （ｘ_ｖ ）＝０ （１０）’
この方法の場合，許容誤差幅δのとり方でＳＮ比（相関ピーク値とバックグラウンドとの相関値の割合）が変化するので，本発明では許容誤差幅δを，以下に説明するようにランダムパターンの変動の大きさと計測雑音の割合を考慮して最適になるように定める。
【００５９】
二次元の相関強度Ｂ（ｘ_ｖ ）は排他的二値化画像｛Ｂ^ｉ （ｘ_ｖ ）｝の平均で求められる。即ち，
Ｂ（ｘ_ｖ ）＝＜Ｂ^ｉ （ｘ_ｖ ）＞ （１１）
この式は，一般的方法の式（２）に相当する。相関ピークはＢ（ｘ_ｖ ）のピークとして現れる。
【００６０】
ＥＢＭでは，相関強度のピークは，強度Ｉ^ｉ _Ｐ が正しくＩ^ｉ _Ｔ に等しいとして認識される確率Ｐ_Ｐ である。雑音強度を除いて，Ｉ^ｉ _Ｐ とＩ^ｉ _Ｔ に含まれるＩ^ｉ _０ は等しい。誤りは，Ｅ^ｉ _ｐ とＥ^ｉ _Ｔ が許容誤差範囲δより大きい時に生じる。従って，この誤差を生じる時の確率は次の式で表される。
【００６１】
【数６】

【００６２】
ここに，Ｎ（０，２σ_Ｅ ^２ ）は平均０，分散２σ_Ｅ ^２ の合成ガウス分布である。
【００６３】
一般的方法と異なり，Ｐ_Ｐ はσ_Ｉ の依存性がない。δ〜σ_Ｅ の時，式（１２）の確率密度のほとんど全部をカバーしている。この時，Ｐ_Ｐ はほぼ１に等しい。その期待値Ｂ_ＰＥと非バイアスの標準偏差σ_ＢＰは次のようになる。
【００６４】
Ｂ_ＰＥ＝Ｐ_Ｐ （１３）
【００６５】
【数７】

【００６６】
一方，非相関信号であるべき強度Ｉ^ｉ _Ｂ が雑音のためにＩ^ｉ _０ 程度で現れ，相関がある信号として処理されることがある。この誤認識する確率Ｐ_Ｂ は，Ｉ^ｉ _Ｂ ＝（Ｉ^ｉ _Ｘ ＋Ｅ^ｉ _Ｂ ）とＩ^ｉ _Ｔ ＝（Ｉ^ｉ _０ ＋Ｅ^ｉ _Ｔ ）の差が許容誤差範囲δの中にある確率と同じである。その時，確率は，Ｉ_０ ，Ｉ_Ｘ ，Ｅ_Ｂ およびＥ_Ｔ の合成確率密度関数を使用して次のように表される。
【００６７】
【数８】

【００６８】
ここに，Ｎ（０，２σ_Ｉ ^２ ＋２σ_Ｅ ^２ ）は，合成ガウス確率密度関数である。平均値は０，偏差は２σ_Ｉ ^２ ＋２σ_Ｅ ^２ である。δ〜σ_Ｅ ≪σ_Ｉ （即ち，δ≪２σ_Ｉ ^２ ＋２σ_Ｅ ^２ である）。式（１５）において，Ｐ_Ｂ は確率密度関数のＩ＝０付近の小さい部分である。そのためＰ_Ｂ ＝０に近い。Ｎ個のデータペアに対するその期待値Ｂ_ＢＥと非バイアス標準偏差σ_ＢＢは次のようになる。
【００６９】
Ｂ_ＢＥ＝Ｐ_Ｂ （１６）
【００７０】
【数９】

【００７１】
式（９）のように，一つのデータ対当たりに規格化されたピークのＳＮ比は次のようになる。
【００７２】
【数１０】

【００７３】
Ｎデータペアに対するＳＮ比は，Ｎが増加するにつれＳＮＲ_Ｂ はＳＮＲ_Ｂ ^＊ ×Ｎ^１／２ に比例する。
【００７４】
ここで重要な点は，Ｐ_Ｐ とＰ_Ｂ がランダムパターンの光強度の変動σ_Ｉ に影響されないことである。これらは，主に計測雑音σ_Ｅ にのみ決められる。式（１２）からＰ_Ｐ はσ_Ｅ と独立である。Ｐ_Ｂ は光の強度変動に影響されるが（式１５），確率はδにより調整できる。計測雑音σ_Ｅ は光強度の変動σ_Ｉ より一般的に小さい。従って，排他的二値化によってＳＮ比を向上させることが可能である。実際の装置にＥＢＭを適用する前に，許容誤差幅δは，ある与えられたσ_Ｉ とσ_Ｅ に関して最良のＳＮ比になるように最適化しておく。ＳＮ比を検討する場合，あらかじめ概算されたパラメータσ_Ｉ ／σ_Ｅ とδ／σ_Ｅ が利用できる。というのは，システムが改良されてσ_Ｅ が減少しても，δ／σ_Ｅ を一定になるようにすることにより，式（１２）のＰ_Ｐ は同じになるからである。最適化は式（１２）と式（１５）を使用して式（１８）を数値計算し，ＳＮＲ_Ｂ ^＊ が大きくなるようにすることにより達成できる。σ_Ｉ ／σ_Ｅ に関してパラメータδ／σ_Ｅ の最適値が図２（ａ）の表の第２欄とに挙げられ，図２（ｂ）にプロットされている。最適値はσ_Ｉ とσ_Ｅ に弱く依存している。３から１００のσ_Ｉ ／σ_Ｅ の現実的な場合に，δ／σ_Ｅ は２．５から４．３に変化する。
【００７５】
（３）一般的方法と排他的二値化法との比較
排他的二値化法と一般的方法とにより計算されたそれぞれのＳＮ比のグラフが図３（ａ）と図２（ａ）の表の第３欄と第４欄に示されている。図３（ａ）のグラフにおいて，ＥＢＭはσ_Ｉ ／σ_Ｅ に対してＳＮ比が最大になるδ／σ_Ｅ を用いたものである。ＥＢＭにより得られる最適ＳＮ比はσ_Ｉ ／σ_Ｅ が増加するにつれて増加するが，一般的な方法のＳＮ比はσ_Ｉ ／σ_Ｅ ＞２で一定（＝１／３^１／２ ）である。現実的な装置のσ_Ｉ ／σ_Ｅ は，２から１００程度あるので，ＳＮＲ_Ｂ ^＊ は０．６から６．０に増加する。図３（ｂ）は排他的二値化法と一般法とでのＳＮ比の割合を示す。その割合はσ_Ｉ ／σ_Ｅ の平方根にほぼ比例して増加するが，データ数には依存しない。σ_Ｉ ／σ_Ｅ ≧２^１／２ のとき，ＳＮＲ_Ｂ ^＊ はＳＮＲ_Ｇ より大きく，例えば，ＳＮＲ_Ｂ ／ＳＮＲ_Ｇ ＝３．３と１０．３（σ_Ｉ ／σ_Ｅ ＝１０と１００）である。しかし，σ_Ｉ ／σ_Ｅ ≦２^１／２ のとき，ＳＮＲ_Ｂ はＳＮＲ_Ｇ より小さい。このことから分かるように，良いＳＮ比で測定するためには，計測雑音を減少させることがＥＢＭでは重要である。
【００７６】
シミュレーションの例
ＥＢＭの有効性を示すために，前出のモデルで説明された場合についてシミュレーションを行なった。次のパラメータが仮定されている。ＴＡＤ画素は１１×１１（Ｍ＝１１）である。画素の大きさはΔｄ＝１ｍｍである。フレーム間隔はｔ_ｍａｘ ＝５ｍｓである。ＳＤのサンプリング間隔はΔｔ＝１ｍｓである。各画素の光強度は統計的に独立であり，ガウス分布しているとする（平均ｍ_Ｉ ＝１０，標準偏差σ_Ｉ ＝１，Ｎ（１０，１））。統計的に独立な計測雑音を，ＴＡＤとＳＤの各画素により得られたデータに付加した。その確率密度関数は平均０のガウス関数である。計測雑音の標準偏差に対して２つの場合が考えられる。σ_Ｅ ＝０．１（すなわち，σ_Ｉ ／σ_Ｅ ＝１０）と他の場合はσ_Ｅ ＝０．０１（即ち，σ_Ｉ ／σ_Ｅ ＝１００）である。雑音は光強度と統計的に独立であると仮定されている。ＥＢＭのパラメータδ／σ_Ｅ はσ_Ｉ ／σ_Ｅ ＝１０およびσ_Ｉ ／σ_Ｅ に対して，図２（ａ）の表と図２（ｂ）を参照して，最適値として３．２と４．１をセットする。ランダムパターンは速度ｖ_ｖ ＝（ｖ_ｘ ，ｖ_ｙ ），ｖ_ｘ ＝−１ｍ／ｓ，ｖ_ｙ ＝−１ｍ／ｓで移動する。ＴＡＤフレームとＳＤデータの２５ペアが同期的に得られ，以下の相関計算に使用された。
【００７７】
図４（ａ）と（ｂ）は，遅延時間０の状態の排他的二値化法による場合と一般的方法による場合の３次元プロットを示している。図４（ａ）は排他的二値化法の場合であり，図４（ｂ）は一般的方法の場合である。σ_Ｉ ／σ_Ｅ ＝１０の場合である。遅延時間０の二次元の相関をＥＢＭの式（１１）と一般的方法の式（２）を使用して計算する。これらの図は次のように規格化されている。相関のある場合は１に規格化した，また，相関のない期待値は０に規格化した。
【００７８】
ＥＢＭにより決定されたピーク高さのＳＮ比は，一般的方法より良くなっている。ＳＮＲ_Ｂ ＝１１であり，ＳＮＲ_Ｇ ＝２．９である（図２（ａ）の表もしくは図２（ｂ）の評価は，９．０と２．８である）。
【００７９】
図５（ａ），（ｂ）は，それぞれ図４（ａ），（ｂ）の遅延時間τ＝３ｍｓの二次元の相関強度である。そのパターンは３ｍｓ後に中心に移動するように現れている。これらの結果は，ＥＢＭでは二次元の相関においてはバックグラウンドレベルの変動が減少している。そのために，一般的方法より明瞭にピーク相関位置の動きを決定することができることを示している。
【００８０】
図６と図７は，図４と図５と同様のものを示していて，σ_Ｉ ／σ_Ｅ ＝１００の場合である。図６（ａ）は排他的二値化方法の初期状態を示し，図７（ａ）は排他的二値化方法の遅延時間τ＝３ｍｓの状態を示す。図６（ｂ）は一般的方法の初期状態を示し，図７（ｂ）は一般的方法の遅延時間τ＝３ｍｓの状態を示す。相関期待値は，Ｂ_ＰＥ＝１，０，Ｂ_ＢＥ＝０．０２である（σ_Ｉ ／σ_Ｅ ＝１００）。ＥＢＭの場合にはＳＮ比が向上している。ＳＮＲ_Ｂ ＝２９であり，ＳＮＲ_Ｇ ＝２．８である。図２（ａ）の表と図２（ｂ）のグラフによる評価では２９と２．８である。
【００８１】
図８は，本発明の装置構成の実施の形態を示す図である。
【００８２】
図８において，５はアレイ素子光検出器（ＴＡＤ）である。６は単素子光検出器（ＳＤ）である。２１は入力インタフェースであって，アレイ素子光検出器（ＴＡＤ）５，単素子光検出器（ＳＤ）６の入力インタフェースである。２２はコンピュータである。２３は入力部であって，キーボード等である。２４は出力部であって，ディスプレイ，プリンタ等である。
【００８３】
３１はパラメータ保持部であって，相関演算に必要な各種パラメータを保持するものである。３２は相関演算部であって，アレイ素子光検出器５および単素子光検出器６のデータおよび各種パラメータにより計測信号Ａと参照信号Ｂの相関を求めるものである。相関演算部３２は二値化されたフレーム画像の画素値を画素毎に累積して保持する二値化相関バッファを備える。３３はデータ保持部，アレイ素子光検出器５および単素子光検出器６から入力されるデータを保持するものである。３４は同期制御部であって，アレイ素子光検出器５および単素子光検出器６から入力されるデータの同期をとるものである。３５は最適パラメータ演算部であって，各種パラメータに従って，ＳＮ比が最大にあるようなδ／σ_Ｅ を求めるものである。
【００８４】
図２の構成において，σ_Ｉ ，σ_Ｅ ，最適なδ等の各種パラメータが入力部２３から入力され，パラメータ保持部３１に保持される。σ_Ｉ とσ_Ｅ は予め計測して求めておくか，あるいはＳＤの遅延時間０の記録データを使用して求める。この場合，データ保持部３３の記録データをもとに，最適パラメータ演算部３５で計算させることができる。同期制御部３４はアレイ素子光検出器５，単素子光検出器６から時系列に沿って入力されるデータの同期をフレーム毎にとる。同期をとって入力されたそれぞれのデータはデータ保持部３３に保持される。相関演算部３２はデータ保持部３３に保持されているアレイ素子光検出器５および単素子光検出器６からの入力データと，パラメータ保持部３１に保持されているパラメータに従って，前述の排他的二値化方法により排他的二値化を行ない，参照信号と計測信号の相関演算を行ない，相関ピークを求める。求められた相関ピークの座標，遅延時間等は出力部２４に出力される。
【００８５】
入力部２３から最適パラメータを求めるのに必要な各種パラメータが入力される。σ_Ｉ とσ_Ｅ は，あらかじめ求めておいたものを入力するか，あるいはデータ保持部３３の記録データをもとに最適パラメータ演算部３５で求めるようにしても良い。最適パラメータ演算部３５は，各種パラメータに従って，σ_Ｉ ／σ_Ｅ に対してＳＮ比が最大になるように最適化されたδ／σ_Ｅ を求める。その演算結果は出力部２４に出力される。
【００８６】
図９は，本発明のデータ記録のフローチャートである。図９は図１（ａ），（ｂ）の構成において，ＴＡＤのフレーム画像とＳＤのデータを同期的に取得する動作のフローチャートである。図１（ａ）のアレイ素子光検出器５（ＴＡＤ）とビームスプリッタ４の間にシャッタがあって，アレイ素子光検出器５の露光が終了した後も，単素子光検出器６（ＳＤ）はデータを取得し続けることができる。
【００８７】
Ｓ１でｉ＝１に初期化する。Ｓ２でＴＡＤの露光，ＳＤのデータ取得を同期して開始する。Ｓ３で露光時間待つ（この間，ＴＡＤの画像データを取得するとともにＳＤの一連のデータをとる）Ｓ４で露光時間を終了する。Ｓ５で所定の時間待つ（τの最大値より大きい時間とる）。この間，ＳＤは一連のデータを録りつづける。Ｓ６で一連のＳＤデータの取得を終了する。Ｓ７でＴＡＤのフレーム画像，ＳＤのデータ列を記録して保存する。Ｓ８で予定数の画像フレームのデータが取得できたか判定し，予定数のデータが取得できていれば，処理を終了し，予定数のデータがとれていなければＳ９でｉを更新し，Ｓ２以降の処理を繰り返す。
【００８８】
図１０は本発明の相関演算のフローチャートである。
【００８９】
Ｓ１ 二値化相関値バッファの画素値を全て０に初期化する（ｉ＝１）。
【００９０】
Ｓ２でｉフレームのＴＡＤとＳＤのデータ対を入力する。Ｓ３でＳＤの記録データとＴＡＤの各画素値を比較する。Ｓ４とＳ５で画素値が（Ｉ^ｉ _ＳＤ−δ）と（Ｉ^ｉ _ＳＤ＋δ）の間にあるか，ないかを判定する。画素値が（Ｉ^ｉ _ＳＤ−δ）と（Ｉ^ｉ _ＳＤ＋δ）の間にあれば，Ｓ６で画素値を１にする。間になければ，Ｓ７で画素値を０にする。Ｓ８で二値化されたＴＡＤのフレーム画像の画素値を二値化相関値バッファに加算する。Ｓ９でさらにデータ対を追加するか判定し，さらに追加するのであれば，Ｓ１０でｉ＝ｉ＋１とし，Ｓ２以降の処理を繰り返す。なお，本発明は，取得したデータを保存してから相関計算する方法だけでなく，データを取得しながら相関計算をすることも可能である。従って，Ｓ９における追加は，データ保持部３３に記録されたデータ対を読み出すか，あるいは新たなデータ対を取得することを意味する。Ｓ９の判断でさらに追加しないのであれば，Ｓ１１で二値化相関値バッファの各画素値をｉで割り，平均をとる。Ｓ１２で相関結果を出力して，処理を終了する。
【００９１】
図１１は本発明の最適パラメータ決定のフローチャートである。
【００９２】
Ｓ１でσ_Ｉ とσ_Ｅ を定める。Ｓ２で，δ／σ_Ｅ の値を変えながら各δ／σ_Ｅ に対してＰ_Ｂ （式１２）とＰ_Ｂ （式１５）を数値計算する。Ｓ３でＳＮ比（式１８）を計算し，σ_Ｉ ／σ_Ｅ 毎にＳＮ比を最大するδ／σ_Ｅ を選ぶ。
【００９３】
Ｓ４で最適パラメータを出力する。
【００９４】
【発明の効果】
本発明によれば，参照信号と計測信号の大きさの違いに基づく二値化を，最適化された誤差許容幅に従って，相関ピークの鮮明さを示すＳＮ比が最大になるように行なうことができる。そのため，高い精度で，相関ピークを決定することができる。そのため，高速度で移動する物体の動き等を高速で正確に測定できる。そのため，地球大気のゆらぎの観測等に利用することも可能で，地球大気のゆらぎ層の風速，風向を正確に観測することが可能になる。
【図面の簡単な説明】
【図１】本発明の実施の形態１を示す図である。
【図２】本発明の結果を示すグラフ１を示す図である。
【図３】本発明の結果を示すグラフ２を示す図である。
【図４】本発明の相関結果の三次元図１である。
【図５】本発明の相関結果の三次元図２である。
【図６】本発明の相関結果の三次元図３である。
【図７】本発明の相関結果の三次元図４である。
【図８】本発明の実施の形態２を示す図である。
【図９】本発明のデータ記録のフローチャートを示す図である。
【図１０】本発明の相関演算のフローチャートを示す図である。
【図１１】本発明の最適パラメータ決定のフローチャートを示す図である。
【図１２】従来の技術の説明図である。
【符号の説明】
１：ランダムパターン
３：集光レンズ
４：ビームスプリッタ
５：アレイ素子光検出器
６：単素子光検出器
７：ピンホール
８：計測光
９：参照光

Claims (9)

1. 参照信号と移動する計測対象の信号である計測信号の相関を計測する方法において，
   参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をし，参照信号と計測信号の相関を計測することを特徴とする最適化排他的二値化相関計測方法。
2. 許容誤差幅は，計測環境に応じて定めるものであり，参照信号と計測信号の差と最適化された該許容誤差幅の大小関係に応じて排他的二値化を行なうものであることを特徴とする請求項１に記載の最適化排他的二値化相関計測方法。
3. 相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように該許容誤差幅を最適化することを特徴とする請求項１または２に記載の最適化排他的二値化相関計測方法。
4. 計測対象は形状を変えずに移動するランダムなパターンであって，計測信号は計測対象の二次元画像の信号であり，参照信号は該計測信号の二次元画像の限られた領域の画像の信号であって計測信号と同期して取得されるものであり，参照信号と計測信号の相関により計測対象の動きを計測するものであることを特徴とする請求項１，２または３に記載の最適化排他的二値化相関計測方法。
5. 計測光を受光するアレイ素子光検出器と計測光の一部を受光する単素子光検出器と，アレイ素子光検出器から出力される計測信号と，単素子光検出器の出力する参照信号の相関を取る相関演算部とを備えた最適化排他的二値化相関計測装置であって，
   該相関演算部は，参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と許容誤差幅との関係により計測信号の排他的二値化を最適化し，参照信号と計測信号の相関を計測することを特徴とする最適化排他的二値化相関計測装置。
6. 許容誤差幅は，計測環境に応じて定めるものであり，参照信号と計測信号の差と最適化された該許容誤差幅の大小関係に応じて排他的二値化を行なうものであることを特徴とする請求項５に記載の最適化排他的二値化相関計測装置。
7. 相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように該許容誤差幅を最適化することを特徴とする請求項５または６に記載の最適化排他的二値化相関計測装置。
8. 参照信号と移動する計測対象の信号である計測信号の相関を計測するプログラムにおいて，
   参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をし，参照信号と計測信号の相関を計測することにより最適化排他的二値化相関計測を実行することを特徴とするプログラム。
9. 最適化排他的二値化相関計算を最適化するプログラムにおいて，参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をするものであって，
   許容誤差幅は，計測環境に応じて定めるものであり，
   相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように最適化された該許容誤差幅を求めることを特徴とするプログラム。

Images