JP4379155B2

JP4379155B2 - 近赤外光による生体光計測方法及び装置

Info

Publication number: JP4379155B2
Application number: JP2004058211A
Authority: JP
Inventors: 理河野
Original assignee: Shimadzu Corp
Current assignee: Shimadzu Corp
Priority date: 2004-03-02
Filing date: 2004-03-02
Publication date: 2009-12-09
Anticipated expiration: 2024-03-02
Also published as: JP2005245636A

Description

本発明は被検体の複数の位置で近赤外光による測定を行なう生体光計測方法と装置に関し、例えばそのようにして得た測定結果を２次元的な画像として表示する近赤外光イメージング装置などに関するものである。

生体光計測方法としては、近赤外光イメージング装置により人体の頭部の複数個所で近赤外光による脳からの透過散乱光を検出し２次元画像として表示して脳の活動状態を研究することが行なわれている。

その際、脳の血液からの信号とそれ以外の信号、例えば皮膚血流からの信号とを峻別する必要がある。
その１つの方法は、研究者の視覚的評価で、タスクに同期している信号のみを脳賦活信号として扱うことである。

他の方法は、視覚的評価によって、ｏｘｙ−Ｈｂ（オキシヘモグロビン）とｄｅｏｘｙ−Ｈｂ（デオキシヘモグロビン）の挙動から、脳内の生理学的変化と考えにくい信号はアーチファクトとして処理する方法である。ランダムノイズに関しては、測定の繰返し回数を多くし、積算処理によって除去していた。

さらに他の方法として、複数の計測位置においてそれぞれ得られたヘモグロビン変化信号から、その統計的有意性を判定する評価量を求め、この評価量又は評価量に基づいて判定した結果をモニタに表示することが提案されている（特許文献１参照。）。評価量は、例えば、ヘモグロビン変化信号からＴ値（検定統計量）及び／又はＰ値（有意確率）を求め、Ｔ値又はＰ値から検定関数を定義し、各計測位置の判定関数を合算したものを評価量として用いている。
特開２００３−２６５４４２号公報 L.Molgedey and H.G.Shuster, Separation of a mixture of independant signals using time delayed correlations, Phys. Rev. Lett., Vol.72,No.23,pp3634-3637, 1994 J.CARDOSO and A.SOULOUMIAC, JACOBI ANGLES FOR SIMULTANEOUS DIAGONALIZATION, SIAM J. MATRIX ANAL. APPL., Vol.17,No.1,pp.161-164, 1996

近赤外光イメージング装置で測定された脳計測データは、脳賦活に伴う血流に基づく信号の他、皮膚血流や血圧変動、脈動・呼吸などの変化に基づく信号も重畳されている可能性がある。特に、歩行などの運動をタスクとして与えた場合における脳活動を測定するとき、心拍数が変化し、皮膚血流もタスクに同期している可能性がある。よって、これらの脳活動とは関係しない信号であっても、タスクと同期する信号は、タスクとの同期に基づいて信号を峻別する従来の方法では除去できない。
本発明は不必要な信号を除去して観測信号を求める方法と装置を提供することを目的とするものである。

本発明の生体光計測方法は、被検体の複数の照射位置からの近赤外光照射により被検体内を透過し又は被検体内で散乱もしくは反射して被検体から出てきた光信号強度を前記照射位置に対応した複数の検出位置で検出した信号を処理する生体光計測方法であって、その信号処理において独立成分分析手法を用い、観測された信号から独立な信号成分と混合行列を求めた後、混合行列のうち、不必要成分に該当する列ベクトルとランダムノイズに該当する列ベクトルの一方又は両方に０を代入し、その後、混合行列と独立成分信号の積を計算することによって、不必要成分に基づく信号とトランダムノイズに基づく信号の一方又は両方が除去された元の観測信号を復元するものである。

すなわち、その信号処理は以下のステップを含むデータ処理を行なう。
（ａ）前記複数の検出位置での光信号強度の検出信号をもとにして独立成分分析によって混合行列と独立成分信号を得るステップ、
（ｂ）得られた混合行列のうち、不必要成分に該当する列ベクトルに０を代入するステップ、及び
（ｃ）上のステップ（ｂ）で得られた混合行列とステップ（ａ）で得られた独立成分信号との積を計算して観測信号を求めるステップ。

独立成分分析とは、下記の式で示される観測信号から、混合行列と独立成分信号を求める手法のことである。

混合行列が既知であれば、式（１）の混合行列の逆行列を求めることによって、独立成分信号は容易に求められるが、一般に、混合行列も未知であるため、それらを求めるために、様々なアルゴリズムが提案されている（非特許文献１，２参照。）。

混合行列の列ベクトルは検出位置における特定の独立成分の重みを表わしている。例えば、脳活動を測定する場合について説明すると、一般に、皮膚血流による空間マップは、脳活動によるマップに比べて、局所的ではない。よって、皮膚血流の独立成分に対応する列ベクトルの標準偏差係数の絶対値は、脳活動の独立成分に対応する列ベクトルの標準偏差係数の絶対値と比較して小さくなる。そこで、混合行列の列ベクトルの標準偏差係数の絶対値αを計算することによって、皮膚血流など局所的でないものに相当する独立成分を推測することが可能となる。

ｋ列目の標準偏差係数の絶対値αｋは以下の式で与えられる。

ここでσｋは標準偏差であり、

と表わすことができる。σｋはまた、

と表わしてもよい。上にバーのついた記号は平均値を表わす。

列ベクトルの標準偏差係数の絶対値αを利用する１つの方法は、閾値との大小関係に基づいて、０を代入する列ベクトルを決定するものである。
そこで、ある閾値ｃを設け、標準偏差係数の絶対値αがｃ以下となった場合の成分が皮膚血流成分など局所的でないものとする。その場合、混合行列から標準偏差係数の絶対値αがｃ以下となった列の列ベクトルを０として、混合行列が変更される。０にされる列ベクトルの数は１つに限らず、複数になることもありうる。なお、閾値との比較における「以下」の表現は厳密なものではなく、「未満」も含む意味で使用している。同様に「以上」の表現も厳密なものではなく、「それより大きい」も含む意味で使用している。

混合行列の列ベクトルの標準偏差係数の絶対値αを利用する他の方法は、その絶対値αが最小となる列ベクトルに０を代入する方法である。
ステップ（ｂ）において０を代入する列ベクトルを決定する他の方法として、不必要成分のみについての光信号強度の検出信号の波形を別途測定し、ステップ（ａ）で得られた独立成分信号の波形のうちで不必要成分の検出信号波形にもっとも近い波形をもつ独立成分に該当する列ベクトルに０を代入することもできる。
このようにして、混合行列の特定の列ベクトルに０が代入された後、式（１）によって、観測信号が計算される。これにより、脳活動に基づく信号から皮膚血流など局所的でないものに基づく信号が除去される。

一方、ランダムノイズは標準偏差係数の絶対値が極めて大きい。そこで、ある閾値ｄを儲け、標準偏差係数の絶対値αｉがｄ以上という条件を設定し、その条件を満足した場合に、それに対応する混合行列の列ベクトルを０にすることによって、ランダムノイズを除くことができる。
ステップ（ｂ）において特定の混合行列の列ベクトルを０にする操作は、局所的でないものに基づく信号を除去するための操作と、ランダムノイズを除去するための操作とがある。これらの操作のいずれかを実行してもよく、両方を同時に実行してもよい。

本発明の生体光計測装置は、近赤外光を被検体の複数の照射位置に照射する光照射部と、被検体内を透過し又は被検体内で散乱もしくは反射して被検体から出てきた光信号強度を前記照射位置に対応した複数の検出位置で検出する光検出部と、前記光検出部が検出した光信号強度を処理するデータ処理装置とを備えたものであり、そのデータ処理装置として本発明の生体光計測方法におけるステップ（ａ）から（ｃ）を含むデータ処理を行なうようにしたものである。

本発明によれば、不必要成分に基づく信号とランダムノイズに基づく信号の一方又は両方が除去された元の観測信号を復元することができる。
不必要成分に基づく信号を除去した場合は、例えば脳活動の測定においては、皮膚血流などの成分を除去して脳活動による信号成分の挙動を精度よく検討できる。更に、マッピング画像を作成した場合は、脳賦活部位を精度良く特定することが可能となる。
また、ランダムノイズ成分を除去した場合は、Ｓ／Ｎ比（信号対ノイズ比）が向上する。

図１は本発明が適用されるマルチチャンネル生体光計測装置の一例を概略的に表わしたものであり、被検体に光を照射する送光端子４と被検体中を透過及び／又は散乱もしくは反射した後に外部に放出される光を受光する受光端子６とが配置面上にそれぞれ複数個ずつ配置されたプローブ１２と、プローブ１２の受光端子６が受光した光をもとに解析するデータ処理装置１８とを備えている。

１４は各送光端子４に測定用の光を供給する送光端子ごとの光源であり、ＬＥＤ（発光ダイオード）やＬＤ（レーザダイオード）などの発光素子を用いることができる。各光源１４からの光はそれぞれ光ファイバなどの導光路を介してそれぞれの送光端子４に導かれる。この場合、送光端子４は導光路の光送光端面となる。

１６は各受光端子６が受光した光を検出する検出器であり、フォトダイオードやフォトトランジスタなどの受光素子のほか、光電子増倍管などを使用することができる。各受光端子６からの光は光ファイバなどの導光路を介して各検出器１６に導かれる。

光源１４と検出器１６は、データを取得しようとする位置の送光端子４につながる光源１４を点灯させ、その送光端子４から照射された光を受光しようとする受光端子６につながる検出器１６を動作させるというように各検出器１６が所定の位置のデータを検出するように、動作のタイミングが制御される。光源１４と検出器１６のそのような動作もデータ処理装置１８により制御される。２０はデータ処理装置１８で処理されたデータを表示する表示装置である。

一実施例の動作を図２のフローチャートを参照して説明する。
まず初めに、空間的に異なった位置において、複数個の信号が観測される。
次に、これらの信号が独立成分分析にかけられ、混合行列と独立成分信号が求められる。
求められた混合行列の各列ベクトルの標準偏差係数の絶対値α１〜αｎが求められる。

それらの標準偏差係数の絶対値α１〜αｎのうち、不要成分に該当する信号を除去するために、閾値ｃ以下（未満でもよい。）となった列ベクトルを０とする。閾値ｃは不要成分に該当するものとして経験的に決定する。ここで、不要成分に該当する信号を除去するために、混合行列の特定の列ベクトルを０とする他の方法として、標準偏差係数の絶対値α１〜αｎのうちの最小の列ベクトルを０としてもよく、不必要成分のみについての光信号強度の検出信号の波形を別途測定し、独立成分信号の波形のうちで不必要成分の検出信号波形にもっとも近い波形をもつ独立成分に該当する列ベクトルを０としてもよい。

さらに、ランダムノイズに基づく信号を除去するために、標準偏差係数の絶対値αが閾値ｄ以上の列ベクトルに０を代入してもよい。
混合行列の特定の１又は複数の列ベクトルを０にする操作の後、（１）式によって観測信号を計算すると、皮膚血流などの不要成分に基づく信号とランダムノイズの一方又は両方が除かれた観測信号が得られる。

より具体的に示すために、人の頭部の測定のにおいて、不必要成分のみについての光信号強度の検出信号の波形を別途測定し、独立成分信号の波形のうちで不必要成分の検出信号波形にもっとも近い波形をもつ独立成分に該当する列ベクトルを０としてデータ処理した場合の結果を示す。

図３は不必要成分としての皮膚血流をレーザー血流計により測定した結果の波形図であり、時間に対する皮膚血流の変化を示している。横軸は時間（秒）である。
この測定と同時に頭部の複数個所でオキシヘモグロビンの測定を行ない、独立成分分析を行なった。そのときの独立成分信号を図４に示す。図４中の数字は独立成分信号番号を表わしている。

皮膚血流波形に最も近いのは独立成分信号番号２の波形である。それらの波形を図５に並べて示す。図５の横軸は時間（秒）である。この２つの波形の複数の時間での信号強度の相関関係を図６に示す。横軸Ｘは皮膚血流波形データ、縦軸Ｙは独立成分信号番号２の波形データである。両者は強い相関があり、その相関係数ｒは０．７２４であった。このことはまた、独立成分信号の重み成分（混合行列の列ベクトル成分）の空間マッピングによっても窺うことができ、独立成分信号番号２の空間マッピングは局在化の程度が低くなる。このことから、信号番号２で表わされる独立成分は皮膚血流であると推定される。

図７は３回目のタスク期間後半（３７５秒〜４２０秒での積算）結果の空間マッピングを示したものである。（ｋ）として示したのは本発明の処理を施さなかった場合、（ｌ）として示したのは同じデータについて混合行列で独立成分信号番号２の列ベクトルに０を代入して観測信号を計算した結果を空間マッピングとして示したものである。混合行列で特定の列ベクトルに０を代入するデータ処理を施すことによって、皮膚血流波形に基づく信号が除去され、脳の活動箇所を明瞭に表示することができることを示している。

以下に、独立成分分析について説明する。
［ＩＣＡ（独立成分分析）問題の設定］
ｎチャンネルのｐサンプル観測信号
ｘ（ｔ）＝（ｘ₁（ｔ），ｘ₂（ｔ），・・，ｘ_n（ｔ））^T：（ｔ＝１，２，・・ｐ）
は、ｎ個の独立な信号源ｓ（ｔ）＝（ｓ₁（ｔ），ｓ₂（ｔ），・・，ｓ_n（ｔ））^Tと線形作用素Ａを用いて、
ｘ（ｔ）＝Ａｓ（ｔ）（３）
と表わせると仮定する。このとき、ＩＣＡ問題とは
ｙ（ｔ）＝Ｗｘ（ｔ）（４）
として、観測信号ｘ（ｔ）からｙ（ｔ）が相互に独立であるという条件のみを利用して、独立成分ｙ（ｔ）と変換行列Ｗを求める問題となる。実際には、Ｗ＝Ａ^-1がわかっておらず、Ｗは要素の順番の入れ違いとその大きさの任意性が残るため、適当な置換行列Ｐと対角行列Ｄを用いて、
ＷＡ＝ＰＤ（５）

を満たすＷを推定する問題となり、時間ｔでのＷの候補をＷ_tとして、ｙ_i＝Ｗ_tｘ_iを観測するたびに、これをＷ_t+1＝Ｗ_t＋ΔＷ_tに変えていくアルゴリズムとして解くことになる。

［ＩＣＡの様々な提案手法］
ＩＣＡは、信号源の独立性の基準としての評価関数として、Kullback-Liebler情報量、最尤法、エントロピー最小化、情報量最大化、高次モーメント、キュムラント、相互相関最小化などを用いる方法がある。また、アルゴリズムとして、勾配法、不動点法、Jacobi法などの方法があり、その評価関数とアルゴリズムによって様々な手法が提案されている（Jutten and Heraut(1991)、Bell and Sejowski(1995)、Amari at al(1996)、Cardoso and Laheld(1996)、Girolami and Fyfe(1996)、Molgedey and Schuster(1994)、Hyvarinen and Oja(1997)、Comon(1994)、Cardoso and Souloumiac(1993)、Matsuoka et al.(1995)）。

（Molgedey and Schusterによる手法）
ここでは、ＩＣＡの例として、Molgedey and Schusterによる手法について説明する。信号ｓ（ｔ）の相互相関関数がゼロであるとすると、（３）、（４）、（５）式より、ｙ（ｔ）の相関関数は、
＜ｙ（ｔ）ｙ（ｔ＋τ）^T＞＝Ｗ＜ｘ（ｔ）ｘ（ｔ＋τ）^T＞Ｗ^T
＝ＷＡ＜ｓ（ｔ）ｓ（ｔ＋τ）^T＞Ａ^TＷ^T
＝ＷＡＲＡ^TＷ^T
＝ＰＤＲＤ^TＰ^T （６）
と表せる。Molgedey and Schusterは、τ＝０とτ＝τ₀の両方において相互相関係数がゼロであることを要求している。ここで、ｓ（ｔ）のいくつかの相互相関関数がゼロになるとすると（６）式の非対角要素の２乗和を表わす関数ｆ（Ｗ）

を最小にする問題を解くことになる。

（７）式の最小化問題を解く方法として、”sphering”、”rotation”と呼ばれる手法を用いる。すなわち、sphering変換行列を

rotation変換行列をＣとして、

を最小にする問題を解く。

Spheringは、適当な変換行列により観測を直交化する処理であり、主成分分析（Principle Component Analysis：ＰＣＡ）の技術を利用できる。観測信号の相関行列である正定行列Ｖは、適当な直交行列Ｓと対角行列Λにより、

と表わすことができる。このとき、Ｓは正規直交基底（観測信号の相関行列の固有ベクトル）を列ベクトルとする行列、Λは対応する基底における分散（観測信号の相関行列の固有値）を対角要素とする行列として求めることができる。この手法はＰＣＡと同一であり、Ｓの列ベクトルは主成分ベクトルと呼ばれ、Λの各要素は対応する主成分の分散を表わす。

ＳとΛを利用して、観測信号を直交化する仕組みは以下の通りである。
今、行列Ｖの平方根を、

と表わすことにする。ただし、

は対角行列Λ^-1の各成分の平方根を要素とする対角行列である。このとき、

により、観測信号ｘ（ｔ）を、

と変換する。すると、ｘ’（ｔ）の相関関数は、

となり（Ｉは単位行列）、ｘ’（ｔ）は直交化している。

Rotationは、直交化されたｘ’（ｔ）を回転し対角化（独立化）する処理である。これは、適当な回転行列Ｃを作用させ、ｘ’（ｔ）の相互相関をゼロにする。すなわち、直交化された相関信号の軸を回転させ、Cross Talkをなくす処理である。このためには、いくつかの時間差における相関行列の非対角成分の和を最小にする問題を解くことになる。すなわち、

を最小にする回転行列Ｃを求める問題を解く。

この問題は、同時対角化（Simultaneous Diagonalization）問題と呼ばれ、例えばCardoso and SouloumiacのJacobi法的アルゴリズムで解くことができる。Cardoso and Souloumiacのアルゴリズムでは、Ｃを２行２列の行列に分解し、

と求める。ただし、ｘ、ｙは下式で示す行列Ｇの最も大きい固有値に対応する固有ベクトルであり、

である。

以上の２つの処理を経て、（４）式の変換行列Ｗは、

と求められ、独立成分ｙ（ｔ）は、

と求めることができる。

本発明は人の頭部にフローブを装着して脳賦活状態を２次元的な画像として表示する近赤外光イメージング装置などに利用することができる。

一実施例の生体光計測装置を概略的に示すブロック図である。一実施例の動作を示すフローチャート図である。レーザー血流計による皮膚血流の変化を示す波形図である。オキシヘモグロビンの独立成分信号を示す波形図である。皮膚血流波形と特定の独立成分信号の波形とを比較して示す波形図である。皮膚血流波形データと特定の独立成分信号の波形データとの相関関係を示すグラフである。観測信号の空間マッピングを示す図であり、（ｋ）は本発明のデータ処理を施さなかった場合、（ｌ）は本発明により皮膚血流と推定される独立成分信号を除去した場合である。

符号の説明

４送光端子
６受光端子
１２プローブ
１４光源
１６検出器
１８データ処理装置
２０表示装置

Claims

被検体の複数の照射位置からの近赤外光照射により被検体内を透過し又は被検体内で散乱もしくは反射して被検体から出てきた光信号強度を前記照射位置に対応した複数の検出位置で検出した信号を処理する生体光計測方法において、
前記信号を処理する方法は、以下のステップを含むデータ処理により不必要成分に基づく信号を除去した観測信号を求める方法。
（ａ）前記複数の検出位置での光信号強度の検出信号をもとにして独立成分分析によって混合行列と独立成分信号を得るステップ、
（ｂ）得られた混合行列のうち、不必要成分に該当する列ベクトルに０を代入するステップ、及び
（ｃ）上のステップ（ｂ）で得られた混合行列とステップ（ａ）で得られた独立成分信号との積を計算して観測信号を求めるステップ。
前記ステップ（ｂ）において、混合行列の列ベクトルの標準偏差係数の絶対値αを求め、その絶対値αが閾値ｃ以下の列ベクトルに０を代入する請求項１に記載の生体光計測方法。
前記ステップ（ｂ）において、混合行列の列ベクトルの標準偏差係数の絶対値αを求め、その絶対値αが最小となる列ベクトルに０を代入する請求項１に記載の生体光計測方法。
不必要成分のみについての光信号強度の検出信号の波形を別途測定し、前記ステップ（ｂ）においてステップ（ａ）で得られた独立成分信号の波形のうちで不必要成分の検出信号波形にもっとも近い波形をもつ独立成分に該当する列ベクトルに０を代入する請求項１に記載の生体光計測方法。
前記ステップ（ｂ）において、混合行列の列ベクトルの標準偏差係数の絶対値αを求め、その絶対値αが閾値ｄ以上の列ベクトルに０を代入してランダムノイズに基づく信号を除去する請求項１から４のいずれかに記載の生体光計測方法。
近赤外光を被検体の複数の照射位置に照射する光照射部と、
被検体内を透過し又は被検体内で散乱もしくは反射して被検体から出てきた光信号強度を前記照射位置に対応した複数の検出位置で検出する光検出部と、
前記光検出部が検出した光信号強度を請求項１から５のいずれかに記載のステップ（ａ）から（ｃ）を含んでデータ処理するデータ処理装置と、
を備えた生体光計測装置。