JP2018004286A - 信号処理装置、信号処理方法、信号処理プログラム、及び磁場計測システム - Google Patents

Abstract

【課題】より少ない計算コストで外来ノイズを除去することのできる信号処理技術を提供する。【解決手段】信号処理装置は、センサから測定データを取得するデータ取得部と、前記センサの感度限界を算出する感度限界決定部と、前記感度限界に基づいて径方向に可変の空間分割を行う空間分割部と、前記空間分割を用いて信号源を推定する推定部と、前記信号源の推定結果を用いてノイズを除去するノイズ除去部と、を有する。【選択図】図５

Description

本発明は信号処理の分野に関し、特に磁場解析の技術に関する。

脳磁計、脊髄計、心磁計など、非襲撃的に生体磁気を測定し評価する装置が普及しつつある。生体磁場は微弱であり、測定時に外来ノイズを除去する技術が求められる。たとえば、脳磁計ではＳＳＳ（Signal Space Separation：信号空間分離）法が用いられている。ＳＳＳ法は、信号処理により外来ノイズを除去する手法である。複数のセンサで測定された磁気信号を媒介変数表示し、線形独立なベクトルの基底において測定信号ベクトルの成分を計算することにより、中心に近い容積部分からの信号と中心から遠い容積部分からの信号を分離して、外来ノイズを除去する（たとえば、特許文献１参照）。

公知のＳＳＳ法では、磁気信号を媒介変数表示してから線形独立である基底を得るために、球面調和関数による級数展開を用いている。球面調和関数による級数展開を用いて線形独立である基底を得るための計算が複雑であり、計算コストが大きい。

そこで、より少ない計算コストで外来ノイズを除去することのできる信号処理技術を提供することを目的とする。

上記の課題を解決するために、実施形態では、センサの感度限界に基づいて空間分割を行うことで、効率的な磁場源の推定とノイズ除去を実現する。

本発明の一態様では、信号処理装置は、
センサから測定データを取得するデータ取得部と、
前記センサの感度限界を算出する感度限界決定部と、
前記感度限界に基づいて径方向に可変の空間分割を行う空間分割部と、
前記空間分割を用いて信号源を推定する推定部と、
前記信号源の推定結果を用いてノイズを除去するノイズ除去部と、
を有する。

より少ない計算コストで外来ノイズを除去する信号処理技術が実現できる。

本発明が適用される磁場計測システムの外観図である。 デュワ内部のセンサの配置例を示す図である。 センサの構成例を示す図である。 第１実施形態の空間分割を説明する図である。 第１実施形態の空間分割を説明する図である。 第１実施形態の空間分割の感度限界の算出に用いたセンサ配置を示す図である。 第１実施形態の空間分割の感度限界の算出結果を示す図である。 第１実施形態の磁場解析方法のフローチャートである。 実施形態の磁場解析装置のハードウェア構成図である。 実施形態の磁場解析装置の機能ブロック図である。 第２実施形態の磁場解析方法のフローチャートである。 図１２のスパース解の算出方法のフローチャートである。 第２実施形態の磁場解析装置の機能ブロック図である。

図１は、本発明が適用される磁場計測システム１の一例としての脳磁計測システムの外観図である。磁場計測システム１は、測定装置３と、データ収録サーバ４２と、信号処理装置の一例としての磁場解析装置２０を含む。測定装置３は、液体ヘリウムを用いた極低温環境の保持容器であるデュワ３０を有し、デュワ３０の内部に多数の磁場センサが配置されている。

デュワ３０の前面に、被測定者の頭部に適合する曲面を有する計測面３１が形成されている。計測時に、被測定者は測定テーブル４に仰向けで横たわり、頭部を計測面３１に接触または近接させる。デュワ３０の計測面３１の内側に配置されるセンサで得られた脳磁信号は、データ収録サーバ４２に収録される。データ収録サーバ４２に収録されたデータは、磁場解析装置２０で解析され、出力される。

図１の例では、データ収録サーバ４２と磁場解析装置２０を個別の装置としているが、一体化して磁場解析装置としてもよい。なお、センサを内蔵するデュワ３０と測定テーブル４は一般的に磁気シールドルーム内に配置されているが、図示の便宜上、磁気シールドルームを省略している。

図２は、デュワ３０の側断面図であり、デュワ３０の内部の磁場検出センサの配置例を示す。計測面３１の内側に配置されているセンサ１０−１〜１０−ｍ（ｍは任意の自然数）の各々は、たとえばＳＱＵＩＤ（Superconduction Quantum Interference Device：超伝導量子干渉素子）を用いた磁場センサである。各センサ１０は、その感度方向が、計測面３１で形成される球状の空間の中心を向くように互いに隣接して配置されている。

図３は、センサ１０の概略図である。センサ１０は、磁束を捕捉するピックアップコイル１５と、ピックアップコイル１５で検出された磁束を電圧に変換するＳＱＵＩＤチップ１４を有する。この例では、ピックアップコイル１５として、互いに逆相に巻かれた（差動構成された）一対のコイル１５ａ、１５ｂを有するグラジオメータを用いる。コイル１５ａは信号検出用のコイルであり、コイル１５ｂは参照用のコイルである。コイル１５ａとコイル１５ｂの間の距離ｈは、頭皮から脳内の信号源までの距離を考慮して、たとえば５０ｍｍに設定されている。コイル１５ａの直径は、有意なダイポール情報を損失なく測定できる空間分解能を考慮して、たとえば１５〜２０ｍｍに設定されている。ピックアップコイル１５は、たとえば円筒形のボビンにニオブ（Ｎｂ）線を同軸状に巻いて形成される。逆相に巻かれたコイル１５ａ、１５ｂにより逆向きの遮蔽電流が誘起され、一様な外来ノイズは相殺されて、信号源からの磁場の差分が検出される。

磁気シールドルームや差動コイルのグラジオメータを用いたとしても、外来ノイズを完全に排除できるわけではない。実施形態では後述するように効果的な空間分割と信号処理により、残存する外来ノイズをさらに低減する。

ピックアップコイル１５の出力はＳＱＵＩＤチップ１４に接続される。ＳＱＵＩＤチップ１４は、たとえばシリコン（Ｓｉ）基板に超伝導材料で形成された入力コイル１２、ＳＱＵＩＤリング１１、及びフィードバックコイル１３を有する。入力コイル１２は、ピックアップコイル１５で感知された磁束をＳＱＵＩＤリング１１に伝達する。ＳＱＵＩＤリング１１は、図３の構成例では２個のジョセフソン接合を有する。ＳＱＵＩＤリング１１には、ジョセフソン接合の臨界電流値付近のバイアス電流が流されており、磁束の出入りに従って電圧が誘起される。磁束変化と出力電圧の関係は線形でないため、外部のＦＬＬ（Flux Locked Loop）回路１８で負帰還を行う。ＦＬＬ回路１８は、ＳＱＵＩＤリング１１の出力電圧を積分し、積分電圧を出力するとともに、フィードバックコイル１３を介して積分電圧を帰還磁束としてＳＱＵＩＤリング１１に返す。これにより、ＳＱＵＩＤリング１１の鎖交磁束はゼロにロックされ、ＦＬＬ回路１８の出力電圧は検知された信号磁束に比例したものとなる。

ＦＬＬ回路１８は、一般に磁気シールドルームの外部に配置され、ＦＬＬ回路１８の出力値が、データ収録サーバ４２に記録される。磁場解析装置２０は、データ収録サーバ４２に記録されたデータを用いて磁場解析を行い、たとえば脳内の磁場分布に基づいて神経細胞の活動、すなわち電流の発生源を推定する。

以下で、具体的な信号処理（磁場解析）の手法を説明する。
＜第１実施形態＞
図４及び図５は、第１実施形態の磁場解析で行われる空間分割を説明する図である。図４は三次元の空間分割を示す図、図５は実施形態の空間分割を平面的に示す図である。

第１実施形態では、計算が複雑で高コストの球面調和関数による級数展開を用いずに、センサ１０の感度に応じた空間分割を行って空間フィルタ法を適用し、外来ノイズを除去する。

センサ１０の感度に応じた空間分割では、三次元空間の解析で一般に用いられる均等な単位空間を用いず、磁場源からセンサ１０までの距離に従って有限個の空間分割を行う。これは、磁場源からセンサ１０までの距離によってセンサ１０の感度が相対的に変わるという知見に基づく。良好な実施態様では、センサ１０の感度限界を算出して、感度限界に基づいて径方向に可変の空間分割を行う。実施形態の空間分割により、無限遠までを含めて、中心に近い容積部分からの信号と、中心から遠い容積部分からの信号を分離して、効率的に外来ノイズを除去することができる。

図４の空間図において、太い実線はセンサ１０が配置されている球面Ｓ１を示す。点線は、センサ１０の配置面よりも内側、すなわち、脳内空間での球面Ｓ２を示す。細い実線は、センサ１０の配置面よりも外側の空間での球面Ｓ３を示す。一般に外来ノイズは、センサ１０の配置面よりも外側の空間から到来する。

図４の容積部分ｐ１、ｐ２で示すように、径方向の長さと立体角によって空間を分割することができる。図４では２個の容積部分ｐ１、ｐ２だけが示されているが、全空間にわたって分割される。しかし、径方向に等距離の分割を行う場合、無限遠まで空間分割するのに無限個の分割が必要となる。そこで、ある容積部分に存在する磁場源（電流）に対するセンサ１０の感度限界を考慮し、センサ１０で感度限界を超える距離では、感度限界から無限遠までを１つの分割空間として扱う。また、センサ１０から遠くなるにしたがって径方向の空間分割を粗くする。

より詳細には、ある磁場源からセンサ１０に届く磁界成分と、無限遠側にある他の磁場源からセンサ１０に届く磁界成分とが、センサ１０において区別できなくなる空間位置から無限遠までをまとめて、一つの分割空間にする。以下の説明では、径方向で異なる位置にある２つの磁場源からの磁界成分を、センサ１０で区別して検出（空間的に分解）し得る限界距離を、センサ１０の「感度限界」と称する。これにより、無限遠までの空間を分割する場合でも有限個の空間分割で済ませることができ、計算コストを現実的な範囲に抑えることができる。

図５は、空間分割の概念を説明する図である。太い実線は、センサ１０の配置位置を示す。点線は、センサ１０の配置位置よりも内側の空間分割境界を示す。細い実線は、センサ１０の配置位置よりも外側の空間分割境界を示す。センサ１０の配置位置から遠くなるほど、Ｌ１、Ｌ２、Ｌ３と空間分割の間隔を粗くしている。センサの感度限界となる距離ｒ_lim以降では、無限遠∞までを径方向にひとつの分割空間としている。

具体的には、センサ１０の配置位置から遠くなるほど、センサ１０からの距離ｒの３乗のオーダーで間隔を広げる。センサ１０で検知される磁界の強さは、１／ｒ³で減衰するからである。これについては、図６及び図７を参照して後述する。

ただし、センサ１０の配置位置の近傍では、距離ｒの３乗のオーダーでの可変分割を適用すると分割空間が細かくなりすぎる。また、センサ１０の近傍では、もともと磁場検出感度が相対的に良いので、空間分割を過度に細かくしなくても良好な検出結果を得ることができる。したがって、センサ１０の近傍の所定の範囲では、等間隔（間隔Ｌ１）で空間分割を行うのが望ましい。所定の範囲は、センサ１０の設計感度によって決定される。

図５のように、センサ１０の感度限界以降の空間領域をひとつの分割空間とすることで磁場解析の計算量が増えることを防止できる。また、磁場源に対するセンサ１０の相対感度が高い領域では空間分割の密度を高くし、相対感度が低い領域では空間分割の間隔を広げることで、トータルの計算量をさらに低減することができる。さらに、センサ１０から所定の範囲では、空間分割を等間隔にすることで、空間分割が必要以上に細かくなって計算量が増えることを防止するとともに、磁場検出の精度を担保する。

図６は、実施形態の空間分割におけるセンサ１０の感度限界の算出に用いた配置モデルを示す。センサ１０−０（センサ０）と電流Ｉ₀の距離をｒ、センサ１０−０とセンサ１０−１（センサ１）の距離をａ、電流Ｉ₀と電流Ｉ₁の間隔をｄとする。また、センサ１０−０とセンサ１０−１は同じ構成であり、センサ１０−０とセンサ１０−１の感度は、磁力線の方向と一致しているものとする。

センサ１０−０及びセンサ１０−１の感度限界は、電流Ｉ₀と電流Ｉ₁による磁場がセンサー１０−０で同じ値となるとき、センサ１０−１で検出される電流Ｉ₀と電流Ｉ₁による磁場の差分ΔＭが、センサの感度以下となるときの限界距離と考えることができる。より詳細には、ΔＭはセンサ１０−１で検出される電流Ｉ₀による磁場の傾き（変化量）と電流Ｉ₁による磁場の傾き（変化量）との差分である。

この磁場の差分ΔＭがセンサの感度よりも小さくなる場合、電流Ｉ₀からの磁場なのか、電流Ｉ₁からの磁場なのか、磁場源を区別できなくなる。

図６のモデルで、ビオ・サバールの法則を用い、電流Ｉ₀と電流Ｉ₁の間隔ｄを無限大にして計算を行うと、センサ１０−１で検出される電流Ｉ₀による磁場の傾きと電流Ｉ₁による磁場の傾きの差分ΔＭは、

で表される。この差分ΔＭがセンサ１０−１の感度よりも小さい場合は、センサ１０−１が無限遠方向に磁場源を区別できない。この意味で、感度限界を、無限遠方向の２つの磁場源からの磁場を分解可能な限界距離と言い換えてもよい。

図７は、ΔＭを磁場源までの距離ｒ［ｍ］の関数としてプロットした図である。センサ１０−０及びセンサ１０−１の感度を水平ラインＡで示す。図７において、

が、空間分割におけるセンサ１０の感度限界を求める条件式となる。図７のシミュレーション結果から、ΔＭは１／ｒ³のオーダーで減少し、センサ１０から比較的近い距離で感度限界になることがわかる。ΔＭがセンサ感度を超える領域を、中心から近い容積部分として用いる。ΔＭがセンサ１０の感度以下となる距離ｒから無限遠までをひとつの分割空間として、中心から遠い容積部分とする。

センサ１０の感度限界を空間分割の指標に用いることで、無限遠までの空間分割を有限個の空間分割で済ませ、計算コストを現実的な範囲に抑えることができる。

また、ΔＭを等間隔にとった場合、対応する横軸の間隔はｒ3のオーダーで広くなり、一つの間隔に含まれるプロット点の数が多くなる。したがって、センサ１０から遠くなるほど空間分割の間隔を広げて間引き処理することで、計算量を少なくすることができる。

さらに、センサ１０から所定の距離内ではもともと磁場源に対するセンサ１０の感度が良いので過度に細かく空間分割しなくても高精度の測定結果が得られるので、所定の距離内では空間分割を等間隔にしてもよい。

上述の可変の空間分割法は、図３のような垂直方向の磁束の差分を測定する軸型グラジオメータだけではなく、センサ１０が平面型グラジオメータを用いたＳＱＵＩＤセンサである場合にも適用可能である。

図８は、第１実施形態の磁場解析方法のフローチャートである。

ステップＳ１１で、センサ１０の感度に応じた空間分割を行う。センサ１０の感度に応じた空間分割は、センサ１０の感度限界の算出と、算出された感度限界に基づく可変の空間分割を含む。センサ１０の感度限界の算出は、図７の条件を満たさなくなる限界距離、すなわち、センサ１０が有する感度において、ある磁場源（電流）からの磁場と、他の磁場源からの磁場を無限遠方向で区別して検出し得る（分解可能な）限界距離ｒ_limの算出である。

感度限界が算出されたならば、感度限界から無限遠までをひとつの分割空間として空間分割を行う。また、感度限界までは距離ｒの３乗のオーダーで径方向の空間分割の間隔を広げる。センサ１０の配置位置から所定の範囲内では径方向に等間隔で分割してもよい。径方向への分割間隔は可変であるが、立体角θは等間隔で分割してもよい。感度限界よりも外側の空間は、立体角θに応じた有限個の空間部分（容積部分）に分割される。

ステップＳ１２で、空間フィルタ法により磁場源の推定を行う。ここで用いる空間フィルタ法は、ステップＳ１１で空間分割した各容積部分の磁場源を推定できればどのような手法でもよい。たとえば、最小ノルムフィルタ法や、ＵＧＭＮ（Unit Gain minimum Norm）空間フィルタ法で磁場源の解析を行う。時刻ｔにおける脳磁信号をＭ個のセンサ１０で検出し、各センサ１０で得られた磁場信号ｂi（ｔ）を要素にもつＭ次元のベクトルをベクトルｂ（ｔ）とする。測定データ、すなわちベクトルｂ（ｔ）に対する関心領域中の空間位置ｒの電流密度への重み行列をｂ（ｔ）に作用させて、位置ｒでの電流密度を算出し、電流密度分布ｊ（ｒ、ｔ）を得る。

ステップＳ１３で、関心空間の中心に近い容積部分からの信号と、中心から遠い容積部分からの信号を分離し、中心から遠い容積部分からの信号を外来ノイズとして除去する。ステップＳ１１で行った空間分割に従って、Ｓ１２で磁場源推定が行われているので、推定結果を用いて中心から遠い容積部分からの信号を外来ノイズとして除去するのは容易である。たとえば、感度限界よりも外側の容積部分からの信号を、外来ノイズとして除去してもよい。

この方法により、少ない計算量で推定精度を維持して磁場源を解析することができる。

図９は、実施形態の磁場解析装置２０のハードウェア構成図である。磁場解析装置２０は、ＣＰＵ（Central Processing Unit：演算装置）２１、ＲＡＭ（Random Access Memory）２２、ＲＯＭ（Read Only Memory）２３、ドライブ２４、入出力インタフェース２５、及びディスプレイ２８を有し、これらがバス２７で接続されている。このような磁場解析装置２０は、汎用コンピュータで実現されてもよい。

ＣＰＵ２１は、磁場解析装置２０の全体を制御する。ＲＡＭ２２は、ＣＰＵ２１のワークエリアとして使用され、磁場解析に用いるデータを一時的に保持する。ＲＯＭ２３は、後述する信号処理プログラムと、基本の入出力プログラムを記憶している。ドライブ２４は、ＨＤＤ（Hard Disk Drive）、ＳＤＤ（Solid State Drive）等の補助記憶装置や、リムーバブルディス等の外部記憶媒体に対して磁場解析データの読み書きを行う。入出力インタフェース２５は、タッチパネル、キーボード、マウス等の操作入力手段と、所定の通信プロトコルによるデータ／ファイル入出力装置と、スピーカ等の出力手段とを含む。データ収録サーバ４２で収集されたセンサ出力は、入出力インタフェース２５により磁場解析装置２０で取得される。ディスプレイ２８は磁場解析結果を画像表示する。

図１０は、磁場解析装置２０の機能ブロック図である。磁場解析装置２０は、データ取得部２０１と、出力制御部２０２と、表示制御部２０３と、演算部２１０を有する。演算部２１０は、感度限界算出部２１１と、空間分割部２１２と、磁場源推定部２１３と、ノイズ除去部２１４を有する。

感度限界算出部２１１は、センサ１０の感度、位置情報、ビオ・サバールの法則等に基づいて導かれる図７の条件式に従って、センサ１０で２つの磁場源からの磁場を無限遠方向に分解可能な限界距離、すなわち感度限界を算出する。

空間分割部２１２は、算出された感度限界にしたがって、有限の空間分割を行う。具体的には、図５のように径方向に感度限界以降の空間領域をひとつの分割空間にまとめる。好ましくは、空間分割部２１２は、感度限界までは、センサ１０の配置位置から遠くなるほど径方向の分割サイズを大きくする可変の空間分割を行う。さらに好ましくは、センサ１０から一定の距離範囲では、径方向に等間隔で空間分割を行う。なお、空間分割のための立体角については等間隔であってもよい。

磁場源推定部２１３は、適切な処理方法、たとえば、最小ノルムフィルタ法や、ＵＧＭＮ空間フィルタ法を用いて、分割された各空間における磁場源を推定し、磁場源の空間分布を推定する。磁場源とは、この例では神経細胞の活動に伴う電流の発生源や、ノイズとなる外部の磁場源である。

ノイズ除去部２１４は、各空間での推定結果を用い、中心から遠い容積部分の磁場源推定値を外来ノイズとして除去する。空間分割部２１２による空間分割にしたがって磁場源が推定されているので、磁場源推定部２１３の演算量が低減され、また、ノイズ除去部２１４によるノイズ除去も容易である。

演算部２１０の各部の機能は、ＣＰＵ２１によって実現することができる。また、データ取得部２０１と、出力制御部２０２は、入出力インタフェース２５により実現される。表示制御部２０３はＣＰＵ２１により実現され、ディスプレイ２８に表示される画像を制御する。

磁場解析装置２０を用いることで、少ない計算コストで精度良く磁場源の推定を行うことができる。磁場解析装置２０は、脳磁計だけでなく、脊髄計や心磁計等、その他の生体磁場の解析に適用することができる。
＜第２実施形態＞
図１１〜図１３を参照して、第２実施形態の磁場解析を説明する。第２実施形態では、第１実施形態での空間分割法を前提とし、演算処理方法を工夫して外来ノイズの除去と磁場源の空間分布の推定を同時に行うことで、計算コストをさらに低減する。

図１１は、第２実施形態の磁場解析方法のフローチャートである。まず、ステップＳ２１で、センサ１０の感度限界、すなわち無限遠方向での空間分解の限界に応じた空間分割を行う。この空間分割は、第１実施形態で説明した空間分割であり、感度限界から無限遠までを径方向への単一の分割空間として扱う空間分割である。これにより、感度限界から無限遠までを、立体角で決まる有限の個数の空間に分け、計算コストを現実的なものとすることができる。また、センサ１０の配置位置から遠くなるにしたがって空間分割を粗くする可変分割として、検出精度の維持と計算量の低減を両立させてもよい。

ステップＳ２２で、分割した各空間の代表点での単位電流素あたりの磁場センサへの影響値を算出する。磁場センサは、たとえば図３に示す軸型グラジオメータとＳＱＵＩＤを用いたセンサ１０である。軸型グラジオメータに替えて平面型グラジオメータを用いてもよいし、軸型と水平型を併用してもよい。

磁場センサへの影響値は、単位電流のある場所からセンサ１０にどのくらいの磁場が影響するかを示す値である。電流素から磁場を求める方法として、Biot-Savartの式、Sarvasの式、境界要素法、有限要素法などがある。実施形態では、一例としてSarvasの式を用いるが、空間領域によって電流素から磁場を求める式を使い分けてもよい。たとえば、遠方の磁場源についてはSarvasの式の仮定が成り立たない場合も多いので、センサ１０の配置位置よりも内部の領域にはSarvasの式を適用し、外部の領域にはBiot-Savartの式を適用してもよい。

センサ１０の数がｍ個、各時刻でのセンサ１０の観測データをｘ（ｔ）、観測データを１回で処理できる長さに分割したときの分割された時間長をＴ、空間分割した各空間の代表点での単位電流素の量をｑ（ｔ）（ただし、立体空間での方向成分が３軸あるので、データ点ｎは代表点の数の３倍になる）とする。分割した各空間の代表点での単位電流素当たりのセンサ１０への影響値は、センサ１０の１個目のコイル１５ａと２個目のコイル（参照コイル）１５ｂを貫通する磁束の差分になり、以下の式（１）で記述される。

ここで、ｌ₁₁−ｌ^ref ₁₁は、センサ１における点１からの磁場の差分（変化量）を示し、ｌ_mn−ｌ^ref _mnは、センサｍにおける点ｎからの磁場の差分（変化量）を示す。式（１）の右辺のｍ×ｎ行列に含まれる各要素または成分が単位電流素当たりのセンサ１０への影響値である。

図１１の処理フローでは、各空間の代表点での単位電流素当たりの各センサ１０への影響値を式（１）から算出しているが、磁場解析時の計算量を減らすため、予め計算を行っておき、計算結果をテーブル等で保存しておいてもよい。この場合、実行中に式（１）を用いた算出処理を省略することができる。

ステップＳ２３で、観測データを指定された時間長で分割する。観測データを１回で処理することができる時間長をＴとすると、次の時間ｔs+1は、式（２）で表される。

一回で処理可能な時間長Ｔは任意に指定できるが、経験的には３〜７秒と指定するのが望ましい。

ステップＳ２４で、最初の時間の観測データを処理対象として設定する。

ステップＳ２５で、各分割空間の代表点のセンサ１０への影響値の集合を辞書と考え、各分割空間の代表点での各時間の電流素の量（強さ）を係数と考えて、式（１）式を用いて、スパース解を求める。一般に、スパース表現は観測信号行列Ｘを、辞書行列Ｄとその係数行列Ｃを用いて、
Ｘ＝ＤＣ＋Ｅ
と表現される。ここでＥは誤差を表わす行列であるが、一般に行列Ｅを各行列要素が平均ゼロの行列とすることができる。式（１）の右辺の左側の行列を辞書行列Ｄ、右側の行列を係数行列Ｃと考える。

スパース解とは、ほとんどの成分の値がゼロで有効値が疎（スパース）である場合の推定値である。解析対象の全空間で見た場合、一般的に磁場源となる電流が存在するのはごく一部であり、磁場源電流（係数行列）がスパースとなる解を求めることで磁場源を推定することができる。スパース解を求める方法としては様々なものが存在するが、一例として、マッチング追跡の手法によりスパースな解を求める。

ステップＳ２６で、指定時間Ｔごとに分割された観測データｘ（ｔ）の処理が全て完了したかを判定する。指定時間で分割された観測データの処理が全て完了していれば（Ｓ２６でＹＥＳ）、処理を終了する。すべての観測データの処理が完了していなければ（Ｓ２６でＮＯ）、ステップＳ２７で、次の時間ｔｓ+1の観測データを処理対象を設定し、すべての処理が完了するまでステップＳ２５〜Ｓ２６を繰り返す。

図１１の例では、観測データの処理を逐次的に行っているが、必ずしも時間順に処理しなくてもよい。並列処理を行う場合は、時刻の異なる観測データ群を同時に処理してもよい。

図１２は、図１１のステップＳ２５の処理の具体例を示すフローチャートである。図１２では、マッチング追跡法を用いる。マッチング追跡法は、残差の低減に寄与する辞書データを順次選択していく貪欲法であり、解の最適性は保証されないが、多くの場合優れた近似を与える。

ステップＳ２５１で、残差Ｒに観測データｘ（ｔ）を設定する。式（１）と同様に、センサ１０の数がｍ個、各時刻でのセンサーの観測データをｘ（ｔ）、観測データを１回で処理できる分割された時間長をＴ、残差をＲとすると、残差Ｒは式（３）で記述される。

ステップＳ２５２で、辞書データの中で、現在の残差から寄与分を引いて求めた新たな残差が最も小さくなるものを選択する。ここではマッチング追跡のアルゴリズムを適用して、式（１）における分割された各空間の代表点での単位電流素当たりの磁場センサーへの影響値（右辺の左側の行列）を辞書データの行列Ｄと考える。行列Ｄは式（４）で記述される。

この辞書データＤのｎ個のベクトルの中で、現在の残差から寄与分を引いて求めた新たな残差が最も小さくなるものを選択する。新たな残差は、式（５）で計算される。

ここで、行列ｄ^* _kは、行列ｄ_kの共役転置行列であり、

は、ｋ番目の辞書データｄｋの共役転置ともとの残差行列Ｒとの行列積を正規化したものである。式（５）は、行列Ｒと寄与分との近さを表わす２乗ノルムである。

辞書データＤのｎ個のベクトルｄ１〜ｄｎの中で、現在の残差から寄与分を引いて求めた新たな残差、すなわち式（５）の値が最も小さくなるｋを選択する。

ステップＳ２５３で、選択された辞書データの行列ｄｋと、正規化された行列積（式（６））を保存する。選択された辞書データの行列ｄｋは位置を表わし、式（６）の正規化された行列積は、電流素の方向と強さを表わすので、これらの演算結果を磁場源の推定で用いるために保存する。

ステップＳ２５４では、式（７）により、現在の残差行列Ｒを、選択した辞書データの寄与分を引き算した新たな残差行列（式（５）の計算結果）に更新する。

現在の残差行列Ｒから選択した辞書データｄｋの寄与分を引き算して、新たな残差を作ることで、残差の低減に寄与する辞書データを順次選択していく。

ステップ２５５で、残差の大きさが指定値以下になったかを判定する。このステップは近似の程度を判定する工程であり、残差が十分に小さくなればうまく近似できていることになる。数式で示せば式（８）で表される。

ここで、εはあらかじめ指定された残差の閾値である。指定する残差の大きさ（閾値）εは、ノイズが大きい環境では大きく、ノイズが小さい環境では小さく設定するのが望ましい。

このようなマッチング追跡法で、残差が一定範囲に収束するまで、寄与分が最も大きい辞書データから順に処理してスパース解を求めることで、外来ノイズの除去と磁場源の空間分布の推定を同時に行うことができる。この処理方法を第１実施形態の空間分割の手法と組み合わせることで、処理時間と演算量をさらに低減することができる。

図１３は、第２実施形態の磁場解析装置２０Ａの機能ブロック図である。第２実施形態の磁場解析装置２０Ａのハードウェア構成は、第１実施形態と同じであり、図９に示されるとおりなので重複する説明を省略する。

磁場解析装置２０Ａは、データ取得部２０１と、出力制御部２０２と、表示制御部２０３と、演算部２１０Ａを有する。演算部２１０Ａは、感度限界算出部２１１と、空間分割部２１２と、磁場源推定部２１３Ａを有する。磁場源推定部２１３Ａは、影響値算出部２２３と、スパース解算出部２２４を含む。

感度限界算出部２１１は、第１実施形態と同様に、センサ１０の感度、位置情報、ビオ・サバールの法則等に基づいて導かれる図７の条件式に従って、センサ１０の感度限界を算出する。

空間分割部２１２は、算出された感度限界にしたがって、有限の空間分割を行う。具体的には、図５のように径方向に感度限界以降の空間領域をひとつの分割空間にまとめて有限個数の分割空間を生成する。好ましくは、空間分割部２１２は、感度限界までは、センサ１０の配置位置から遠くなるほど径方向の分割サイズを大きくする可変の空間分割を行う。さらに好ましくは、センサ１０から一定の距離範囲では、径方向に等間隔で空間分割を行う。なお、空間分割のための立体角については等間隔であってもよい。

磁場源推定部２１３の影響値算出部２２３は、空間分割部２１２により分割された各空間の代表点での単位電流素当たりのセンサ１０への影響値（磁場）を算出する。単位電流素当たりのセンサ１０への影響値は、Biot-Savartの式、Sarvasの式、境界要素法、有限要素法など、適切なアルゴリズムを用いて算出され得る。

スパース解算出部２２４は、データ取得部２０１が取得観測データを１回で処理可能な長さに分割し、分割された時間ごとに、順次または並列で、磁場源の分布を表わすスパース解を求める。スパース解の算出方法として、図１３のようなマッチング追跡法を用いた場合は、外来ノイズの除去と磁場源の推定を同時に行うことができ、処理時間と計算量を短縮することができる。

磁場解析装置２０Ａを用いることで、第１実施形態よりもさらに少ない計算コストで精度良く磁場源の推定を行うことができる。磁場解析装置２０Ａも、脳磁計だけでなく、脊髄計や心磁計等、その他の生体磁場の解析に適用することができる。

なお、第２実施形態では、外来ノイズの除去と磁場源の空間分布の推定を同時に行うことで公知の磁場解析技術よりも精度が高くかつ早い処理速度でノイズ除去が行えるので、センサ１０の配置面よりも外側の磁場解析と外来ノイズの除去に第２実施形態の手法を用い、内部領域については従来の磁場源推定方法を用いることも可能である。

第１実施形態及び第２実施形態を通して、本願の構成と手法は脳磁計だけではなく、脊髄系、心磁計等の任意の生体磁場計測システムや、地磁気の計測と解析にも適用可能である。また、センサ１０の種類は軸型のグラジオメータを用いたセンサに限定されず、平面型のグラジオメータを用いたＳＱＵＩＤセンサにも適用可能である。いずれの場合も、センサ１０の感度限界を考慮し、感度限界から無限遠までを径方向で１個の分割空間とすることで、有限個の空間分割で計算コストを現実的な範囲に抑えることができる。また、センサ１０の配置位置から遠くなるに従って空間分割を粗くすることで、計算コストの低減と検出精度の維持を両立させることができる。

さらに、実施形態の空間分割を用いた信号処理の手法は、生体磁場や地磁気などの解析だけではなく、画像処理等の分野にも適用可能である。その場合は、実施形態の空間分割を用い、センサからの測定データに基づいて信号源の分布を解析することができる。

第１実施形態と第２実施形態において、ＲＯＭ２３に格納された信号処理プログラムをＣＰＵ２１で実行することによって、磁場解析装置２０及び２０Ａの各機能を実現することが可能である。この場合、ＣＰＵ２１は、信号処理プログラムにより、
（ａ）センサから測定データを取得する手順と、
（ｂ）前記センサの感度限界を算出する手順と、
（ｃ）前記感度限界に基づいて径方向に可変の空間分割を行う手順と、
（ｄ）前記空間分割を用いて信号源を推定する手順と、
（ｅ）前記信号源の推定結果を用いてノイズを除去する手順と、
を実行する。

広く信号処理に適用可能であり、特に磁場計測の分野に好適に適用される。

１ 磁場計測システム
３ 測定装置
１０ センサ
２０、２０Ａ 磁場解析装置（信号処理装置）
２１ ＣＰＵ（プロセッサ）
２２ ＲＡＭ（メモリ）
２３ ＲＯＭ（メモリ）
２４ ドライブ
２６ 入出力インタフェース
２８ ディスプレイ
３０ デュワ
４２ データ収録サーバ
２０１ データ取得部２０１
２０２ 出力制御部
２０３ 表示制御部
２１０ 演算部
２１１ 感度限界算出部２１１
２１２ 空間分割部
２１３ 磁場源推定部
２１４ ノイズ除去部
２２３ 影響値算出部
２２４ スパース解算出部

特許第５０１４７８３号

Claims (10)

1. センサから測定データを取得するデータ取得部と、
   前記センサの感度限界を算出する感度限界決定部と、
   前記感度限界に基づいて径方向に可変の空間分割を行う空間分割部と、
   前記空間分割を用いて信号源を推定する推定部と、
   前記信号源の推定結果を用いてノイズを除去するノイズ除去部と、
   を有することを特徴とする信号処理装置。
2. 前記空間分割部は、前記感度限界から無限遠までをひとつの分割空間として径方向の空間分割を行うことを特徴とする請求項１に記載の信号処理装置。
3. 前記空間分割部は、前記センサから遠くなるにしたがって、空間分割を粗くすることを特徴とする請求項１または２に記載の信号処理装置。
4. 前記空間分割部は、前記センサから一定の距離範囲では径方向に等間隔に空間分割を行い、前記一定の距離範囲を超える空間では、前記センサから遠くなるにしたがって空間分割を粗くすることを特徴とする請求項１〜３のいずれか１項に記載の信号処理装置。
5. 前記空間分割部は、前記センサからの距離の３乗にしたがって空間分割を粗くすることを特徴とする請求項３または４に記載の信号処理装置。
6. 前記推定部は、各分割空間における単位電流素が前記センサに与える影響値を算出し、前記影響値の集合を辞書行列とし、前記各分割空間での電流素の強さを係数行列としてスパース解を求めることで前記信号源を推定することを特徴とする請求項１〜５のいずれか１項に記載の信号処理装置。
7. 前記推定部は、前記影響値を算出する際に、前記センサよりも外側の空間と前記センサよりも内側の空間で異なる信号算出法を用いて、前記単位電流素からの影響値を算出することを特徴とする請求項６に記載の信号処理装置。
8. 磁気センサにより磁場を測定する測定装置と、
   前記測定装置で測定されたデータを解析する解析装置と、
   を有し、
   前記解析装置は、
   前記測定装置から測定データを取得するデータ取得部と、
   前記磁気センサの感度限界を算出する感度限界決定部と、
   前記感度限界に基づいて径方向に可変の空間分割を行う空間分割部と、
   前記空間分割を用いて磁場源を推定する推定部と、
   前記磁場源の推定結果を用いてノイズを除去するノイズ除去部と、
   を有することを特徴とする磁場計測システム。
9. センサから測定データを取得し、
   前記センサの感度限界を算出し、
   前記感度限界に基づいて径方向に可変の空間分割を行い、
   前記空間分割を用いて信号源を推定し、
   前記信号源の推定結果を用いてノイズを除去する、
   ことを特徴とする信号処理方法。
10. コンピュータに、
    センサから測定データを取得する手順と、
    前記センサの感度限界を算出する手順と、
    前記感度限界に基づいて径方向に可変の空間分割を行う手順と、
    前記空間分割を用いて信号源を推定する手順と、
    前記信号源の推定結果を用いてノイズを除去する手順と、
    を実行させることを特徴とする信号処理プログラム。

Images