(57)【特許請求の範囲】
【請求項１】
参照信号と移動する計測対象の信号である計測信号の相関を計測する方法において，
参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をし，参照信号と計測信号の相関を計測することを特徴とする最適化排他的二値化相関計測方法。
【請求項２】
許容誤差幅は，計測環境に応じて定めるものであり，参照信号と計測信号の差と最適化された該許容誤差幅の大小関係に応じて排他的二値化を行なうものであることを特徴とする請求項１に記載の最適化排他的二値化相関計測方法。
【請求項３】
相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように該許容誤差幅を最適化することを特徴とする請求項１または２に記載の最適化排他的二値化相関計測方法。
【請求項４】
計測対象は形状を変えずに移動するランダムなパターンであって，計測信号は計測対象の二次元画像の信号であり，参照信号は該計測信号の二次元画像の限られた領域の画像の信号であって計測信号と同期して取得されるものであり，参照信号と計測信号の相関により計測対象の動きを計測するものであることを特徴とする請求項１，２または３に記載の最適化排他的二値化相関計測方法。
【請求項５】
計測光を受光するアレイ素子光検出器と計測光の一部を受光する単素子光検出器と，アレイ素子光検出器から出力される計測信号と，単素子光検出器の出力する参照信号の相関を取る相関演算部とを備えた最適化排他的二値化相関計測装置であって，
該相関演算部は，参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と許容誤差幅との関係により計測信号の排他的二値化を最適化し，参照信号と計測信号の相関を計測することを特徴とする最適化排他的二値化相関計測装置。
【請求項６】
許容誤差幅は，計測環境に応じて定めるものであり，参照信号と計測信号の差と最適化された該許容誤差幅の大小関係に応じて排他的二値化を行なうものであることを特徴とする請求項５に記載の最適化排他的二値化相関計測装置。
【請求項７】
相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように該許容誤差幅を最適化することを特徴とする請求項５または６に記載の最適化排他的二値化相関計測装置。
【請求項８】
参照信号と移動する計測対象の信号である計測信号の相関を計測するプログラムにおいて，
参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をし，参照信号と計測信号の相関を計測することにより最適化排他的二値化相関計測を実行することを特徴とするプログラム。
【請求項９】
最適化排他的二値化相関計算を最適化するプログラムにおいて，参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をするものであって，
許容誤差幅は，計測環境に応じて定めるものであり，
相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように最適化された該許容誤差幅を求めることを特徴とするプログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は，最適化排他的二値化相関計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラムに関するものである。特に，移動する計測対象の信号である計測信号と参照信号の相関を計測するものである。高速に移動する物体の移動速度，移動方向等を観測することを可能にするものであり，地球大気のゆらぎの観測等にも利用できるものである。
【０００２】
星のまたたきは，星からの光が地球大気圏の屈折率ゆらぎによって波面（位相）を乱されることによって生じる。この星の瞬きが地表に生成するランダムパターンを観測することによりゆらぎ層における風向，風速，ゆらぎ層の高度，ゆらぎの強度等を計測することができる。このランダムパターンは高速に移動するので，高速な画像処理を必要とする。
【０００３】
【従来の技術】
このような，ランダムパターンの観測には，ランダムパターンの二次元画像の一部の画像信号を参照信号として高速に取得し，ランダムパターンの二次元画像を計測信号として，両者の相関をとることによりランダムパターンを計測する。
【０００４】
図１２（ａ）は計測信号を表し，図１２（ｂ）は参照信号を表す。
Ａは計測対象の計測信号，Ｂは参照信号である。
【０００５】
計測信号Ａと参照信号Ｂの相関をとり，両者が同一の信号であるかどうかを判別する。計測信号Ａと参照信号Ｂは，Ａ＝＜Ａ＞＋α，Ｂ＝＜Ｂ＞＋βで表すことができる。
【０００６】
【数１】
0003554822.tif 000002
【０００７】
【数２】
0003554822.tif 000003
【０００８】
であり，＜Ａ＞，＜Ｂ＞は，それぞれの信号強度の平均値を表す。αは計測信号Ａの平均値からの変動量，βは参照信号Ｂの平均値からの変動を表す。
【０００９】
計測信号Ａと参照信号Ｂとの相関Ｃは次の式で計算する。
【００１０】
0003554822.tif 000004
ここに，＜＞はＡとＢの平均を表す。ＡおよびＢがランダムに変動するものであると，＜α＞＝０，＜β＞＝０である。
【００１１】
また，計測信号Ａと参照信号Ｂが同一でないとき時（Ａ≠Ｂ）
＜αβ＞＝０
計測信号Ａと参照信号Ｂが同一である時（Ａ＝Ｂの時）
＜αβ＞＝σ２  となる。σ２  は計測信号Ａと参照信号の分散である。
【００１２】
従って，Ｃ＝＜Ａ＞２    （Ａ≠Ｂの時）
Ｃ＝＜Ａ＞２  ＋σ２  （Ａ＝Ｂの時）
で表される。
【００１３】
例えば，図１（ａ）の装置でランダムに変動するパターンが，パターン形状を変えないで移動する場合の移動速度を計測すると考える（図１（ａ）の詳細は後述する）。
【００１４】
アレイ素子光検出器５と単素子光検出器６は，同期して複数（ｉ＝１，２，・・・・，Ｎ）のデータ対を取得する。
【００１５】
計測信号Ａはランダムパターン１がアレイ素子光検出器５の任意の一つの画素で観測される信号Ａｉ  の集合である。参照信号Ｂはピンホール７を介して単素子光検出器６で観測される信号で，対となっているアレイ素子光検出器５のデータ取得時刻から遅延時間τだけ遅れて取得された信号Ｂｉ  の集合である。
【００１６】
図１（ａ）の装置で，参照信号と計測信号の一般的な相関Ｃは上記の式のようになる。
【００１７】
その結果からランダムパターン１が形状を変えずに移動している場合には，この相関Ｃを計算することにより，ランダムパターン１の移動する方向と速度を
求めることができる。
【００１８】
この一般的な相関Ｃの計算に対して，相関を計算する前に計測信号Ａのｉ番目の要素と参照信号Ｂのｉ番目の要素を比較して，それらが異なる場合は両方の要素の値をゼロにしてから排他的相関Ｃ’を
Ｃ’＝＜Ａ×Ｂ＞
を計算することにより
Ｃ’＝０  （Ａ≠Ｂ）
Ｃ’＝＜Ａ＞２  ＋σ２  （Ａ＝Ｂの時）      （２０）                  のように相関の無い部分の値を減らすことで，ピークをより鮮明にすることができる。その結果として少ない信号要素数で相関を計算することができる。
【００１９】
【発明が解決しようとする課題】
上記の方法は，後述するように，計測信号か参照信号，あるいはその両方に計測雑音が含まれる場合には，Ａ＝Ｂであった信号が，Ａ＝Ｂ＋ε（εは計測雑音）と変化してしまうためそのまま適用すると上記式（２０）で示されたピークの高さが大幅に減少してしまう。そのため，Ａ＝Ｂ±δの範囲にある場合は，Ａ＝Ｂと見なす許容誤差幅δを導入する必要がある。この許容誤差幅δの値をどのように設定するかは，以下に説明するように非常に重要である。
【００２０】
このδのを適切に設定しないと，例えば，δが小さすぎる場合は，本来はＡ＝Ｂであった信号がＡ≠Ｂと判定される確率が高くなり，相関のピークが減少する。また，δが大きすぎる場合は，本来はＡ≠Ｂであった信号がＡ＝Ｂと判定される確率が高くなり，相関の無い部分の値が増加し，ピークを明瞭に判定できなくなる。そこで，本発明では信号変動幅と計測雑音を考慮して最適な許容誤差幅δを求め，排他的二値化を最適化して相関計測を行なう最適化排他的二値化相関計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラムを提供することを目的とする。
【００２１】
【課題を解決するための手段】
本発明は，参照信号と移動する計測対象の信号である計測信号の相関を計測する最適化排他的二値化計測方法および最適化排他的二値化計測装置およびその方法を実行するプログラムにおいて，参照信号と計測信号のそれぞれの信号値を比較し，参照信号と計測信号の値が同じか否かに応じて計測信号を排他的に二値化するものであり，参照信号と計測信号の値が同じか否かを判定する基準として最適化された許容誤差幅をもたせ，参照信号と計測信号と該許容誤差幅との関係により計測信号に対して排他的二値化をし，参照信号と計測信号の相関を計測するようにした。
【００２２】
そして，許容誤差幅は，計測環境に応じて定めるものである。許容誤差幅は計測環境中の計測雑音に応じて定まるものであり，参照信号と計測信号の差と最適化された該許容誤差幅の大小関係に応じて排他的二値化を行なうようにした。
【００２３】
さらに，相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量をＳＮ比とした時，ＳＮ比が最大になるように該許容誤差幅を最適化するものである。
【００２４】
ここでＳＮ比とは，相関ピークの期待値と相関が無い部分の期待値の差を，これらの期待値のゆらぎの大きさで割った量で，ＳＮ比が大きいほど相関ピークが鮮明であることを示すものである。
【００２５】
本発明は，相関をとる二次元画像の計測信号と参照信号の信号強度が同じであるかないかを判断を，ＳＮ比ができるかぎり良くなるように求めた許容誤差幅に基づいて排他的に二値化して計測信号と参照信号の相関をとるようにしたので，ＳＮ比が良い相関が得られる。また，信号強度の異なる部分は画素値を，例えば，０にし，信号強度が同じ部分の画素値を，例えば，１にするように二値化してから，相関をとるようにしたので，保存データ量も少なくてすみ，簡単な計算で高速に相関処理を行なうことができる。
【００２６】
【発明の実施の形態】
図１は本発明の実施の形態１を示す図である。図１（ａ）は観測装置を示し，図１（ｂ）は観測方法の説明図である。
【００２７】
図１（ａ）において，１は計測対象の光パターンであり，ランダムパターンである。ランダムパターン１は，例えば，高速で移動する物体から放射される光，あるいは星が地表面につくるランダムパターン等である。３は集光レンズであり，ランダムパターン１から放射される光を集光するものである。集光レンズ３は望遠鏡でも良い。４はビームスプリッタであり，ランダムパターン１から放射される光を計測光８と参照光９に分離するものである。５はアレイ素子光検出器であり，例えば，ＣＣＤのような二次元の光検出器である。６は単素子光検出器であり，ピンホール７を通過する狭い視野の光を検出するものである。単素子光検出器６はアレイ素子光検出器５に比べて高速にサンプリングできるものである（例えば，光電子増倍管等を使用する）。７はピンホールである。８は計測光であり，アレイ素子光検出器５に入射する光である。９は参照光であって，単素子光検出器６に入射する光である。
【００２８】
図１（ａ）の構成において，ランダムパターン１からの放射光は集光レンズ３で集光され，アレイ素子光検出器５に入射される。また，その放射光の一部はビームスプリッタ４で分離され，ピンホール７で視野を狭められて単素子光検出器６に入射される。
【００２９】
アレイ素子光検出器５はシャッターを有しており，ランダムパターン１の移動速度に対して十分な短時間露光が可能である。ピンホール７の位置は，ランダムパターン１のうち，アレイ素子光検出器５の中心のアレイ素子で観測される部分が単素子光検出器でも観測されるように調整されている。アレイ素子（アレイ素子光検出器５のアレイ素子）と単素子（単素子光検出器６の単素子）が同期してデータを取得する。一回の露光で，時間毎のアレイ素子のフレーム画像と遅延時間τの異なる単素子の一組のデータ対ができる。
【００３０】
図１（ｂ）を参照して，図１（ａ）の観測装置によりランダムパターン１を観測する方法を説明する。図１（ｂ）はアレイ素子と単素子が同期して取得した一組のデータ対を表している。図１（ｂ）において，１５はアレイ素子光検出器５の一回の露光で得られたフレーム画像で，１６はフレーム画像１５と同期して取得された単素子光検出器のデータ列である。
【００３１】
図１（ｂ）のフレーム画像１５に対して，次のフレーム画像を取得するまでの間に実際に二次元アレイに投影される画像はたえず時間とともに変化している。単素子光検出器６はその間のフレーム画像の中心部の光強度変化を記録し続け，一連のデータ列を取得する。図１（ｂ）はこのフレーム画像１５と同期して取得された単素子光検出器のデータ列との一組のデータ対の関係を示している。
【００３２】
２５は計測光の像Ａ０  であり，ランダムパターン１の一部Ｒ０  が，アレイ素子光検出器５上に投影されたものである。
【００３３】
Ｒ０  は，遅延時間τ経過した後には移動して，ピンホール７を介して単素子光検出器６によって参照光Ｂ０  として検出される。
【００３４】
ｘｖｐはアレイ素子光検出器５上の中心アレイから測ったアレイ素子の位置を示すベクトルである（以下添字のｖはその参照記号がベクトル量を表していることを示す）。
【００３５】
アレイ素子光検出器のフレーム画像１５は，遅延時間τ＝０の画像である。遅延時間τにおける単素子光検出器６のデータＢ０  （光強度）と遅延時間０におけるアレイ素子光検出器５の各アレイの各画素との相関をとると，データＢ０  の光強度と像Ａ０  がもっとも相関が高いとする。このことは，計測光の像Ａ０  （座標位置ｘｖｐ）はτ秒後にはピンホールの位置，即ちアレイ素子光検出器５の中心のアレイ素子に移動することを表している。従って，ランダムパターン１が形状を変えずに移動するものであれば，フレーム画像１５中の計測光の像Ａ０  と，データ列１６の遅延時間τにおける参照光のデータＢ０  とは同一の値が記録されている。そこで，このようなデータ対を多数取得して相関計算をすると，フレーム画像１５における計測光の像Ａ０  の位置に相関ピークが現れる。
【００３６】
このように，パターンの一部Ｒ０  は，アレイ素子光検出器５の露光時には，座標位置ｘｖｐにあり，遅延時間τ経過した後は，座標位置０ｖ  に移動することを観測できる。従って，図１（ｂ）の観測方法により，アレイ素子光検出器５に投影されるランダムパターンの移動速度と移動方向を観測することができる。
【００３７】
以下において，ＳＤ（Ｓｉｎｇｌｅ  Ｄｅｔｅｃｔｏｒ）は単素子光検出器６，ＴＡＤ（Ｔｗｏ  Ｄｉｍｅｎｓｉｏｎａｌ  Ａｒｒａｙ  Ｄｅｔｅｃｔｏｒ）はアレイ素子光検出器５を表す。ランダムパターン１は形を変えずに速度ｖｖ  （速さｖのベクトル）で視線に垂直な面を移動するものとする。ランダムパターンの強度分布はガウス分布（Ｎ（ｍＩ  ，σＩ  ２  ））であるとする。ここでｍＩ  は平均，σＩ  は標準偏差である。パターン（前記のランダムパターン，以下同様）は，ＴＡＤ上に焦点が合わされ，隣接するピクセルと統計的に独立しているとする。入射光の半分はビームスプリッタ４で分離されてＳＤに投射される。ＳＤに入射される光の視野は，ＴＡＤの中央画素の領域に照射されるのと同じになるように制限される。
【００３８】
ＴＡＤはＭ画素×Ｍ画素であり，画素の幅はΔｄである。測定可能な最大変移はｄｍａｘ  ＝（Ｍ／２）Δｄである。画素間の感度の差はあらかじめ補正しておく。アレイ素子光検出器５および単素子光検出器６の露光時間Δｔは十分に短く，その間のランダムパターンの移動距離は１画素の幅以下であるようにする。フレーム間隔ｔｍａｘ  は十分に長く，アレイ素子光検出器５はその移動を検出できるものとする。ＳＤはＴＡＤの露光の開始に同期してデータサンプリングを開始し，次のフレームの開始までサンプリングを継続する。ＳＤもＴＡＤもパターンの強度変化を検出するのに必要な分解能を備えている。ここで，最小測定速度および最高測定速度はそれぞれ，ｖｖｍｉｎ＝Δｄ／ｔｍａｘ  ，およびｖｖｍａｘ＝ｄｍａｘ  ／で表される。速度分解能Δｖｖ  ／ｖｖ  はΔｔ／ｔｍａｘ  （ｖｖ  ＜Δｄ／Δｔの時）およびΔｄ／ｄｍａｘ  （ｖｖ  ＞Δｄ／Δｔの時）である。
【００３９】
ＴＡＤの画像フレーム（以下フレームと称する）の画素とＳＤのデータ列のデータ対の関係は，図１（ｂ）で前述した通りである。参照光（参照信号Ａに同じ）と計測光（計測信号Ｂに同じ）の平均をとるのに必要な十分な数の参照光と計測光のデータ対が測定される。画面の位置ｘｖ  の光強度をＩ（ｘｖ  ，ｔ）であらわす（ｘｖ  はベクトルである）。ｘｖ  はＴＡＤの画素の座標に対応する。座標ｘｖ  の原点はＴＡＤの中心とする。ＳＤはＴＡＤの画素のアレイの原点の光強度を観測する。ＴＡＤのｉ番目のフレームの座標ｘｖ  ，時刻ｔｉ  の画素値をＩｉ  （ｘｖ  ，ｔｉ  ）で表す。ｉ番目のフレームが開始されてからτ秒後のＳＤで記録されるデータの値をＩｉ  ＳＤ（τ）で表す。ｉ番目のデータ対（計測光の各アレイ素子と参照光とのデータ対のことであって，以下同じ）は，Ｉｉ  ＳＤ（τ）＝Ｉ（０，ｔｉ  ＋τ）であり，ここに，ｔｉ  は，ｉ番目のフレームの開始時刻である。τはｉフレームが開始されてからの遅延時間である。測定を実行した結果として，ＴＡＤの各画素の強度変化はＩＴＡＤ  （ｘｖ  ）＝｛Ｉｉ  ＴＡＤ  ｝（ｉ＝１，２，・・・のＩＴＡＤ  の集合），およびＳＤにおいてＩＳＤ＝｛Ｉｉ  ＳＤ（τ）｝（ｉ＝１，２，・・・）の各フレームにおける遅延時間τのＩＳＤの集合）が得られる。位置ｘｖｐに対して，ＩＴＡＤ  （ｘｖｐ）＝ＩＳＤ（τ）があるτに対して求まると，速度ｖｖ  は，ｖｖ  ＝−ｘｖ  ／τで求まる。ｘｖｐを特定するためにＩＴＡＤ  （ｘｖ  ）とＩＳＤ（τ）の相関を計算する。実際の測定では，必要なデータに計測雑音が含まれる。簡単化のために，検出器ＳＤとＴＡＤの各画素が独立に雑音を生成し，その分布は平均０，分散σＥ  ２  のガウス分布Ｎ（０，σＥ  ２  ）であるとする。
【００４０】
画素ｘｖｐの光強度Ｉ０  ＝ＩＴＡＤ  （ｘｖｐ）とし，他の画素の光強度ＩＸ  ＝ＩＴＡＤ  （ｘｖ  ≠ｘｖｐ）とする。パターンは形状を変えないで移動しているとする。装置に雑音がなければ，Ｉ０  は対になったＳＤのデータ列の中で遅延時間τの光強度ＩＳＤ（τ）に等しい。しかし，計測雑音のために，ＴＡＤによるＩ０  とＳＤの値とに違いを生じる。相関のある画素ｘｖｐの強度をＩＰ  とする。ＴＡＤのそれ以外の画素での強度をＩＢ  とする（バックグラウンドレベル）。遅延時間τのＳＤのデータをＩＴ  とする。それぞれの計測雑音をＥｐ  ，ＥＢ  およびＥＴ  とする。この時，これらは次のように表される。
【００４１】
ＩＰ  ＝Ｉ０  ＋Ｅｐ  
ＩＢ  ＝ＩＸ  ＋ＥＢ                  （１）
ＩＴ  ＝Ｉ０  ＋ＥＴ  
式（１）において，各変数のｉ番目のデータは，Ａｉ  のように，上付き添字で表す。
【００４２】
上記において，相関のある座標での強度（ピークレベル）とバックグラウンドでの強度では差がある。ピークレベルはＩＰ  とＩＴ  の相関であり，バックグラウンドレベルはＩＢ  とＩＴ  の相関である。ピークレベルとバックグラウンドの大きさの雑音による変動は，相関強度の決定にエラーを生じる。
【００４３】
本発明の計算方法について，二次元モデルを例として相関の求め方について説明する。
【００４４】
（１）一般的方法（従来の方法）
二次元の相関は次の式のようになる。
【００４５】
0003554822.tif 000005
相関強度のピークは次の式で表される。
【００４６】
ここでは，あとの計算の単純化のために，前述した従来技術の説明におけるＣの計算をする前に平均値をひいてある。
【００４７】
Ｃｐ  ＝＜（Ｉｉ  Ｐ  −＜Ｉｉ  Ｐ  ＞）×（Ｉｉ  Ｔ  −＜Ｉｉ  Ｔ  ＞）＞    （３）
Ｅｐ  およびＥＴ  は分散σＥ  ２  のガウス分布であるとする。また，計測対象のランダムパターンは分散σＩ  ２  のガウス分布であるとすると，式（１）および雑音がガウス分布であることから，ＩＰ  とＩＴ  の分布はガウス分布Ｎ（ｍＩ  ，σＩ  ２  ＋σＥ  ２  ）で表される。ここにｍＩ  は平均であり，σＩ  ２  ＋σＥ  ２  は分散である。二つの信号ＩＰ  とＩＴ  は式（１）の共通項Ｉ０  によって統計的に依存している。Ｃｐ  の期待値ＣｐＥと非バイアス標準偏差σＣＰはＮデータペアに対して
ＣｐＥ＝σＩ  ２                                      （４）
【００４８】
【数３】
0003554822.tif 000006
【００４９】
である。
【００５０】
一方，バックグラウンドレベルに対しては，
ＣＢ  ＝＜（Ｉｉ  Ｂ  −＜Ｉｉ  Ｂ  ＞）×（Ｉｉ  Ｔ  −＜Ｉｉ  Ｔ  ＞）＞    （６）
ＩＢ  とＩＴ  は独立であるから，ＣＢ  の期待値ＣＢＥと非バイアス標準偏差σＣＢはＮ個のデータペアに対して次のようになる。
【００５１】
ＣＢＥ＝０                                                  （７）
【００５２】
【数４】
0003554822.tif 000007
【００５３】
結果的に，一つのデータ対当たりに規格化されたピークのＳＮ比は，次のようになる。
【００５４】
【数５】
0003554822.tif 000008
【００５５】
ここに，＊  は規格化されたことを表す。
【００５６】
σＥ  ＜σＩ  の時，ＳＮＲＧ  は，光強度の変動が大きくなるにつれ，１／３１／２  で一定になる（後述する）。ＳＮ比を向上させるにはＮを大きくする必要がある。その時，ＳＮＲＧ  は，ほぼＳＮＲＧ  Ｎ１／２  に比例するように修正される。
（２）  排他的二値化法
本発明の最適化排他的二値化法（ＥＢＭ）による相関計算方法に従えば，上記のような一般的方法でのランダムパターンのＳＮＲＧ  の限界を越えることが可能になる。ＥＢＭは，ＴＡＤの画素のデータのうち，ペアになるＳＤの値と等しくないもの（相関のない画素）のデータを除外する（画素値を０にする）。各ＴＡＤのフレームでのＩｉ  ＴＡＤ  （ｘｖ  ）は，ペアのＳＤの記録データであるＩｉ  ＳＤ（τ）と比較され，次のように二値化される。即ち，両者が同じ値の時，画素値を１し，両者の値が異なる時，画素値を０にする。実際は雑音の影響があるので，両者が全く同じであることはないので，本発明では，両者の大きさがある許容誤差範囲δの範囲にある時に両者は等しいとする。そのため，本発明の排他的二値化画像は次のように表される。
【００５７】
Ｉｉ  ＳＤ（τ）−δ≦Ｉｉ  ＴＡＤ  （ｘｖ  ）≦Ｉｉ  ＳＤ（τ）＋δの時
Ｂｉ  （ｘｖ  ）＝１                                          （１０）
ここで，δ＝０として，Ｂｉ  （ｘｖ  ）＝１の代わりに，Ｂｉ  （ｘｖ  ）＝Ａ２  とすると，前述した従来の技術の説明における排他的相関Ｃ’（式２０）と同じになる。
【００５８】
それ以外の時
Ｂｉ  （ｘｖ  ）＝０                                          （１０）’
この方法の場合，許容誤差幅δのとり方でＳＮ比（相関ピーク値とバックグラウンドとの相関値の割合）が変化するので，本発明では許容誤差幅δを，以下に説明するようにランダムパターンの変動の大きさと計測雑音の割合を考慮して最適になるように定める。
【００５９】
二次元の相関強度Ｂ（ｘｖ  ）は排他的二値化画像｛Ｂｉ  （ｘｖ  ）｝の平均で求められる。即ち，
Ｂ（ｘｖ  ）＝＜Ｂｉ  （ｘｖ  ）＞                                （１１）
この式は，一般的方法の式（２）に相当する。相関ピークはＢ（ｘｖ  ）のピークとして現れる。
【００６０】
ＥＢＭでは，相関強度のピークは，強度Ｉｉ  Ｐ  が正しくＩｉ  Ｔ  に等しいとして認識される確率ＰＰ  である。雑音強度を除いて，Ｉｉ  Ｐ  とＩｉ  Ｔ  に含まれるＩｉ  ０  は等しい。誤りは，Ｅｉ  ｐ  とＥｉ  Ｔ  が許容誤差範囲δより大きい時に生じる。従って，この誤差を生じる時の確率は次の式で表される。
【００６１】
【数６】
0003554822.tif 000009
【００６２】
ここに，Ｎ（０，２σＥ  ２  ）は平均０，分散２σＥ  ２  の合成ガウス分布である。
【００６３】
一般的方法と異なり，ＰＰ  はσＩ  の依存性がない。δ〜σＥ  の時，式（１２）の確率密度のほとんど全部をカバーしている。この時，ＰＰ  はほぼ１に等しい。その期待値ＢＰＥと非バイアスの標準偏差σＢＰは次のようになる。
【００６４】
ＢＰＥ＝ＰＰ                                                    （１３）
【００６５】
【数７】
0003554822.tif 000010
【００６６】
一方，非相関信号であるべき強度Ｉｉ  Ｂ  が雑音のためにＩｉ  ０  程度で現れ，相関がある信号として処理されることがある。この誤認識する確率ＰＢ  は，Ｉｉ  Ｂ  ＝（Ｉｉ  Ｘ  ＋Ｅｉ  Ｂ  ）とＩｉ  Ｔ  ＝（Ｉｉ  ０  ＋Ｅｉ  Ｔ  ）の差が許容誤差範囲δの中にある確率と同じである。その時，確率は，Ｉ０  ，ＩＸ  ，ＥＢ  およびＥＴ  の合成確率密度関数を使用して次のように表される。
【００６７】
【数８】
0003554822.tif 000011
【００６８】
ここに，Ｎ（０，２σＩ  ２  ＋２σＥ  ２  ）は，合成ガウス確率密度関数である。平均値は０，偏差は２σＩ  ２  ＋２σＥ  ２  である。δ〜σＥ  ≪σＩ  （即ち，δ≪２σＩ  ２  ＋２σＥ  ２  である）。式（１５）において，ＰＢ  は確率密度関数のＩ＝０付近の小さい部分である。そのためＰＢ  ＝０に近い。Ｎ個のデータペアに対するその期待値ＢＢＥと非バイアス標準偏差σＢＢは次のようになる。
【００６９】
ＢＢＥ＝ＰＢ                                                    （１６）
【００７０】
【数９】
0003554822.tif 000012
【００７１】
式（９）のように，一つのデータ対当たりに規格化されたピークのＳＮ比は次のようになる。
【００７２】
【数１０】
0003554822.tif 000013
【００７３】
Ｎデータペアに対するＳＮ比は，Ｎが増加するにつれＳＮＲＢ  はＳＮＲＢ  ＊  ×Ｎ１／２  に比例する。
【００７４】
ここで重要な点は，ＰＰ  とＰＢ  がランダムパターンの光強度の変動σＩ  に影響されないことである。これらは，主に計測雑音σＥ  にのみ決められる。式（１２）からＰＰ  はσＥ  と独立である。ＰＢ  は光の強度変動に影響されるが（式１５），確率はδにより調整できる。計測雑音σＥ  は光強度の変動σＩ  より一般的に小さい。従って，排他的二値化によってＳＮ比を向上させることが可能である。実際の装置にＥＢＭを適用する前に，許容誤差幅δは，ある与えられたσＩ  とσＥ  に関して最良のＳＮ比になるように最適化しておく。ＳＮ比を検討する場合，あらかじめ概算されたパラメータσＩ  ／σＥ  とδ／σＥ  が利用できる。というのは，システムが改良されてσＥ  が減少しても，δ／σＥ  を一定になるようにすることにより，式（１２）のＰＰ  は同じになるからである。最適化は式（１２）と式（１５）を使用して式（１８）を数値計算し，ＳＮＲＢ  ＊  が大きくなるようにすることにより達成できる。σＩ  ／σＥ  に関してパラメータδ／σＥ  の最適値が図２（ａ）の表の第２欄とに挙げられ，図２（ｂ）にプロットされている。最適値はσＩ  とσＥ  に弱く依存している。３から１００のσＩ  ／σＥ  の現実的な場合に，δ／σＥ  は２．５から４．３に変化する。
【００７５】
（３）一般的方法と排他的二値化法との比較
排他的二値化法と一般的方法とにより計算されたそれぞれのＳＮ比のグラフが図３（ａ）と図２（ａ）の表の第３欄と第４欄に示されている。図３（ａ）のグラフにおいて，ＥＢＭはσＩ  ／σＥ  に対してＳＮ比が最大になるδ／σＥ  を用いたものである。ＥＢＭにより得られる最適ＳＮ比はσＩ  ／σＥ  が増加するにつれて増加するが，一般的な方法のＳＮ比はσＩ  ／σＥ  ＞２で一定（＝１／３１／２  ）である。現実的な装置のσＩ  ／σＥ  は，２から１００程度あるので，ＳＮＲＢ  ＊  は０．６から６．０に増加する。図３（ｂ）は排他的二値化法と一般法とでのＳＮ比の割合を示す。その割合はσＩ  ／σＥ  の平方根にほぼ比例して増加するが，データ数には依存しない。σＩ  ／σＥ  ≧２１／２  のとき，ＳＮＲＢ  ＊  はＳＮＲＧ  より大きく，例えば，ＳＮＲＢ  ／ＳＮＲＧ  ＝３．３と１０．３（σＩ  ／σＥ  ＝１０と１００）である。しかし，σＩ  ／σＥ  ≦２１／２  のとき，ＳＮＲＢ  はＳＮＲＧ  より小さい。このことから分かるように，良いＳＮ比で測定するためには，計測雑音を減少させることがＥＢＭでは重要である。
【００７６】
シミュレーションの例
ＥＢＭの有効性を示すために，前出のモデルで説明された場合についてシミュレーションを行なった。次のパラメータが仮定されている。ＴＡＤ画素は１１×１１（Ｍ＝１１）である。画素の大きさはΔｄ＝１ｍｍである。フレーム間隔はｔｍａｘ  ＝５ｍｓである。ＳＤのサンプリング間隔はΔｔ＝１ｍｓである。各画素の光強度は統計的に独立であり，ガウス分布しているとする（平均ｍＩ  ＝１０，標準偏差σＩ  ＝１，Ｎ（１０，１））。統計的に独立な計測雑音を，ＴＡＤとＳＤの各画素により得られたデータに付加した。その確率密度関数は平均０のガウス関数である。計測雑音の標準偏差に対して２つの場合が考えられる。σＥ  ＝０．１（すなわち，σＩ  ／σＥ  ＝１０）と他の場合はσＥ  ＝０．０１（即ち，σＩ  ／σＥ  ＝１００）である。雑音は光強度と統計的に独立であると仮定されている。ＥＢＭのパラメータδ／σＥ  はσＩ  ／σＥ  ＝１０およびσＩ  ／σＥ  に対して，図２（ａ）の表と図２（ｂ）を参照して，最適値として３．２と４．１をセットする。ランダムパターンは速度ｖｖ  ＝（ｖｘ  ，ｖｙ  ），ｖｘ  ＝−１ｍ／ｓ，ｖｙ  ＝−１ｍ／ｓで移動する。ＴＡＤフレームとＳＤデータの２５ペアが同期的に得られ，以下の相関計算に使用された。
【００７７】
図４（ａ）と（ｂ）は，遅延時間０の状態の排他的二値化法による場合と一般的方法による場合の３次元プロットを示している。図４（ａ）は排他的二値化法の場合であり，図４（ｂ）は一般的方法の場合である。σＩ  ／σＥ  ＝１０の場合である。遅延時間０の二次元の相関をＥＢＭの式（１１）と一般的方法の式（２）を使用して計算する。これらの図は次のように規格化されている。相関のある場合は１に規格化した，また，相関のない期待値は０に規格化した。
【００７８】
ＥＢＭにより決定されたピーク高さのＳＮ比は，一般的方法より良くなっている。ＳＮＲＢ  ＝１１であり，ＳＮＲＧ  ＝２．９である（図２（ａ）の表もしくは図２（ｂ）の評価は，９．０と２．８である）。
【００７９】
図５（ａ），（ｂ）は，それぞれ図４（ａ），（ｂ）の遅延時間τ＝３ｍｓの二次元の相関強度である。そのパターンは３ｍｓ後に中心に移動するように現れている。これらの結果は，ＥＢＭでは二次元の相関においてはバックグラウンドレベルの変動が減少している。そのために，一般的方法より明瞭にピーク相関位置の動きを決定することができることを示している。
【００８０】
図６と図７は，図４と図５と同様のものを示していて，σＩ  ／σＥ  ＝１００の場合である。図６（ａ）は排他的二値化方法の初期状態を示し，図７（ａ）は排他的二値化方法の遅延時間τ＝３ｍｓの状態を示す。図６（ｂ）は一般的方法の初期状態を示し，図７（ｂ）は一般的方法の遅延時間τ＝３ｍｓの状態を示す。相関期待値は，ＢＰＥ＝１，０，ＢＢＥ＝０．０２である（σＩ  ／σＥ  ＝１００）。ＥＢＭの場合にはＳＮ比が向上している。ＳＮＲＢ  ＝２９であり，ＳＮＲＧ  ＝２．８である。図２（ａ）の表と図２（ｂ）のグラフによる評価では２９と２．８である。
【００８１】
図８は，本発明の装置構成の実施の形態を示す図である。
【００８２】
図８において，５はアレイ素子光検出器（ＴＡＤ）である。６は単素子光検出器（ＳＤ）である。２１は入力インタフェースであって，アレイ素子光検出器（ＴＡＤ）５，単素子光検出器（ＳＤ）６の入力インタフェースである。２２はコンピュータである。２３は入力部であって，キーボード等である。２４は出力部であって，ディスプレイ，プリンタ等である。
【００８３】
３１はパラメータ保持部であって，相関演算に必要な各種パラメータを保持するものである。３２は相関演算部であって，アレイ素子光検出器５および単素子光検出器６のデータおよび各種パラメータにより計測信号Ａと参照信号Ｂの相関を求めるものである。相関演算部３２は二値化されたフレーム画像の画素値を画素毎に累積して保持する二値化相関バッファを備える。３３はデータ保持部，アレイ素子光検出器５および単素子光検出器６から入力されるデータを保持するものである。３４は同期制御部であって，アレイ素子光検出器５および単素子光検出器６から入力されるデータの同期をとるものである。３５は最適パラメータ演算部であって，各種パラメータに従って，ＳＮ比が最大にあるようなδ／σＥ  を求めるものである。
【００８４】
図２の構成において，σＩ  ，σＥ  ，最適なδ等の各種パラメータが入力部２３から入力され，パラメータ保持部３１に保持される。σＩ  とσＥ  は予め計測して求めておくか，あるいはＳＤの遅延時間０の記録データを使用して求める。この場合，データ保持部３３の記録データをもとに，最適パラメータ演算部３５で計算させることができる。同期制御部３４はアレイ素子光検出器５，単素子光検出器６から時系列に沿って入力されるデータの同期をフレーム毎にとる。同期をとって入力されたそれぞれのデータはデータ保持部３３に保持される。相関演算部３２はデータ保持部３３に保持されているアレイ素子光検出器５および単素子光検出器６からの入力データと，パラメータ保持部３１に保持されているパラメータに従って，前述の排他的二値化方法により排他的二値化を行ない，参照信号と計測信号の相関演算を行ない，相関ピークを求める。求められた相関ピークの座標，遅延時間等は出力部２４に出力される。
【００８５】
入力部２３から最適パラメータを求めるのに必要な各種パラメータが入力される。σＩ  とσＥ  は，あらかじめ求めておいたものを入力するか，あるいはデータ保持部３３の記録データをもとに最適パラメータ演算部３５で求めるようにしても良い。最適パラメータ演算部３５は，各種パラメータに従って，σＩ  ／σＥ  に対してＳＮ比が最大になるように最適化されたδ／σＥ  を求める。その演算結果は出力部２４に出力される。
【００８６】
図９は，本発明のデータ記録のフローチャートである。図９は図１（ａ），（ｂ）の構成において，ＴＡＤのフレーム画像とＳＤのデータを同期的に取得する動作のフローチャートである。図１（ａ）のアレイ素子光検出器５（ＴＡＤ）とビームスプリッタ４の間にシャッタがあって，アレイ素子光検出器５の露光が終了した後も，単素子光検出器６（ＳＤ）はデータを取得し続けることができる。
【００８７】
Ｓ１でｉ＝１に初期化する。Ｓ２でＴＡＤの露光，ＳＤのデータ取得を同期して開始する。Ｓ３で露光時間待つ（この間，ＴＡＤの画像データを取得するとともにＳＤの一連のデータをとる）Ｓ４で露光時間を終了する。Ｓ５で所定の時間待つ（τの最大値より大きい時間とる）。この間，ＳＤは一連のデータを録りつづける。Ｓ６で一連のＳＤデータの取得を終了する。Ｓ７でＴＡＤのフレーム画像，ＳＤのデータ列を記録して保存する。Ｓ８で予定数の画像フレームのデータが取得できたか判定し，予定数のデータが取得できていれば，処理を終了し，予定数のデータがとれていなければＳ９でｉを更新し，Ｓ２以降の処理を繰り返す。
【００８８】
図１０は本発明の相関演算のフローチャートである。
【００８９】
Ｓ１  二値化相関値バッファの画素値を全て０に初期化する（ｉ＝１）。
【００９０】
Ｓ２でｉフレームのＴＡＤとＳＤのデータ対を入力する。Ｓ３でＳＤの記録データとＴＡＤの各画素値を比較する。Ｓ４とＳ５で画素値が（Ｉｉ  ＳＤ−δ）と（Ｉｉ  ＳＤ＋δ）の間にあるか，ないかを判定する。画素値が（Ｉｉ  ＳＤ−δ）と（Ｉｉ  ＳＤ＋δ）の間にあれば，Ｓ６で画素値を１にする。間になければ，Ｓ７で画素値を０にする。Ｓ８で二値化されたＴＡＤのフレーム画像の画素値を二値化相関値バッファに加算する。Ｓ９でさらにデータ対を追加するか判定し，さらに追加するのであれば，Ｓ１０でｉ＝ｉ＋１とし，Ｓ２以降の処理を繰り返す。なお，本発明は，取得したデータを保存してから相関計算する方法だけでなく，データを取得しながら相関計算をすることも可能である。従って，Ｓ９における追加は，データ保持部３３に記録されたデータ対を読み出すか，あるいは新たなデータ対を取得することを意味する。Ｓ９の判断でさらに追加しないのであれば，Ｓ１１で二値化相関値バッファの各画素値をｉで割り，平均をとる。Ｓ１２で相関結果を出力して，処理を終了する。
【００９１】
図１１は本発明の最適パラメータ決定のフローチャートである。
【００９２】
Ｓ１でσＩ  とσＥ  を定める。Ｓ２で，δ／σＥ  の値を変えながら各δ／σＥ  に対してＰＢ  （式１２）とＰＢ  （式１５）を数値計算する。Ｓ３でＳＮ比（式１８）を計算し，σＩ  ／σＥ  毎にＳＮ比を最大するδ／σＥ  を選ぶ。
【００９３】
Ｓ４で最適パラメータを出力する。
【００９４】
【発明の効果】
本発明によれば，参照信号と計測信号の大きさの違いに基づく二値化を，最適化された誤差許容幅に従って，相関ピークの鮮明さを示すＳＮ比が最大になるように行なうことができる。そのため，高い精度で，相関ピークを決定することができる。そのため，高速度で移動する物体の動き等を高速で正確に測定できる。そのため，地球大気のゆらぎの観測等に利用することも可能で，地球大気のゆらぎ層の風速，風向を正確に観測することが可能になる。
【図面の簡単な説明】
【図１】本発明の実施の形態１を示す図である。
【図２】本発明の結果を示すグラフ１を示す図である。
【図３】本発明の結果を示すグラフ２を示す図である。
【図４】本発明の相関結果の三次元図１である。
【図５】本発明の相関結果の三次元図２である。
【図６】本発明の相関結果の三次元図３である。
【図７】本発明の相関結果の三次元図４である。
【図８】本発明の実施の形態２を示す図である。
【図９】本発明のデータ記録のフローチャートを示す図である。
【図１０】本発明の相関演算のフローチャートを示す図である。
【図１１】本発明の最適パラメータ決定のフローチャートを示す図である。
【図１２】従来の技術の説明図である。
【符号の説明】
１：ランダムパターン
３：集光レンズ
４：ビームスプリッタ
５：アレイ素子光検出器
６：単素子光検出器
７：ピンホール
８：計測光
９：参照光
【図１】
0003554822.tif 000014
【図２】
0003554822.tif 000015
【図３】
0003554822.tif 000016
【図４】
0003554822.tif 000017
【図５】
0003554822.tif 000018
【図６】
0003554822.tif 000019
【図７】
0003554822.tif 000020
【図８】
0003554822.tif 000021
【図９】
0003554822.tif 000022
【図１０】
0003554822.tif 000023
【図１１】
0003554822.tif 000024
【図１２】
0003554822.tif 000025