JP2022042885A - 測定装置および測定方法 - Google Patents

Abstract

【課題】光学式距離計といった測定装置において、測定精度の低下を抑制しつつ測定速度を向上させる測定装置を提供する。【解決手段】周波数変調レーザ光を出力するレーザ装置と、周波数変調レーザ光の一部の参照光と、残りの少なくとも一部の測定光とに分岐させる分岐部と、測定光を計測対象物に照射して反射された反射光と参照光とを混合してビート信号を発生させるビート信号発生部１５０と、ビート信号の時系列データを予め定められた多項式に当てはめた場合のモデルパラメータを推定する推定部３２０と、モデルパラメータに基づいて、多項式に対応する特性方程式が極値となる解を算出する第１算出部３３０と、算出された解のうち、絶対値が最大となる解を探索し、探索した解に対応するビート周波数を算出する探索部３４０と、ビート周波数に基づき、参照光と測定光との伝搬距離の差を算出する第２算出部３５０とを備える、測定装置。【選択図】図４

Description

本発明は、測定装置および測定方法に関する。

共振器内に周波数シフタが設けられ、時間の経過とともに発振周波数が線形に変化する複数の縦モードレーザを出力する周波数シフト帰還レーザ（ＦＳＦＬ：Frequency Shifted Feedback Laser）が知られている。また、このような周波数シフト帰還レーザを用いた光学式の距離計が知られている（例えば、特許文献１および非特許文献１を参照）。

特許第３５８３９０６号明細書

原武文，「ＦＳＬレーザによる距離センシングとその応用」，オプトニューズ，Ｖｏｌ．７，Ｎｏ．３，２０１２年，ｐｐ．２５－３１ 足立修一，「システム同定の基礎」，東京電機大学出版局，ｐｐ．１７０－１８０

周波数シフト帰還レーザを用いた光学式距離計は、非接触で大量の三次元情報を取得可能であり、例えば、設計および生産現場等で用いられてきた。このような光学式距離計は、スループットの低減を抑制しつつ、より高精度に測定できることが望まれていた。

そこで、本発明はこれらの点に鑑みてなされたものであり、光学式距離計といった測定装置において、測定精度の低下を抑制しつつ測定速度を向上させるようにすることを目的とする。

本発明の第１の態様においては、周波数変調レーザ光を出力するレーザ装置と、前記レーザ装置が出力する前記周波数変調レーザ光を、前記周波数変調レーザ光の一部の参照光と、残りの少なくとも一部の測定光とに分岐させる分岐部と、前記測定光を計測対象物に照射して反射された反射光と前記参照光とを混合してビート信号を発生させるビート信号発生部と、前記ビート信号の時系列データを予め定められた多項式に当てはめた場合のモデルパラメータを推定する推定部と、前記モデルパラメータに基づいて、前記多項式に対応する特性方程式が極値となる解を算出する第１算出部と、算出された解のうち、絶対値が最大となる解を探索し、探索した解に対応するビート周波数を算出する探索部と、前記ビート周波数に基づき、前記参照光と前記測定光との伝搬距離の差を算出する第２算出部とを備える、測定装置を提供する。

前記推定部は、前記ビート信号の前記時系列データを自己回帰モデルで表現した場合の前記モデルパラメータを推定してもよい。

前記推定部は、最小２乗法、ユールウォーカ法、バーグ法、および共分散法のうちいずれか１つの方法に基づき、前記モデルパラメータを推定してもよい。

前記推定部は、前記ビート信号の前記時系列データから予め定められた漸化式を用いて前記モデルパラメータを推定してもよい。

前記推定部は、逐次最小２乗法に基づき、前記モデルパラメータを推定してもよい。

前記第１算出部は、前記多項式にＺ変換を適用したことに基づいて算出される伝達関数に対応する前記特性方程式を用いてもよい。

前記第１算出部は、前記伝達関数に基づく前記特性方程式を因数分解することにより、前記特性方程式が極値となる解を算出してもよい。

前記探索部は、探索した解の偏角に基づき、探索した解に対応する前記ビート周波数を算出してもよい。

本発明の第２の態様においては、レーザ装置から複数のモードの周波数変調レーザ光を出力するステップと、前記レーザ装置が出力する前記周波数変調レーザ光を、前記周波数変調レーザ光の一部の参照光と、残りの少なくとも一部の測定光とに分岐させるステップと、前記測定光を計測対象物に照射して反射された反射光と、前記参照光とを混合してビート信号を発生させるステップと、前記ビート信号の時系列データを予め定められた多項式に当てはめた場合のモデルパラメータを推定するステップと、前記モデルパラメータに基づいて、前記多項式に対応する特性方程式が極値となる解を算出するステップと、算出された解のうち、絶対値が最大となる解を探索し、探索した解に対応するビート周波数を算出するステップと、前記ビート周波数に基づき、前記参照光と前記測定光との伝搬距離の差を算出するステップとを備える、測定方法を提供する。

本発明によれば、光学式距離計といった測定装置において、測定精度の低下を抑制しつつ測定速度を向上できるという効果を奏する。

本実施形態に係る測定装置１００の構成例を計測対象物１０と共に示す。 本実施形態に係るレーザ装置１１０の構成例を示す。 本実施形態に係るレーザ装置１１０が出力するレーザ光の一例を示す。 本実施形態に係るビート信号発生部１５０および周波数解析部１６０の構成例を示す。

［測定装置１００の構成例］
図１は、本実施形態に係る測定装置１００の構成例を計測対象物１０と共に示す。測定装置１００は、当該測定装置１００および計測対象物１０の間の距離を光学的に測定する。また、測定装置１００は、計測対象物１０に照射するレーザ光の位置を走査して、計測対象物１０の三次元的な形状を計測してもよい。測定装置１００は、レーザ装置１１０と、分岐部１２０と、光サーキュレータ１３０と、光ヘッド部１４０と、ビート信号発生部１５０と、周波数解析部１６０と、表示部１７０を備える。

レーザ装置１１０は、レーザ共振器を有し、複数のモードの周波数変調レーザ光を出力する。レーザ装置１１０は、共振器内に周波数シフタが設けられ、時間の経過とともに発振周波数が線形に変化する複数の縦モードレーザを出力する。レーザ装置１１０は、一例として、周波数シフト帰還レーザである。周波数シフト帰還レーザについては後述する。

分岐部１２０は、レーザ装置１１０が出力する周波数変調レーザ光を、周波数変調レーザ光の一部の参照光と、残りの少なくとも一部の測定光とに分岐させる。分岐部１２０は、一例として、光ファイバ型の光カプラである。図１の例において、分岐部１２０は、測定光を光サーキュレータ１３０に供給し、参照光をビート信号発生部１５０に供給する。

光サーキュレータ１３０は、複数の入出力ポートを有する。光サーキュレータ１３０は、例えば、一のポートに入力した光を次のポートから出力させ、当該次のポートから入力する光を更に次のポートから出力させる。図１は、光サーキュレータ１３０が３つの入出力ポートを有する例を示す。この場合、光サーキュレータ１３０は、分岐部１２０から供給される測定光を光ヘッド部１４０に出力する。また、光サーキュレータ１３０は、光ヘッド部１４０から入力する光をビート信号発生部１５０へと出力する。

光ヘッド部１４０は、光サーキュレータ１３０から入力する光を計測対象物１０に向けて照射する。光ヘッド部１４０は、一例として、コリメータレンズを有する。この場合、光ヘッド部１４０は、光ファイバを介して光サーキュレータ１３０から入力する光をコリメータレンズでビーム状に調節してから出力する。

また、光ヘッド部１４０は、計測対象物１０に照射した測定光の反射光を受光する。光ヘッド部１４０は、受光した反射光をコリメータレンズが光ファイバに集光して光サーキュレータ１３０に供給する。この場合、光ヘッド部１４０は、共通の１つのコリメータレンズを有し、当該コリメータレンズが、測定光を計測対象物１０に照射し、また、計測対象物１０からの反射光を受光してよい。なお、光ヘッド部１４０および計測対象物１０の間の距離をｄとする。

これに代えて、光ヘッド部１４０は、集光レンズを有してもよい。この場合、光ヘッド部１４０は、光ファイバを介して光サーキュレータ１３０から入力する光を計測対象物１０の表面に集光する。そして、光ヘッド部１４０は、計測対象物１０の表面で反射した反射光の少なくとも一部を受光する。光ヘッド部１４０は、受光した反射光を集光レンズが光ファイバに集光して光サーキュレータ１３０に供給する。この場合においても、光ヘッド部１４０は、共通の１つの集光レンズを有し、当該集光レンズが、測定光を計測対象物１０に照射し、また、計測対象物１０からの反射光を受光してよい。

ビート信号発生部１５０は、測定光を計測対象物１０に照射して反射された反射光を光サーキュレータ１３０から受けとる。また、ビート信号発生部１５０は、分岐部１２０から参照光を受けとる。ビート信号発生部１５０は、反射光および参照光を混合してビート信号を発生させる。ビート信号発生部１５０は、例えば、光電変換素子を有し、ビート信号を電気信号に変換して出力する。

ここで、反射光は、光ヘッド部１４０から計測対象物１０までの距離を往復しているので、参照光と比較して少なくとも距離２ｄに応じた伝搬距離の差が生じることになる。レーザ装置１１０が出力する光は、時間の経過とともに発振周波数が線形に変化するので、参照光および反射光の発振周波数は、当該伝搬距離の差に対応する伝搬遅延に応じた周波数差が生じる。ビート信号発生部１５０は、このような周波数差に対応するビート信号を発生させる。

周波数解析部１６０は、ビート信号発生部１５０が発生させたビート信号を周波数解析して、参照光と測定光との伝搬距離の差を検出する。周波数解析部１６０は、予め定められたサンプリング周波数でビート信号をサンプリングして、ビート信号を周波数解析する。周波数解析部１６０の周波数解析については後述する。

表示部１７０は、周波数解析部１６０の解析結果を表示する。表示部１７０は、ディスプレイ等を有し、検出結果を表示してよい。また、表示部１７０は、記憶部等に解析結果を記憶させてもよい。

以上の測定装置１００は、計測対象物１０に照射した測定光の反射光と、参照光との間の周波数差を解析することにより、測定装置１００および計測対象物１０の間の距離を測定可能とする。即ち、測定装置１００は、非接触および非破壊の光学式距離計を構成できる。測定装置１００のより詳細な構成について次に説明する。

［レーザ装置１１０の構成例］
図２は、本実施形態に係るレーザ装置１１０の構成例を示す。図２のレーザ装置１１０は、周波数シフト帰還レーザの一例を示す。レーザ装置１１０は、レーザ共振器を有し、当該レーザ共振器内でレーザ光を発振させる。レーザ装置１１０のレーザ共振器は、周波数シフタ１１２と、増幅媒体１１４と、ＷＤＭカプラ１１６と、ポンプ光源１１７と、出力カプラ１１８とを含むレーザ共振器を有する。

周波数シフタ１１２は、入力する光の周波数を略一定の周波数だけシフトする。周波数シフタ１１２は、一例として、音響光学素子を有するＡＯＦＳ（Acousto-Optic Frequency Shifter）である。ここで、周波数シフタ１１２による周波数シフト量を＋ν_ｓとする。即ち、周波数シフタ１１２は、共振器を周回する光の周波数を、１周回毎にν_ｓだけ周波数が増加するようにシフトさせる。

増幅媒体１１４は、ポンプ光が供給され、入力光を増幅する。増幅媒体１１４は、一例として、不純物が添加された光ファイバである。不純物は、例えば、エルビウム、ネオジウム、イッテルビウム、テルビウム、ツリウム等の希土類元素である。また、増幅媒体１１４は、ＷＤＭカプラ１１６を介してポンプ光源１１７からポンプ光が供給される。出力カプラ１１８は、共振器内でレーザ発振した光の一部を外部に出力する。

即ち、図２に示すレーザ装置１１０は、共振器内に周波数シフタ１１２を有するファイバリングレーザを構成する。レーザ装置１１０は、共振器内にアイソレータを更に有することが望ましい。また、レーザ装置１１０は、予め定められた波長帯域の光を通過させる光バンドパスフィルタを共振器内に有してもよい。このようなレーザ装置１１０が出力するレーザ光の周波数特性について次に説明する。

図３は、本実施形態に係るレーザ装置１１０が出力するレーザ光の一例を示す。図３は、時刻ｔ_０においてレーザ装置１１０が出力するレーザ光の光スペクトルを左側に示す。当該光スペクトルにおいては、横軸が光強度、縦軸が光の周波数を示す。また、光スペクトルの複数の縦モードを番号ｑで示す。複数の縦モードの周波数は、略一定の周波数間隔で並ぶ。ここで、光が共振器を１周する時間をτ_ＲＴ（＝１／ν_Ｃ）とすると、複数の縦モードは、次式のように１／τ_ＲＴ（＝ν_Ｃ）間隔で並ぶことになる。なお、ν_０は、時刻ｔ_０における光スペクトルの初期周波数とする。また、ν_Ｃは、レーザ共振器の共振周波数νｃである。

図３は、レーザ装置１１０が出力する複数の縦モードの時間経過にともなう周波数の変化を右側に示す。図３の右側においては、横軸が時間、縦軸が周波数を示す。即ち、図３は、レーザ装置１１０が出力するレーザ光の周波数の時間的な変化を右側に示し、当該レーザ光の時刻ｔ_０における瞬時周波数を左側に示したものである。

レーザ装置１１０は、共振器内の光が共振器を１周する毎に、周波数シフタ１１２が周回する光の周波数をν_ｓだけ増加させる。即ち、時間がτ_ＲＴ経過する毎に、各モードの周波数はν_ｓだけ増加するので、周波数の時間変化ｄν／ｄｔは、ν_ｓ／τ_ＲＴと略等しくなる。したがって、（数１）式で示した複数の縦モードは、時間ｔの経過に伴って、次式のように変化する。

［距離測定処理の詳細］
本実施形態に係る測定装置１００は、（数２）式で示すような周波数成分を出力するレーザ装置１１０を用いて、例えば、光ヘッド部１４０および計測対象物１０の間の距離ｄを測定する。ここで、参照光および反射光の間の光路差が、距離ｄを往復した距離２ｄだけであり、距離２ｄに対応する伝搬遅延をΔｔとする。即ち、時刻ｔにおいて、測定光が計測対象物１０から反射して戻ってきた場合、戻ってきた反射光の周波数は、時刻ｔよりも時間Δｔだけ過去の周波数と略一致するので、次式で示すことができる。

一方、時刻ｔにおける参照光は、（数２）式と同様に次式で示すことができる。ここで、参照光をν_ｑ’（ｔ）とした。

ビート信号発生部１５０は、このような反射光および参照光を重畳させるので、（数３）式の複数の縦モードと（数４）式で示す複数の縦モードとの間の複数のビート信号が発生することになる。このようなビート信号の周波数をν_Ｂ（Ｍ，ｄ）とすると、ν_Ｂ（Ｍ，ｄ）は、（数３）式および（数４）式より次式で示すことができる。なお、Ｍを縦モード番号の間隔（＝ｑ－ｑ’）とし、Δｔ＝２ｄ／ｃとした。

（数５）式より、距離ｄは、次式のように示される。ここで、１／τ_ＲＴ＝ν_ｃとした。

（数６）式より、縦モード番号の間隔Ｍを判別すれば、ビート信号の周波数観測結果から距離ｄを算出できることがわかる。なお、間隔Ｍは、レーザ装置１１０の周波数シフト量ν_ｓを変化させた場合のビート信号の変化を検出することで、判別することができる。このような間隔Ｍの判別方法は、特許文献１等に記載されているように既知であるから、ここでは詳細な説明を省略する。

（数６）式より、ビート信号の信号周波数ν_Ｂ（Ｍ，ｄ）の測定精度が距離ｄの測定精度に与える影響は、ν_ｓおよびν_ｃの値が大きいほど小さくなる。また、ν_ｃにも測定誤差が含まれる場合があるので、Ｍの値は小さい方が好ましい。なお、距離ｄが同じであれば、ν_ｓが小さいほどＭの値は小さくなる傾向にある。したがって、ν_ｓの値を大きくすることと比較して、ν_ｃの値を大きくした方が距離ｄの測定精度を向上させるためには好ましい。しかしながら、ν_ｃは観測帯域においてビート信号が最大となる周波数なので、ν_ｃの値を大きくすることは、サンプリング周波数を大きくすることに相当する。

また、ν_Ｂ（Ｍ，ｄ）およびν_ｃの測定精度は、サンプリング時間の逆数である周波数分解能に大きく影響を受ける。即ち、ν_Ｂ（Ｍ，ｄ）およびν_ｃの測定精度を向上させるためには、サンプリング時間を長くすることが好ましい。

このように、距離ｄの測定精度の向上には、より高いサンプリング周波数で、より長いサンプリング時間を用いてビート信号をサンプリングすることが好ましいので、サンプルデータ数が増加してしまうことになる。サンプルデータ数が増加すると、周波数解析にかかる時間が長くなるので、ビート信号による距離測定は、スループットが低下してしまうことがある。

ここで、周波数解析は、従来、デジタルフーリエ変換を用いていた。デジタルフーリエ変換に基づく周波数解析は、例えば、時系列に並ぶデータ列の先頭のデータと最後のデータとを連続したデータとするために窓関数を各データに乗じる窓関数処理、周波数に変換する変換処理、変換した周波数データからピーク周波数を探索する探索処理等を含む。サンプルしたデータの全てにこのような処理を実行すると、処理時間がかかってしまい、測定精度の低下を抑制しつつ測定速度を向上させることは困難であった。

そこで、本実施形態に係る測定装置１００は、時系列データをモデルパラメータで表現し、当該モデルパラメータに対応する伝達関数を用いることにより、周波数解析を高速に実行する。このような周波数解析を実行する測定装置１００の周波数解析部１６０について、次に説明する。

［ビート信号発生部１５０および周波数解析部１６０の構成例］
図４は、本実施形態に係るビート信号発生部１５０および周波数解析部１６０の構成例を示す。ビート信号発生部１５０は、反射光および参照光のビート信号ν_Ｂ（Ｍ，ｄ）を発生させて電気信号に変換する。ビート信号発生部１５０は、光カプラ１５２と光電変換部１５４とを有する。光カプラ１５２は、反射光と参照光とを合波してビート信号ν_Ｂ（Ｍ，ｄ）を発生させる。光電変換部１５４は、合波した反射光および参照光を受光して電気信号に変換する。光電変換部１５４は、フォトダイオード等でよい。光電変換部１５４は、一例として、バランス型フォトダイオードである。

周波数解析部１６０は、ビート信号発生部１５０からビート信号ν_Ｂ（Ｍ，ｄ）を受けとって、周波数解析する。周波数解析部１６０は、フィルタ部１６２、ＡＤ変換部２０２、クロック信号供給部２１０、および信号処理部２２０を有する。

フィルタ部１６２は、ユーザ等が周波数解析したい周波数帯域とは異なる周波数帯域の信号成分を低減させる。ここで、ユーザ等が周波数解析したい周波数帯域を０からν_ｃとする。フィルタ部１６２は、例えば、周波数ν_ｃ以下の信号成分を通過させるローパスフィルタである。この場合、フィルタ部１６２は、周波数ν_ｃよりも高い周波数の信号成分を低減させたビート信号をＡＤ変換部２０２に供給する。

ＡＤ変換部２０２は、入力するアナログ信号をデジタル信号に変換する。ＡＤ変換部２０２はビート信号をデジタル信号に変換する。クロック信号供給部２１０は、ＡＤ変換部２０２にクロック信号を供給する。これにより、ＡＤ変換部２０２は、例えば、受け取ったクロック信号の周波数と略同一のサンプリングレートでアナログ信号をデジタル信号に変換する。

ここで、観測帯域を０からν_ｃとすると、ビート信号の周波数は、最大でもレーザ共振器の共振周波数ν_ｃである。したがって、クロック信号供給部２１０が、レーザ共振器の共振周波数ν_ｃの２倍以上の周波数のクロック信号をＡＤ変換部２０２に供給することで、ビート信号を観測することができる。このように、周波数解析部１６０は、１つのビート信号をサンプリングして生成したビート信号をサンプリングデータとして周波数解析する。

信号処理部２２０は、ビート信号を周波数データに変換する。信号処理部２２０は、例えば、集積回路等で構成されている。この場合、信号処理部２２０は、ＦＰＧＡ（Field Programmable Gate Array）、ＤＳＰ（Digital Signal Processor）、および／またはＣＰＵ（Central Processing Unit）を含む。信号処理部２２０は、記憶部３１０と、推定部３２０と、第１算出部３３０と、探索部３４０と、第２算出部３５０とを備える。

記憶部３１０は、ＡＤ変換部２０２が出力するビート信号の時系列データを記憶する。また、記憶部３１０は、測定装置１００が動作の過程で生成する（または利用する）中間データ、算出結果、閾値、およびパラメータ等をそれぞれ記憶してもよい。また、記憶部３１０は、測定装置１００内の各部の要求に応じて、記憶したデータを要求元に供給してもよい。

記憶部３１０は、例えば、コンピュータ等が信号処理部２２０として機能する場合、コンピュータ等を機能させるＯＳ（Operating System）、およびプログラム等の情報を格納してもよい。また、記憶部３１０は、当該プログラムの実行時に参照されるデータベースを含む種々の情報を格納してもよい。例えば、コンピュータは、記憶部３１０に記憶されたプログラムを実行することによって、記憶部３１０、推定部３２０、第１算出部３３０、探索部３４０、および第２算出部３５０の少なくとも一部として機能する。

記憶部３１０は、例えば、コンピュータ等のＢＩＯＳ（Basic Input Output System）等を格納するＲＯＭ（Read Only Memory）、および作業領域となるＲＡＭ（Random Access Memory）を含む。また、記憶部３１０は、ＨＤＤ（Hard Disk Drive）および／またはＳＳＤ（Solid State Drive）等の大容量記憶装置を含んでもよい。また、コンピュータは、ＧＰＵ（Graphics Processing Unit）等を更に備えてもよい。

［周波数解析の詳細］
推定部３２０は、ビート信号の時系列データを予め定められた多項式に当てはめた場合のモデルパラメータを推定する。推定部３２０は、例えば、ビート信号の時系列データに基づき、時系列データを自己回帰モデルで表現した場合のモデルパラメータを推定する。自己回帰モデルは、例えば、ある時刻ｔの値を、時刻ｔよりも過去のデータを使って表現するモデルである。自己回帰モデルは、一例として、次式のように示すことができる。

ここで、ｙ_ｎは、時系列データにおける時点ｎのデータ値であり、ｙ_ｎ－ｉは、時系列データにおける時点ｎよりも過去の方向のｉ番目のデータ値である。また、ａ_ｉは、回帰係数であり、ν_ｎは、白色雑音である。ｋは、時系列データにおいて過去のどこまでの値を用いて時点ｎのデータ値を算出するかを決定する値であり、次数と呼ぶ。次数ｋは、予め定められた値でよく、また、次数ｋは、自己回帰モデルによる推定結果に基づき、調整されてもよい。

時系列データに対応するａ_ｉ、ν_ｎ、ｋ等の自己回帰モデルのモデルパラメータは、既知のアルゴリズム等によって算出することができる。推定部３２０は、例えば、最小２乗法、ユールウォーカ法、バーグ法、および共分散法のうちいずれか１つの方法に基づき、モデルパラメータを推定する。

第１算出部３３０は、モデルパラメータに基づいて、（数７）式に示したような多項式に対応する特性方程式が極値となる解を算出する。第１算出部３３０が用いる特性方程式を算出するために、自己回帰モデルにＺ変換を適用すると、次式が得られる。

（数８）式に基づき、伝達関数Ｈ（ｚ）を次式のように定義する。

（数９）式より、伝達関数Ｈ（ｚ）の分母を０とした特性方程式を次式のように生成する。

第１算出部３３０は、以上のように、自己回帰モデルの多項式にＺ変換を適用したことに基づいて算出される伝達関数に対応する特性方程式を用いる。そして、第１算出部３３０は、伝達関数に基づく特性方程式を因数分解することにより、特性方程式が極値となる解を算出する。（数１０）式の特性方程式を因数分解することにより、次式を得る。

第１算出部３３０は、（数１１）式の解となる（ｋ＋１）個のα_ｊを算出する。第１算出部３３０が算出する解α_ｊの絶対値は、時系列データを周波数軸に変換した場合の周波数スペクトルのピークレベルに相当する。また、解α_ｊの偏角∠α_ｊは、周波数スペクトルのピークに対応する周波数に相当する。

探索部３４０は、第１算出部３３０によって算出された解α_ｊのうち、絶対値｜α_ｊ｜が最大となる解｜α_ｊ｜_ｍａｘを探索する。そして、探索部３４０は、探索した解｜α_ｊ｜_ｍａｘに対応する周波数をビート周波数ｆ_{ｊ＿ｍａｘ}として算出する。探索部３４０は、探索した解｜α_ｊ｜_ｍａｘの偏角∠｜α_ｊ｜_ｍａｘに基づき、探索した解｜α_ｊ｜_ｍａｘに対応するビート周波数を算出する。なお、サンプリング周波数をＦ_Ｓとした場合の、ビート周波数ｆ_{ｊ＿ｍａｘ}は、次式のように算出できる。

第２算出部３５０は、ビート周波数ｆ_{ｊ＿ｍａｘ}に基づき、参照光と測定光との伝搬距離の差を算出する。第２算出部３５０は、例えば、ビート周波数ｆ_{ｊ＿ｍａｘ}を、（数６）式のν_Ｂ（Ｍ，ｄ）とすることで、光ヘッド部１４０および計測対象物１０の間の距離ｄを算出する。

以上のように、本実施形態に係る周波数解析部１６０は、まず、ビート信号の時系列データを自己回帰モデルで表現した場合のモデルパラメータを推定する。ビート信号は、一定の周波数信号を有するビート信号なので、特定の周波数帯域に変換するＦＦＴと比較して、より少ない時系列データで自己回帰モデルのモデルパラメータを算出することができる。例えば、略同一の距離分解能の距離ｄを算出する場合、第１算出部３３０は、ＦＦＴと比較して、半分程度以下の数の時系列データでモデルパラメータを算出できる。

また、ＦＦＴによる処理を実行する場合、時系列データの数を２^ｎ個にしなければならず、距離分解能等との関係を保つためには設計自由度が制限される。これに対して、モデルパラメータを算出するための時系列データの数には、ＦＦＴのような制限はないので、距離分解能等に対応した時系列データの数を容易に設定できる。

そして、周波数解析部１６０は、自己回帰モデルに基づく特性方程式が極値となる解を算出し、そして、算出された解のうち、絶対値が最大となる解を探索する。ここで、ＦＦＴによる処理は、時間軸のデータを周波数軸のデータに変換するだけなので、変換した周波数軸のデータの全ての中から信号強度がピークレベルとなるデータを探索しなければならない。これに対し、特性方程式を用いた場合、次数ｋに対応するｋ個の解の中から、絶対値が最大となる解を探索するので、例えば、１／１００程度以下から、１／１０００程度以下に処理時間を短縮することができる。

そして、周波数解析部１６０は、探索した解に対応するビート周波数を算出し、当該ビート周波数に基づき、参照光と測定光との伝搬距離の差を算出する。ビート周波数から距離ｄを算出する処理は、既知の処理と同様に実行することができる。したがって、本実施形態に係る測定装置１００によれば、測定精度の低下を抑制しつつ測定速度を向上させることができる。

以上の本実施形態に係る測定装置１００において、推定部３２０が既知のアルゴリズムを用いてビート信号の時系列データを多項式に当てはめた場合のモデルパラメータを推定する例を説明したが、これに限定されることはない。推定部３２０は、例えば、ビート信号の時系列データから予め定められた漸化式を用いてモデルパラメータを推定してもよい。このような推定部３２０の動作について、次に説明する。

推定部３２０は、まず、忘却係数λを設定する。忘却係数λは、過去のデータを複数回計算に用いることを防止するための係数であり、１以下の正の数である。忘却係数λは、例えば、０．９７以上０．９９５以下の値が用いられる。一例として、忘却係数λは次式のように設定される。ここで、τは、推定部３２０が計算に用いる時系列データの最大のサンプル数である。

また、推定部３２０は、共分散行列Ｐ（Ｎ）の初期値Ｐ（０）を単位行列Ｉとする。

ここで、モデルパラメータを要素とするベクトルａとし、推定部３２０は、モデルパラメータの初期値ａ（０）の全ての要素を０とする。なお、本実施例において、算出すべきモデルパラメータの数をｍ個とする。

推定部３２０は、ビート信号の時系列データのうちｍ個のデータを用いてモデルパラメータを推定する。ここで、ｋ番目の時系列データをｙ（ｋ）とし、推定部３２０は、時系列データｙ（ｋ－１）から時系列データｙ（ｋ－ｍ）までの、過去のｍ個のデータを要素とするベクトルｚ（ｋ）を生成する。

ここで、ベクトルｚ（ｋ）を用いて共分散行列Ｐ（Ｎ）を次式のように示す。

このような共分散行列Ｐ（Ｎ）は、時点Ｎのベクトルｚ（Ｎ）と時点Ｎよりも１サンプル前の時点Ｎ－１の共分散行列Ｐ（Ｎ－１）とに基づく、次式のような漸化式を用いて算出できる。

ここで、時点Ｎの時系列データｙ（N）の推定値をｙ’（N）とすると、ｙ’（N）は次式のように示される。

また、ｙ（N）の推定残差をε（Ｎ）とすると、次式のように示される。

このような推定残差ε（Ｎ）を用いることにより、時点Ｎのモデルパラメータの推定値ａ’（N）は、次式のように算出することができる。

推定部３２０は、（数２１）式に示すような漸化式を用いることにより、（数７）式に示す時系列データを自己回帰モデルで表現した多項式のモデルパラメータを算出することができる。なお、このような漸化式によるモデルパラメータの推定は、逐次最小２乗法に基づくものであり、例えば、非特許文献２等に詳細な内容が記載されている。

以上のように、推定部３２０は漸化式を用いてもモデルパラメータを算出できるので、第１算出部３３０は、上述のように、モデルパラメータに基づいて、多項式に対応する特性方程式が極値となる解を算出できる。なお、推定部３２０が逐次最小２乗法に基づいてモデルパラメータを推定する場合、時系列データのデータ毎にモデルパラメータを算出することができる。この場合、推定部３２０は、一の時点の時系列データに基づいて算出したモデルパラメータに、次の時点の時系列データに基づいて算出したモデルパラメータの差分を追加してモデルパラメータを更新する。したがって、推定部３２０はモデルパラメータの差分を順次計算するだけなので、計算量を低減させてより高速にモデルパラメータを算出できる。

この場合、第１算出部３３０、探索部３４０、および第２算出部３５０は、推定部３２０と同様に、時系列データのデータ毎に動作することができる。これにより、測定装置１００は、当該測定装置１００から計測対象物１０までの距離ｄを時系列データのデータ毎に出力することができる。

以上、本発明を実施の形態を用いて説明したが、本発明の技術的範囲は上記実施の形態に記載の範囲には限定されず、その要旨の範囲内で種々の変形及び変更が可能である。例えば、装置の全部又は一部は、任意の単位で機能的又は物理的に分散・統合して構成することができる。また、複数の実施の形態の任意の組み合わせによって生じる新たな実施の形態も、本発明の実施の形態に含まれる。組み合わせによって生じる新たな実施の形態の効果は、もとの実施の形態の効果を併せ持つ。

１０ 計測対象物
１００ 測定装置
１１０ レーザ装置
１１２ 周波数シフタ
１１４ 増幅媒体
１１６ ＷＤＭカプラ
１１７ ポンプ光源
１１８ 出力カプラ
１２０ 分岐部
１３０ 光サーキュレータ
１４０ 光ヘッド部
１５０ ビート信号発生部
１５２ 光カプラ
１５４ 光電変換部
１６０ 周波数解析部
１６２ フィルタ部
１７０ 表示部
２０２ ＡＤ変換部
２１０ クロック信号供給部
２２０ 信号処理部
３１０ 記憶部
３２０ 推定部
３３０ 第１算出部
３４０ 探索部
３５０ 第２算出部

Claims (9)

1. 周波数変調レーザ光を出力するレーザ装置と、
   前記レーザ装置が出力する前記周波数変調レーザ光を、前記周波数変調レーザ光の一部の参照光と、残りの少なくとも一部の測定光とに分岐させる分岐部と、
   前記測定光を計測対象物に照射して反射された反射光と前記参照光とを混合してビート信号を発生させるビート信号発生部と、
   前記ビート信号の時系列データを予め定められた多項式に当てはめた場合のモデルパラメータを推定する推定部と、
   前記モデルパラメータに基づいて、前記多項式に対応する特性方程式が極値となる解を算出する第１算出部と、
   算出された解のうち、絶対値が最大となる解を探索し、探索した解に対応するビート周波数を算出する探索部と、
   前記ビート周波数に基づき、前記参照光と前記測定光との伝搬距離の差を算出する第２算出部と
   を備える、測定装置。
2. 前記推定部は、前記ビート信号の前記時系列データを自己回帰モデルで表現した場合の前記モデルパラメータを推定する、請求項１に記載の測定装置。
3. 前記推定部は、最小２乗法、ユールウォーカ法、バーグ法、および共分散法のうちいずれか１つの方法に基づき、前記モデルパラメータを推定する、請求項２に記載の測定装置。
4. 前記推定部は、前記ビート信号の前記時系列データから予め定められた漸化式を用いて前記モデルパラメータを推定する、請求項１または２に記載の測定装置。
5. 前記推定部は、逐次最小２乗法に基づき、前記モデルパラメータを推定する、請求項４に記載の測定装置。
6. 前記第１算出部は、前記多項式にＺ変換を適用したことに基づいて算出される伝達関数に対応する前記特性方程式を用いる、請求項１から５のいずれか一項に記載の測定装置。
7. 前記第１算出部は、前記伝達関数に基づく前記特性方程式を因数分解することにより、前記特性方程式が極値となる解を算出する、請求項６に記載の測定装置。
8. 前記探索部は、探索した解の偏角に基づき、探索した解に対応する前記ビート周波数を算出する、請求項１から７のいずれか一項に記載の測定装置。
9. レーザ装置から複数のモードの周波数変調レーザ光を出力するステップと、
   前記レーザ装置が出力する前記周波数変調レーザ光を、前記周波数変調レーザ光の一部の参照光と、残りの少なくとも一部の測定光とに分岐させるステップと、
   前記測定光を計測対象物に照射して反射された反射光と、前記参照光とを混合してビート信号を発生させるステップと、
   前記ビート信号の時系列データを予め定められた多項式に当てはめた場合のモデルパラメータを推定するステップと、
   前記モデルパラメータに基づいて、前記多項式に対応する特性方程式が極値となる解を算出するステップと、
   算出された解のうち、絶対値が最大となる解を探索し、探索した解に対応するビート周波数を算出するステップと、
   前記ビート周波数に基づき、前記参照光と前記測定光との伝搬距離の差を算出するステップと
   を備える、測定方法。
   

Images