特許 7584350 信号ベクトル導出装置、方法、プログラム、記録媒体

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-11-07

(45)【発行日】2024-11-15

(54)【発明の名称】信号ベクトル導出装置、方法、プログラム、記録媒体

(51)【国際特許分類】

   G01R 33/02 20060101AFI20241108BHJP

   A61B 5/242 20210101ALI20241108BHJP

【ＦＩ】

G01R33/02 R

A61B5/242

【請求項の数】 16

(21)【出願番号】P 2021074991

(22)【出願日】2021-04-27

(65)【公開番号】P2022169141

(43)【公開日】2022-11-09

【審査請求日】2024-02-09

(73)【特許権者】

【識別番号】390005175

【氏名又は名称】株式会社アドバンテスト

(74)【代理人】

【識別番号】100097490

【弁理士】

【氏名又は名称】細田 益稔

(74)【代理人】

【識別番号】100113354

【弁理士】

【氏名又は名称】石井 総

(72)【発明者】

【氏名】緒方 祐史

(72)【発明者】

【氏名】柳田 朋則

【審査官】島田 保

(56)【参考文献】

【文献】米国特許出願公開第２００９／００９３９６４（ＵＳ，Ａ１）

【文献】米国特許出願公開第２００６／０２５１３０３（ＵＳ，Ａ１）

【文献】中国特許出願公開第１０６９７０３４８（ＣＮ，Ａ）

【文献】特開２００３－３８４５５（ＪＰ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｒ ３３／０２

Ａ６１Ｂ ５／２４２

(57)【特許請求の範囲】

【請求項1】

所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出装置であり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例する信号ベクトル導出装置であって、
前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出部と、
前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出部と、
を備えた信号ベクトル導出装置。

【請求項2】

請求項１に記載の信号ベクトル導出装置であって、
前記ベクトルが、磁気双極子モーメントまたは電気双極子モーメントである、
信号ベクトル導出装置。

【請求項3】

請求項１に記載の信号ベクトル導出装置であって、
前記ベクトルが、電気双極子モーメントであり、
前記測定結果の成分と同じ方向の前記ベクトルの成分が０である、
信号ベクトル導出装置。

【請求項4】

請求項１に記載の信号ベクトル導出装置であって、
前記第一係数が、前記ボクセルの各々と前記センサの各々との間の位置関係に基づき定められる、
信号ベクトル導出装置。

【請求項5】

請求項１に記載の信号ベクトル導出装置であって、
前記第二係数のいずれか一つまたは二つが０である、
信号ベクトル導出装置。

【請求項6】

請求項１に記載の信号ベクトル導出装置であって、
前記スペクトルが、ＭＵＳＩＣ法に則して求められる、
信号ベクトル導出装置。

【請求項7】

請求項６に記載の信号ベクトル導出装置であって、
前記スペクトル導出部が、前記センサの測定結果から求められたノイズ部分空間の固有ベクトルに基づき前記スペクトルを導出する、
信号ベクトル導出装置。

【請求項8】

請求項６に記載の信号ベクトル導出装置であって、
前記第一係数に第二係数を乗じて得られた値の和が、ＭＵＳＩＣ法における伝達関数である、
信号ベクトル導出装置。

【請求項9】

請求項１に記載の信号ベクトル導出装置であって、
前記極大値が２つ以上存在する、
信号ベクトル導出装置。

【請求項10】

請求項１に記載の信号ベクトル導出装置であって、
前記スペクトルに基づき、前記信号源の存在するボクセルの位置を導出する位置導出部を備えた信号ベクトル導出装置。

【請求項11】

請求項１０に記載の信号ベクトル導出装置であって、
前記位置導出部が、各ボクセルにおいて前記スペクトルの各々がとる値の最大値に基づき、前記信号源の存在するボクセルの位置を導出する信号ベクトル導出装置。

【請求項12】

請求項１１に記載の信号ベクトル導出装置であって、
前記位置導出部が、
前記最大値の内でも最大の値から所定の範囲内の値をとるボクセルの重心を、前記所定の範囲を増やしながら、前記重心が所定量を超えて変化する回数に１を加えた値が前記信号源の個数になるまで、求め、
前記重心を求めた前記ボクセルを、前記信号源の個数だけクラスタリングし、
クラスタリングされた前記ボクセルの内で、前記スペクトルが最大の値をとるものを前記信号源の存在するボクセルの位置とする、
信号ベクトル導出装置。

【請求項13】

請求項１０に記載の信号ベクトル導出装置であって、
前記位置導出部が既に導出した前記信号源の存在するボクセルの位置に基づき、さらに、ボクセルの大きさを小さくして、前記信号源の存在するボクセルの位置を導出する、
信号ベクトル導出装置。

【請求項14】

所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出方法であり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例する信号ベクトル導出方法であって、
前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出工程と、
前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出工程と、
を備えた信号ベクトル導出方法。

【請求項15】

所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出処理をコンピュータに実行させるためのプログラムであり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例するプログラムであって、
前記信号ベクトル導出処理が、
前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出工程と、
前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出工程と、
を備えたプログラム。

【請求項16】

所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出処理をコンピュータに実行させるためのプログラムを記録したコンピュータによって読み取り可能な記録媒体であり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例する記録媒体であって、
前記信号ベクトル導出処理が、
前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出工程と、
前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出工程と、
を備えた記録媒体。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、磁場などの信号の測定に関する。

【背景技術】

【0002】

従来より、磁場の測定結果から信号源の位置を推定する方法（例えば、部分空間法(MUSIC法, SF, WSF等)）が知られている（特許文献１、非特許文献１および非特許文献２を参照）。また、生体に関する磁場を測定することも知られている（特許文献２、特許文献３および特許文献４を参照）。

【先行技術文献】

【特許文献】

【0003】

【文献】特表２００９－５３４１０３号公報

【文献】特開２００９－１７２０８８号公報

【文献】特開２００７－２０５９４号公報

【文献】特開２００２－２８１４３号公報

【文献】Ishii H., Niki N., Nakaya Y., Nishitani H., Kang Y.M. (2000)“Algorithm for Magnetocardiography Analysis Based on Multi Channel SQUID andThoracic MR Images” In: Aine C.J., Stroink G., Wood C.C., Okada Y., SwithenbyS.J. (eds) Biomag 96. Springer, New York, NY. pp. 261-262.

【文献】H. Kudo; T. Maemura; T. Saito “Multiple Signal Source Localizationfrom Spatio-Temporal Magnetocardiogram” Proceedings of 1994 IEEE NuclearScience Symposium - NSS'94 pp. 1832-1836

【発明の概要】

【発明が解決しようとする課題】

【0004】

本発明は、磁場などの信号の測定精度の向上を課題とする。

【課題を解決するための手段】

【0005】

本発明にかかる信号ベクトル導出装置は、所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出するものであり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例し、前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出部と、前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出部とを備えるように構成される。

【0006】

上記のように構成された信号ベクトル導出装置は、所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出するものである。前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例する。スペクトル導出部が、前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出する。向き導出部が、前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する。

【0007】

なお、本発明にかかる信号ベクトル導出装置は、前記ベクトルが、磁気双極子モーメントまたは電気双極子モーメントであるようにしてもよい。

【0008】

なお、本発明にかかる信号ベクトル導出装置は、前記ベクトルが、電気双極子モーメントであり、前記測定結果の成分と同じ方向の前記ベクトルの成分が０であるようにしてもよい。

【0009】

なお、本発明にかかる信号ベクトル導出装置は、前記第一係数が、前記ボクセルの各々と前記センサの各々との間の位置関係に基づき定められるようにしてもよい。

【0010】

なお、本発明にかかる信号ベクトル導出装置は、前記第二係数のいずれか一つまたは二つが０であるようにしてもよい。

【0011】

なお、本発明にかかる信号ベクトル導出装置は、前記スペクトルが、ＭＵＳＩＣ法に則して求められるようにしてもよい。

【0012】

なお、本発明にかかる信号ベクトル導出装置は、前記スペクトル導出部が、前記センサの測定結果から求められたノイズ部分空間の固有ベクトルに基づき前記スペクトルを導出するようにしてもよい。

【0013】

なお、本発明にかかる信号ベクトル導出装置は、前記第一係数に第二係数を乗じて得られた値の和が、ＭＵＳＩＣ法における伝達関数であるようにしてもよい。

【0014】

なお、本発明にかかる信号ベクトル導出装置は、前記極大値が２つ以上存在するようにしてもよい。

【0015】

なお、本発明にかかる信号ベクトル導出装置は、前記スペクトルに基づき、前記信号源の存在するボクセルの位置を導出する位置導出部を備えるようにしてもよい。

【0016】

なお、本発明にかかる信号ベクトル導出装置は、前記位置導出部が、各ボクセルにおいて前記スペクトルの各々がとる値の最大値に基づき、前記信号源の存在するボクセルの位置を導出するようにしてもよい。

【0017】

なお、本発明にかかる信号ベクトル導出装置は、前記位置導出部が、前記最大値の内でも最大の値から所定の範囲内の値をとるボクセルの重心を、前記所定の範囲を増やしながら、前記重心が所定量を超えて変化する回数に１を加えた値が前記信号源の個数になるまで、求め、前記重心を求めた前記ボクセルを、前記信号源の個数だけクラスタリングし、クラスタリングされた前記ボクセルの内で、前記スペクトルが最大の値をとるものを前記信号源の存在するボクセルの位置とするようにしてもよい。

【0018】

なお、本発明にかかる信号ベクトル導出装置は、前記位置導出部が既に導出した前記信号源の存在するボクセルの位置に基づき、さらに、ボクセルの大きさを小さくして、前記信号源の存在するボクセルの位置を導出するようにしてもよい。

【0019】

本発明は、所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出方法であり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例する信号ベクトル導出方法であって、前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出工程と、前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出工程とを備えた信号ベクトル導出方法である。

【0020】

本発明は、所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出処理をコンピュータに実行させるためのプログラムであり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例するプログラムであって、前記信号ベクトル導出処理が、前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出工程と、前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出工程とを備えたプログラムである。

【0021】

本発明は、所定の方向を有するベクトルにより表される信号を受け、互いに直交する３軸の成分を測定する複数のセンサの測定結果を受けて前記ベクトルの向きを導出する信号ベクトル導出処理をコンピュータに実行させるためのプログラムを記録したコンピュータによって読み取り可能な記録媒体であり、前記センサの測定結果が、前記ベクトルの前記３軸の成分の各々に第一係数を乗じたものの和に比例する記録媒体であって、前記信号ベクトル導出処理が、前記センサの測定結果と、前記第一係数に第二係数を乗じて得られた値の和とに基づき得られるスペクトルであって、前記信号を出力する信号源の存在するボクセルで極大値をとるスペクトルを導出するスペクトル導出工程と、前記スペクトルを得るために用いられた前記第二係数に基づき、前記ベクトルの向きを導出する向き導出工程とを備えた記録媒体である。

【図面の簡単な説明】

【0022】

【図1】本発明の実施形態にかかるボクセルＶおよび磁気センサＭＳの斜視図である。

【図2】本発明の実施形態にかかる信号ベクトル導出装置１の構成を示す機能ブロック図である。

【図3】最大値Ｐのグラフの一例である。

【発明を実施するための形態】

【0023】

以下、本発明の実施形態を図面を参照しながら説明する。

【0024】

図１は、本発明の実施形態にかかるボクセルＶおよび磁気センサＭＳの斜視図である。図２は、本発明の実施形態にかかる信号ベクトル導出装置１の構成を示す機能ブロック図である。

【0025】

図１を参照して、信号源Ｓ１および信号源Ｓ２が信号を出力する。信号は、所定の方向を有するベクトルｍにより表される。ベクトルｍは、例えば、磁気双極子モーメントである。なお、信号源の個数は、例えば２個であるが、磁気センサＭＳの個数未満であれば３個以上であってもよい。ただし、各信号源の出力する信号は、互いに周波数または位相が異なっているものとする。

【0026】

また、信号源Ｓ１および信号源Ｓ２が存在する空間内の位置は、ボクセルＶ（例えば、１０×１０×１０＝１０００個のボクセル）により表される。信号源Ｓ１および信号源Ｓ２は、それぞれ異なるボクセルＶに位置している。なお、１０００個のボクセルＶの各々を、Ｖ１～Ｖ１０００と表記する。

【0027】

複数（例えば、８行８列の６４個）の磁気センサＭＳは、信号（例えば、磁気双極子モーメント）を受け、互いに直交する３軸Ｘ、Ｙ、Ｚの成分Ｂｘ、Ｂｙ、Ｂｚを測定する。なお、６４個の磁気センサＭＳの各々を、ＭＳ１～ＭＳ６４と表記する。

【0028】

ここで、ベクトルｒを、信号源（磁気双極子）から磁気センサＭＳまでの方向ベクトルとすると、磁気センサＭＳにより測定される磁束密度Ｂ（ベクトルｒの関数）は、ビオサバールの法則により、式（１）のように表される。ただし、μ₀は、磁気定数である。また、ベクトルｒは、ボクセルＶの各々（Ｖ１～Ｖ１０００）と、磁気センサＭＳの各々ＭＳ１～ＭＳ６４との間の位置関係といえる。

【0029】

【数1】

式（１）より、Ｂｘは、以下の式（２）のように表される。ただし、ｒｘ、ｒｙ、ｒｚは、それぞれ、ベクトルｒのＸ成分、Ｙ成分、Ｚ成分である。また、ｍｘ、ｍｙ、ｍｚは、それぞれ、ベクトルｍのＸ成分、Ｙ成分、Ｚ成分である。

【0030】

【数2】

ここで、式（２）のｍｘの係数をｖｘ１、ｍｙの係数をｖｘ２、ｍｚの係数をｖｘ３とすると、式（２）は、式（２’）のように表される。すると、磁気センサＭＳの測定結果Ｂｘは、ベクトルｍの３軸Ｘ、Ｙ、Ｚの成分ｍｘ、ｍｙ、ｍｚの各々にｖｘ１、ｖｘ２、ｖｘ３（第一係数）を乗じたものの和（ｖｘ１ｍｘ＋ｖｘ２ｍｙ＋ｖｘ３ｍｚ）に比例することなる。

【0031】

式（１）より、Ｂｙは、以下の式（３）のように表される。

【0032】

【数3】

ここで、式（３）のｍｘの係数をｖｙ１、ｍｙの係数をｖｙ２、ｍｚの係数をｖｙ３とすると、式（３）は、式（３’）のように表される。すると、磁気センサＭＳの測定結果Ｂｙは、ベクトルｍの３軸Ｘ、Ｙ、Ｚの成分ｍｘ、ｍｙ、ｍｚの各々にｖｙ１、ｖｙ２、ｖｙ３（第一係数）を乗じたものの和（ｖｙ１ｍｘ＋ｖｙ２ｍｙ＋ｖｙ３ｍｚ）に比例することなる。

【0033】

式（１）より、Ｂｚは、以下の式（４）のように表される。

【0034】

【数4】

ここで、式（４）のｍｘの係数をｖｚ１、ｍｙの係数をｖｚ２、ｍｚの係数をｖｚ３とすると、式（４）は、式（４’）のように表される。すると、磁気センサＭＳの測定結果Ｂｚは、ベクトルｍの３軸Ｘ、Ｙ、Ｚの成分ｍｘ、ｍｙ、ｍｚの各々にｖｚ１、ｖｚ２、ｖｚ３（第一係数）を乗じたものの和（ｖｚ１ｍｘ＋ｖｚ２ｍｙ＋ｖｚ３ｍｚ）に比例することなる。

【0035】

なお、ｖｘ１、ｖｘ２、ｖｘ３、ｖｙ１、ｖｙ２、ｖｙ３、ｖｚ１、ｖｚ２、ｖｚ３（第一係数）は、式（２）～（４）および式（２’）～（４’）を参照して、ベクトルｒに基づき定められている。

【0036】

図２を参照して、本発明の実施形態にかかる信号ベクトル導出装置１は、相対位置記録部１１、第一係数導出部１２、伝達関数導出部１３、雑音固有ベクトル導出部１４、スペクトル導出部１６、向き導出部１８、位置導出部１９を備える。

【0037】

信号ベクトル導出装置１は、複数のセンサＭＳ１～ＭＳ６４の測定結果を受けてベクトルｍの向きを導出する。

【0038】

相対位置記録部１１は、ボクセルＶの各々１０００個と、磁気センサＭＳの各々ＭＳ１～ＭＳ６４との間の相対位置であるベクトルｒを記録する。

【0039】

第一係数導出部１２は、相対位置記録部１１からベクトルｒを読み出し、第一係数ｖｘ１、ｖｘ２、ｖｘ３、ｖｙ１、ｖｙ２、ｖｙ３、ｖｚ１、ｖｚ２、ｖｚ３を導出する（式（２）～（４）および式（２’）～（４’）を参照）。

【0040】

例えば、第一係数ｖｘ１は、以下の式（５）のように表される。

【0041】

【数5】

ベクトルｒは、ボクセルＶの位置と磁気センサＭＳの位置とによって定まるので、１０００×６４通りの値をとる。よって、第一係数ｖｘ１も、１０００×６４通りの値をとる。式（５）においては、１行目に磁気センサＭＳ１に関するｖｘ１、２行目に磁気センサＭＳ２に関するｖｘ１、…、６４行目に磁気センサＭＳ６４に関するｖｘ１を表記している。さらに、式（５）においては、１列目にボクセルＶ１に関するｖｘ１、２列目にボクセルＶ２に関するｖｘ１、…、１０００列目にボクセルＶ１０００に関するｖｘ１を表記している。例えば、式（５）の１行１０００列目の要素ｖｘ１（１，１０００）は、磁気センサＭＳ１およびボクセルＶ１０００に関するｖｘ１を意味する。すなわち、ベクトルｒを、ボクセルＶ１０００から磁気センサＭＳ１までの方向ベクトルとして、式（２）のｍｘの係数に代入すると、ｖｘ１（１，１０００）を求めることができる。

【0042】

他の第一係数ｖｘ２、ｖｘ３、ｖｙ１、ｖｙ２、ｖｙ３、ｖｚ１、ｖｚ２、ｖｚ３も同様に、１０００×６４通りの値をとる。

【0043】

なお、第一係数ｖｘ１は、そのまま用いることも考えられるが、以下の式（６）のように正規化して、以後の処理に使用する。

【0044】

【数6】

ただし、ｈは行、ｎは列を示す。すなわち、ｈ行ｎ列目の第一係数ｖｘ１を、１、２、…６４行ｎ列目の第一係数ｖｘ１を２乗して合計し１／２乗したもので割ったものを新たにｈ行ｎ列目の第一係数ｖｘ１とする。

【0045】

他の第一係数ｖｘ２、ｖｘ３、ｖｙ１、ｖｙ２、ｖｙ３、ｖｚ１、ｖｚ２、ｖｚ３も同様に、正規化する。

【0046】

第一係数導出部１２は、上記のようにして正規化した第一係数を出力する。

【0047】

雑音固有ベクトル導出部１４は、磁気センサＭＳの測定結果Ｂｘ、Ｂｙ、Ｂｚから、ＭＵＳＩＣ法に則して、ノイズ部分空間の固有ベクトルを求める。

【0048】

まず、磁気センサＭＳの測定結果Ｂｘから以下の式（７）のように、X(t)xを求める。ただし、ｔは測定を行った時間である。Ｔは転置を意味する。

【0049】

【数7】

X(t)xは、１行目に時間ｔ１に測定されたＢｘ、２行目に時間ｔ２に測定されたＢｘ、…、Ｎ行目に時間ｔＮに測定されたＢｘを記載し、かつ、１列目に磁気センサＭ１により測定されたＢｘ、２列目に磁気センサＭ２により測定されたＢｘ、…、６４列目に磁気センサＭ６４により測定されたＢｘを記載した行列を転置した行列である。

【0050】

X(t)xを用いて、以下の式（８）のように相関行列を求める。

【0051】

【数8】

ただし、Ｅはアンサンブル平均を意味する。式（８）により、６４行６４列の行列が得られる。式（８）により得られた相関行列の固有値および固有ベクトルを求める。このようにして得られた固有値のうち、信号源の個数（２個）の分は、大きな値となるが、残りの固有値（６４－２＝６２個）は小さな値となる。そこで、小さな値の固有値に対応する固有ベクトルを求め、ノイズ部分空間の固有ベクトルｅｘとする。ノイズ部分空間の固有ベクトルｅｘは、６４行１列のベクトルである。ノイズ部分空間の固有ベクトルｅｘは、小さな値の固有値に対応して６２個存在する。

【0052】

なお、ノイズ部分空間の固有ベクトルｅｙも同様に求めることができる。まず、式（７）のＢｘをＢｙに置き換え、式（７）および式（８）のX(t)xをX(t)yに置き換えて、式（８）により相関行列を求める。あとは、同様に、小さな値の固有値に対応する固有ベクトルを求め、ノイズ部分空間の固有ベクトルｅｙとする。ノイズ部分空間の固有ベクトルｅｙは、６４行１列のベクトルである。ノイズ部分空間の固有ベクトルｅｙは、小さな値の固有値に対応して６２個存在する。

【0053】

また、ノイズ部分空間の固有ベクトルｅｚも同様に求めることができる。まず、式（７）のＢｘをＢｚに置き換え、式（７）および式（８）のX(t)xをX(t)zに置き換えて、式（８）により相関行列を求める。あとは、同様に、小さな値の固有値に対応する固有ベクトルを求め、ノイズ部分空間の固有ベクトルｅｚとする。ノイズ部分空間の固有ベクトルｅｚは、６４行１列のベクトルである。ノイズ部分空間の固有ベクトルｅｚは、小さな値の固有値に対応して６２個存在する。

【0054】

伝達関数導出部１３は、以下の式（９）、式（１０）および式（１１）のように、伝達関数ｖｘ、ｖｙおよびｖｚを導出する。第一係数ｖｘ１、ｖｘ２、ｖｘ３に、それぞれ第二係数ａｘ、ｂｘ、ｃｘを乗じて得られた値の和を導出する（式（９）参照）。導出結果を、伝達関数ｖｘという。第一係数ｖｙ１、ｖｙ２、ｖｙ３に、それぞれ第二係数ａｙ、ｂｙ、ｃｙを乗じて得られた値の和を導出する（式（１０）参照）。導出結果を、伝達関数ｖｙという。第一係数ｖｚ１、ｖｚ２、ｖｚ３に、それぞれ第二係数ａｚ、ｂｚ、ｃｚを乗じて得られた値の和を導出する（式（１１）参照）。導出結果を、伝達関数ｖｚという。

【0055】

【数9】

なお、伝達関数ｖｘ、ｖｙおよびｖｚは、ＭＵＳＩＣ法における伝達関数である
また、第二係数は、いずれも０以外の値であってもよい。例えば、ａｘ＝ｂｘ＝ｃｘ＝１（すなわち、ｖｘ＝ｖｘ１＋ｖｘ２＋ｖｘ３）でもよいし、ａｙ＝ｂｙ＝１、ｃｙ＝－１（すなわち、ｖｙ＝ｖｙ１＋ｖｙ２－ｖｙ３）でもよいし、ａｚ＝１、ｂｚ＝－１、ｃｚ＝１（すなわち、ｖｚ＝ｖｚ１－ｖｚ２＋ｖｚ３）でもよい。

【0056】

ただし、第二係数のいずれか一つまたは二つが０であってもよい。例えば、ａｘ＝１、ｂｘ＝ｃｘ＝０（すなわち、ｖｘ＝ｖｘ１）でもよいし、ａｘ＝ｂｘ＝１、ｃｘ＝０（すなわち、ｖｘ＝ｖｘ１＋ｖｘ２）でもよい。

【0057】

ここで、（ａｋ、ｂｋ、ｃｋ）（ただし、ｋ＝ｘ、ｙ、ｚ）が、（１，０，０）、（０，１，０）、（０，０，１）、（１，１，０）、（１，－１，０）、（０，１，１）、（０，１，－１）、（１，０，１）、（－１，０，１）、（１，１，１）、（－１，１，１）、（１，－１，１）、（１，１，－１）の１３種類あるものとする。これらに対応するｖｋをそれぞれ、ｖｋ１、ｖｋ２、…、ｖｋ１３とする。

【0058】

例えば、（ａｘ、ｂｘ、ｃｘ）が、（１，０，０）である場合、ｖｘ＝ｖｘ１である。（ａｘ、ｂｘ、ｃｘ）が、（１，１，０）である場合、ｖｘ＝ｖｘ４＝ｖｘ１＋ｖｘ２である。（ａｘ、ｂｘ、ｃｘ）が、（１，１，－１）である場合、ｖｘ＝ｖｘ１３＝ｖｘ１＋ｖｘ２－ｖｘ３である。

【0059】

例えば、（ａｙ、ｂｙ、ｃｙ）が、（１，０，０）である場合、ｖｙ＝ｖｙ１である。（ａｙ、ｂｙ、ｃｙ）が、（１，１，０）である場合、ｖｙ＝ｖｙ４＝ｖｙ１＋ｖｙ２である。（ａｙ、ｂｙ、ｃｙ）が、（１，１，－１）である場合、ｖｙ＝ｖｙ１３＝ｖｙ１＋ｖｙ２－ｖｙ３である。

【0060】

例えば、（ａｚ、ｂｚ、ｃｚ）が、（１，０，０）である場合、ｖｚ＝ｖｚ１である。（ａｚ、ｂｚ、ｃｚ）が、（１，１，０）である場合、ｖｚ＝ｖｚ４＝ｖｚ１＋ｖｚ２である。（ａｚ、ｂｚ、ｃｚ）が、（１，１，－１）である場合、ｖｚ＝ｖｚ１３＝ｖｚ１＋ｖｚ２－ｖｚ３である。

【0061】

スペクトル導出部１６は、信号源Ｓ１、Ｓ２の存在するボクセルＶで極大値をとるスペクトルを導出する。このようなスペクトルは、ＭＵＳＩＣ法に則して求められる。スペクトルの極大値は、信号源の個数に対応して、２つある。なお、信号源の個数が３つ以上あれば、スペクトルの極大値も３つ以上ある。

【0062】

スペクトルは、磁気センサＭＳの測定結果Ｂｘ、Ｂｙ、Ｂｚ（から求められたノイズ部分空間の固有ベクトルｅｘ、ｅｙ、ｅｚ）と、第一係数に第二係数を乗じて得られた値の和（すなわち、伝達関数ｖｘ、ｖｙ、ｖｚ）（式（９）、（１０）、（１１））とに基づき、スペクトル導出部１６により導出される。スペクトル導出部１６は、伝達関数導出部１３の出力する伝達関数ｖｘ、ｖｙ、ｖｚと、雑音固有ベクトル導出部１４の出力するノイズ部分空間の固有ベクトルｅｘ、ｅｙ、ｅｚとに基づきスペクトルを導出する。

【0063】

スペクトル導出部１６は、以下のようにして、スペクトルＰｘ１を導出する。
（１）ノイズ部分空間の固有ベクトルｅｘ（６４行１列のベクトル）は、６２個存在するが、これらのベクトルｅｘを６２列並べて、ｅｘを６４行６２列の行列とする。
（２）行列ｅｘを転置し、伝達関数ｖｘ１（６４行１０００列の行列）に乗じる。すなわち、ｅｘ^Ｔｖｘ１を求める。これは、６２行１０００列の行列となる。
（３）（２）で得た行列の各成分を２乗する。
（４）（３）で得た行列の成分を列ごとに合計し、１行に並べて、１行１０００列の行列を得る。例えば、（３）で得た行列の１行Ｑ列目＋２行Ｑ列目＋…＋６２行Ｑ列目が、（４）で得られる１行１０００列の行列の１行Ｑ列目の成分となる（ただし、Ｑは１～１０００の整数）。
（５）（４）で得た行列の各成分の逆数をとると、スペクトルＰｘ１（１行１０００列の行列）が得られる。

【0064】

なお、スペクトルＰｘ１の各列は、ボクセルＶ１～Ｖ１０００に対応する。他のスペクトルも同様である。

【0065】

スペクトル導出部１６は、スペクトルＰｘ２、Ｐｘ３、…、Ｐｘ１３も導出する。スペクトルＰｘ２、Ｐｘ３、…、Ｐｘ１３は、上記（２）における伝達関数ｖｘ１を、ｖｘ２、ｖｘ３、…、ｖｘ１３に置き換えれば、導出することができる。

【0066】

スペクトル導出部１６は、スペクトルＰｙ１、Ｐｙ２、Ｐｙ３、…、Ｐｙ１３を導出する。上記（１）におけるノイズ部分空間の固有ベクトルｅｘをｅｙに置き換え、上記（２）における伝達関数ｖｘ１を、ｖｙ１、ｖｙ２、ｖｙ３、…、ｖｙ１３に置き換えれば、スペクトルＰｙ１、Ｐｙ２、Ｐｙ３、…、Ｐｙ１３を導出することができる。

【0067】

スペクトル導出部１６は、スペクトルＰｚ１、Ｐｚ２、Ｐｚ３、…、Ｐｚ１３を導出する。上記（１）におけるノイズ部分空間の固有ベクトルｅｘをｅｚに置き換え、上記（２）における伝達関数ｖｘ１を、ｖｚ１、ｖｚ２、ｖｚ３、…、ｖｚ１３に置き換えれば、スペクトルＰｚ１、Ｐｚ２、Ｐｚ３、…、Ｐｚ１３を導出することができる。

【0068】

位置導出部１９は、スペクトルＰｘ１、Ｐｘ２、Ｐｘ３、…、Ｐｘ１３、Ｐｙ１、Ｐｙ２、Ｐｙ３、…、Ｐｙ１３、Ｐｚ１、Ｐｚ２、Ｐｚ３、…、Ｐｚ１３に基づき、信号源Ｓ１、Ｓ２の存在するボクセルＶの位置を導出する。

【0069】

スペクトル導出部１６の出力するスペクトルは、以下の式（１２）のように表される。

【0070】

【数10】

各ボクセルにおいてスペクトルの各々がとる値の最大値Ｐ（すなわち、式（１２）の各列における最大値）を求める（式（１３）を参照）。

【0071】

【数11】

Ｐにおいて、極大値をとる列（ボクセルに対応）が信号源の個数（２つ）現れるので、その列に対応するボクセルが、信号源Ｓ１、Ｓ２の存在するボクセルである。極大値をとる列の検出法を、以下に説明する。

【0072】

図３は、最大値Ｐのグラフの一例である。図３において、縦軸はスペクトルの値、横軸はボクセル（Ｖ１～Ｖ１０００）を示す。

【0073】

図３を参照して、ボクセルＶ７５０において最大値Ｐ（のスペクトル）の値がＳＰ１（極大値）、ボクセルＶ２５０において最大値Ｐ（のスペクトル）の値がＳＰ２（極大値）をとるものとする。ただし、ＳＰ１がＳＰ２よりも大きいものとする。

【0074】

位置導出部１９は、最大値Ｐの内でも最大の値ＳＰ１から所定の範囲内（例えば、最大値Ｐの値が０．９５ＳＰ１以上）の値をとるボクセルの重心を、所定の範囲を増やしながら（例えば、最大値Ｐの値が０．９５ＳＰ１以上→０．９０ＳＰ１以上→０．８５ＳＰ１以上→…というように０．０５ＳＰ１ずつ所定の範囲を広げる）、重心が所定量を超えて変化する回数に１を加えた値が信号源の個数（２個）になるまで、求める。

【0075】

最大値Ｐの内でも最大の値ＳＰ１から近傍におけるボクセルの重心は、おおむねボクセルＶ７５０である。しかし、最大の値ＳＰ１から所定の範囲が広がって、ＳＰ２を含むようになると、ボクセルの重心がボクセルＶ７５０から、かなり小さくなる。すると、重心が所定量を超えて変化するので、その回数（１回）に１を加えた値が信号源の個数（２個）になるので、ボクセルの重心を求めることを終了する。

【0076】

次に、重心を求めたボクセルを、信号源の個数（２個）だけクラスタリングする。例えば、教師なし機械学習のKmeansクラスタリングを行い、信号源の数だけラベリングする。

【0077】

最後に、クラスタリングされたボクセルの内で、スペクトルが最大の値をとるものを信号源の存在するボクセルの位置とする。

【0078】

なお、位置導出部１９は、このようにして既に導出した信号源の存在するボクセルの位置に基づき、さらに、ボクセルの大きさを小さくして、信号源の存在するボクセルの位置を導出するようにしてもよい。このようにして、信号源の存在するボクセルの位置を、高精度かつ高速に計算することができる。

【0079】

向き導出部１８は、位置導出部１９から信号源Ｓ１、Ｓ２に対応する（すなわち、Ｐにおいて極大値をとる）Ｐｋｊ（ただし、ｋ＝ｘ、ｙ、ｚかつｊ＝１、２、３、…）を受ける。さらに、向き導出部１８は、信号源Ｓ１、Ｓ２に対応するＰｋｊを得るために用いられた第二係数に基づき、ベクトルｍの向きを導出する。

【0080】

例えば、向き導出部１８に、位置導出部１９から信号源Ｓ１およびＳ２に対応するスペクトルとして、７５０列目（ボクセルＶ７５０）のＰｘ１３（Ｐｙ１３またはＰｚ１３）および２５０列目（ボクセルＶ２５０）のＰｘ４（Ｐｘ４またはＰｚ４）が与えられたとする。

【0081】

すると、向き導出部１８は、ボクセルＶ７５０の信号源Ｓ１におけるベクトルｍの向きとして、Ｐｘ１３（Ｐｙ１３またはＰｚ１３）を得るために用いられた第二係数（ａｋ、ｂｋ、ｃｋ）（ただし、ｋ＝ｘ、ｙ、ｚ）は、（１，１，－１）である。よって、向き導出部１８は、ボクセルＶ７５０の信号源Ｓ１におけるベクトルｍの向きを、ベクトル（１，１，－１）に平行であると導出する。ただし、ベクトル（１，１，－１）は、Ｘ成分１、Ｙ成分１、Ｚ成分－１のベクトルである。

【0082】

さらに、向き導出部１８は、ボクセルＶ２５０の信号源Ｓ２におけるベクトルｍの向きとして、Ｐｘ４（Ｐｘ４またはＰｚ４）を得るために用いられた第二係数（ａｋ、ｂｋ、ｃｋ）（ただし、ｋ＝ｘ、ｙ、ｚ）は、（１，１，０）である。よって、向き導出部１８は、ボクセルＶ２５０の信号源Ｓ２におけるベクトルｍの向きを、ベクトル（１，１，０）に平行であると導出する。ただし、ベクトル（１，１，０）は、Ｘ成分１、Ｙ成分１、Ｚ成分０のベクトルである。

【0083】

次に、本発明の実施形態の動作を説明する。

【0084】

第一係数導出部１２により、相対位置記録部１１からベクトルｒが読み出され、第一係数ｖｘ１、ｖｘ２、ｖｘ３、ｖｙ１、ｖｙ２、ｖｙ３、ｖｚ１、ｖｚ２、ｖｚ３が導出される（式（２）～（４）および式（２’）～（４’）を参照）。

【0085】

なお、第一係数は、１０００×６４通りの値をとり（式（５）参照）、正規化されて（式（６）参照）、伝達関数導出部１３に与えられる。

【0086】

伝達関数導出部１３によって、第一係数および第二係数ａｘ、ｂｘ、ｃｘ、ａｙ、ｂｙ、ｃｙ、ａｚ、ｂｚ、ｃｚに基づき、伝達関数ｖｘ、ｖｙ、ｖｚが導出される（式（９）、式（１０）および式（１１）参照）。

【0087】

雑音固有ベクトル導出部１４によって、磁気センサＭＳの測定結果Ｂｘ、Ｂｙ、Ｂｚから、ＭＵＳＩＣ法に則して、ノイズ部分空間の固有ベクトルｅｘ、ｅｙ、ｅｚが導出される。

【0088】

スペクトル導出部１６により、伝達関数ｖｘ、ｖｙ、ｖｚおよびノイズ部分空間の固有ベクトルｅｘ、ｅｙ、ｅｚに基づき、スペクトルＰｘ１、Ｐｘ２、Ｐｘ３、…、Ｐｘ１３、Ｐｙ１、Ｐｙ２、Ｐｙ３、…、Ｐｙ１３、Ｐｚ１、Ｐｚ２、Ｐｚ３、…、Ｐｚ１３が導出される（式（１２）参照）。

【0089】

位置導出部１９により、各ボクセルにおいてスペクトルの各々がとる値の最大値Ｐ（すなわち、式（１２）の各列における最大値）が求められる（式（１３）および図３を参照）。最大値Ｐに基づき、信号源Ｓ１、Ｓ２の存在するボクセル２５０、７５０が導出される。

【0090】

向き導出部１８により、信号源Ｓ１、Ｓ２に対応するＰｋｊを得るために用いられた第二係数に基づき、ベクトルｍの向きが導出される。

【0091】

本発明の実施形態によれば、磁場などの信号の測定精度が向上する。

【0092】

例えば、伝達関数ｖｋが、ｖｋ１、ｖｋ２、ｖｋ３のみである場合、ベクトルｍの向きが、Ｘ方向、Ｙ方向またはＺ方向に平行な場合しか、ベクトルｍの向きを測定できない。ベクトルｍの向きが、それ以外の向き、例えば、ベクトル（１，１，０）（すなわち、Ｘ成分１、Ｙ成分１、Ｚ成分０のベクトル）に平行である場合は、ベクトルｍの向きを測定できない。

【0093】

しかし、本発明の実施形態によれば、伝達関数ｖｋが、ｖｋ１、ｖｋ２、ｖｋ３、…、ｖｋ１３と多くの種類があるので、ベクトルｍの向きが、Ｘ方向、Ｙ方向およびＺ方向に平行ではない場合であっても、ベクトルｍの向きを測定できる。

【0094】

なお、本発明の実施形態においては、信号のベクトルを磁気双極子モーメントであるとしてきたが、信号のベクトルは磁気双極子モーメントに限定されない。信号のベクトルは、例えば、電気双極子モーメント（ベクトルｐ）であってもよい。

【0095】

磁気センサＭＳにより測定される磁束密度Ｂ（ベクトルｒの関数）は、式（１４）のように表される。

【0096】

【数12】

式（１４）より、Ｂｘは、以下の式（１５）のように表される。ただし、ｐｘ、ｐｙ、ｐｚは、それぞれ、ベクトルｐのＸ成分、Ｙ成分、Ｚ成分である。

【0097】

【数13】

ここで、式（１５）のｐｘの係数をｖｘ１、ｐｙの係数をｖｘ２、ｐｚの係数をｖｘ３とすると、式（１５）は、式（１５’）のように表される。すると、磁気センサＭＳの測定結果Ｂｘは、ベクトルｐの３軸Ｘ、Ｙ、Ｚの成分ｐｘ、ｐｙ、ｐｚの各々にｖｘ１、ｖｘ２、ｖｘ３（第一係数）を乗じたものの和（ｖｘ１ｐｘ＋ｖｘ２ｐｙ＋ｖｘ３ｐｚ）に比例することなる。ただし、測定結果Ｂｘの成分と同じ方向（Ｘ方向）のベクトルの成分（ｐｘ）は０であり、それに乗じられる第一係数ｖｘ１は１である。

【0098】

式（１４）より、Ｂｙは、以下の式（１６）のように表される。

【0099】

【数14】

ここで、式（１６）のｐｘの係数をｖｙ１、ｐｙの係数をｖｙ２、ｐｚの係数をｖｙ３とすると、式（１６）は、式（１６’）のように表される。すると、磁気センサＭＳの測定結果Ｂｙは、ベクトルｐの３軸Ｘ、Ｙ、Ｚの成分ｐｘ、ｐｙ、ｐｚの各々にｖｙ１、ｖｙ２、ｖｙ３（第一係数）を乗じたものの和（ｖｙ１ｐｘ＋ｖｙ２ｐｙ＋ｖｙ３ｐｚ）に比例することなる。ただし、測定結果Ｂｙの成分と同じ方向（Ｙ方向）のベクトルの成分（ｐｙ）は０であり、それに乗じられる第一係数ｖｙ２は１である。

【0100】

式（１４）より、Ｂｚは、以下の式（１７）のように表される。

【0101】

【数15】

ここで、式（１７）のｐｘの係数をｖｚ１、ｐｙの係数をｖｚ２、ｐｚの係数をｖｚ３とすると、式（１７）は、式（１７’）のように表される。すると、磁気センサＭＳの測定結果Ｂｘは、ベクトルｐの３軸Ｘ、Ｙ、Ｚの成分ｐｘ、ｐｙ、ｐｚの各々にｖｚ１、ｖｚ２、ｖｚ３（第一係数）を乗じたものの和（ｖｚ１ｐｘ＋ｖｚ２ｐｙ＋ｖｚ３ｐｚ）に比例することなる。ただし、測定結果Ｂｚの成分と同じ方向（Ｚ方向）のベクトルの成分（ｐｚ）は０であり、それに乗じられる第一係数ｖｚ３は１である。

【0102】

信号ベクトル導出装置１の構成および動作は、信号のベクトルが磁気双極子モーメント（ベクトルｍ）の場合と同様であり、説明を省略する。

【0103】

また、上記の実施形態は、以下のようにして実現できる。ＣＰＵ、ハードディスク、メディア（ＵＳＢメモリ、ＣＤ－ＲＯＭなど）読み取り装置を備えたコンピュータに、上記の各部分、例えば相対位置記録部１１、第一係数導出部１２、伝達関数導出部１３、雑音固有ベクトル導出部１４、スペクトル導出部１６、向き導出部１８および位置導出部１９を実現するプログラムを記録したメディアを読み取らせて、ハードディスクにインストールする。このような方法でも、上記の機能を実現できる。

【符号の説明】

【0104】

１ 信号ベクトル導出装置
１１ 相対位置記録部
１２ 第一係数導出部
１３ 伝達関数導出部
１４ 雑音固有ベクトル導出部
１６ スペクトル導出部
１８ 向き導出部
１９ 位置導出部
ＭＳ 磁気センサ
Ｖ ボクセル
Ｂ 磁束密度
ｖｘ１、ｖｘ２、ｖｘ３、ｖｙ１、ｖｙ２、ｖｙ３、ｖｚ１、ｖｚ２、ｖｚ３ 第一係数
ａｋ、ｂｋ、ｃｋ 第二係数
ｖｋ１、ｖｋ２、…、ｖｋ１３ 伝達関数
Ｓ１、Ｓ２ 信号源
ｍ ベクトル（磁気双極子モーメント）
ｅｘ、ｅｙ、ｅｚ ノイズ部分空間の固有ベクトル
Ｐ 最大値

【図1】

【図2】

【図3】