JP6501432B2 - 血流動態画像化診断装置及び診断方法 - Google Patents

Description

本発明は、血球を有する生体組織にレーザー光を照射し、その生体組織から反射されたスペックル信号に基づき、血流速度を測定し画像化するに際して、関心領域毎に血流波形の拍動（ビート）の強弱を数値化して画像診断する血流動態画像化診断装置及び診断方法に関する。

従来の血流速度測定装置は、眼底や皮膚など血球を有する生体組織にレーザー光を照射して、その血球からの反射光が干渉した結果形成されるランダムな斑点模様の画像、いわゆるスペックル画像が、固体撮像装置（ＣＣＤやＣＭＯＳ）等のイメージセンサー上に導かれる。このスペックル画像は、連続的に所定時間間隔で多数枚取り込んで、記憶される。その記憶された多数の画像は、所定枚数選択して、各画像の各画素における出力の時間変動速度を反映した値を算出し、この値から血球の速度（血流速度）を算出する。この種の血流速度測定装置では、各画素の出力変動速度の値が血球の移動速度に対応するので、この算出された各画素の出力変動値に基づき、生体組織での血流分布を二次元画像（血流マップ）として、モニター画面上にカラー表示することができる。実際に観測される血流マップは、毎秒３０フレーム程度で算出される一連の血流マップ（以下、元マップと呼ぶ）で構成され、動画として表示することもできるので、眼底や皮膚の血行動態を観測する装置として実用化されている（特許文献１〜７参照）。

また本発明者らは、数秒間の血流測定で得られた一連の元マップに対して、心臓の拍動（心拍）による周期性を利用し、元マップを各心拍の先頭から対応する元マップ同士を平均化し、一心拍の平均血流マップ（以下、心拍マップと呼ぶ）を作成し、一心拍における拍動内の血流変化を観測視野内の各部位において解析し、動脈性の鋭い立ち上がり波形を有する部位と、静脈性の緩やかに上下する波形を有する部位を区別できる数値、即ち歪度等を導入し、血流マップ上に動脈性の拍動部分と静脈性の拍動部分を表示することができる血流速度画像化装置も提案した（特許文献８参照）。

さらに本発明者らは、従来の装置に新たな血流画像診断機能を付加し、演算部に一心拍以上の複数の血流マップデータから、平均化処理を行い、時間情報を持たないすべての元マップを平均化した血流マップ（以下、合成マップと呼ぶ）を作成し、この合成マップから生体組織の観察領域の表層血管内の血流と、その周囲の背景血流を分離し、表示部の合成マップ上にそれぞれの血流が区別して表示される機能を加え、分離された各部位の血流波形を特徴付ける種々の変数を定義し、これらを比較して臨床診断に利用する手法も提案している（特許文献８参照）。以下の文中では、このような機能が付加された装置を、血流動態画像化診断装置と呼ぶことにする。

以下、このような従来の血流動態画像化診断装置について、図１５〜図２４を参照して説明する。図１５は、従来技術に基づき構成した血流動態画像化診断装置の概略全体構成図である（特許文献８参照）。レーザー光照射系２は、生体組織１（例えば、足背）にハーフミラー３を介してレーザー光を照射する。受光部４は、レーザー反射光を固体撮像装置に結像する。タイミングパルス６に基づき駆動される受光部４は、受光レンズにより結像された生体組織像を、タイミングパルス６に基づいて電気信号に変換し、フレーム蓄積方式で信号電荷を読み出して、映像信号として増幅して出力する。

画像取込部５は、出力された映像信号に対して、連続的に複数の画像を取り込む。画像記憶部７は、この取り込まれた画像を記憶する。演算部８は、記憶された複数画像の対応する各画素の出力信号の時間変化から、生体組織内の血流速度を演算する。表示部９では演算部８から得られた血流マップの情報に加え、波形数値等の他の情報を重畳させた二次元分布情報をマップ表示できる。

図１６は、従来技術による演算部８の詳細を示す図である。血流解析部１１は、画像記憶部７から取り込まれた画像を解析して、生体組織での血流分布を二次元画像表示した血流マップを作成する。血流マップは、取り込んだ複数枚の画像上の各画素において光強度の変化率を演算し、その変化率をマップ化したものであり、血球の動きを反映する。血流マップ上の各画素は血流値を有し、血球の動きが速い領域は数値が高く、逆に血球の動きがゆっくりとした場合には低値になる。上述したように、画像記憶部７は、一心拍または二心拍以上の所定時間間隔（例示の場合は、１／３０秒間隔）で、連続的に複数の画像を取り込むが、この記憶された複数画像の対応する各画素の出力信号の時間変化から、生体内の血流速度を演算する。この演算した血流値を、図１７において、血流の時間変化として例示している。血流値は、通常脈動するが、図中上方の矢印ａ（４個）は、各々の拍動における最低血流マップの時間的な位置を示している。

心拍解析部１２では、血流解析部１１からの血流マップを解析する。即ち、心拍解析部１２は、時系列に並んだ血流マップ毎に、関心領域内の平均化した血流値に対して、各拍動ごとに血流値が最低となる血流マップを探索して、各心拍の最低血流マップを決定する。この最低となる画像を、図１８において、先頭マップｅとして表示している。決定された血流マップは、心拍マップ作成部１３に送られる。

図１７は心拍検知に成功した場合の血流の時間変化を示している。例示の場合、三心拍の拍動ｂ、ｃ、ｄが観察される。図１８に示す一連のマップｆは、拍動ｂの時間における複数の血流マップであり、一連のマップｇは、拍動ｃに、また一連のマップｈは、拍動ｄに対応する複数の血流マップである。図１８は、拍動に対応する複数の血流マップから合成マップを作成する方法を説明する図である。

心拍マップ作成部１３は、最初に、一つまたは二つ以上の拍動の先頭マップｅ同士を各画素について平均化し、平均化した血流マップｉを作成する。その後、心拍マップ作成部１３は、順次各拍動において２番目、３番目の血流マップｊ、ｋ、・・・同士を平均化し、拍動の終わりまで平均化した血流マップｍを作成する。出来上がった心拍の先頭から時間軸に並べた平均化した複数の血流マップｎを、心拍マップと呼ぶ。さらに心拍マップ作成部１３で作成された心拍マップｎは、合成マップ作成部１４及び波形数値作成部１５に送られる。合成マップ作成部１４では、心拍マップ作成部１３で作成された心拍マップｎを平均化して、一枚の合成マップ（図２０）を作成する。

心拍マップ上の平均血流値を時系列順に並べると、図１９に示すような拍動による血流変化が確認できる。図１９は、一心拍における血流の時間変化を示すグラフである。この血流変化を波形数値計算の基とする。波形数値作成部１５は、この拍動による血流変化から、関心領域内の歪度などの血流動態を数値化した波形数値を作成する。波形数値は、波形マップ作成部１６と共に、後述の表層・背景分離部１７にも送られる。波形マップ作成部１６は、波形数値を画像上に配置してマップ化した波形マップを作成する。

波形数値は、血流が多い少ないというような量としての指標値とは異なり、小動脈の血管から末梢血管までの血管抵抗を表している。数値化することにより個体間で比較することができ、血管抵抗性に代表される血流の質ともいうべき流れの状態を比較できる指標である。また波形数値から血管が動脈性の拍動か、静脈性の拍動かなど識別することができ、マップ内の動静脈の血管の判別にも有用であることがわかっている。

表層・背景分離部１７では、作成した合成マップ（図２０）を基に、関心領域において血流の数値、例えば関心領域内の血流値のヒストグラムを算出する。この血流分布のヒストグラムから、表層と背景領域を分離する血流値の閾値を決定し、基準とした閾値を境に、表層と背景を分離する。また、波形数値作成部１５で作成した歪度などの波形数値を基に、表層血流と背景血流に分離する。例えば、大きな血管があるとした場合、その周囲には組織を養っている末梢血管があり、これらを組織血管と呼んでいる。大きな血管は手前側の表層にあるので、これを表層血流と称し、また、組織血管は大きな血管よりもさらに奥まで血管が走っており、背景までも含んだ血管でもあるので、これを背景血流と称している。

分離マップ作成部１８は、分離した表層と背景をそれぞれ二値化した分離マップを作成する。例えば分離マップ表層にある大きな血管がある領域を白色に、また背景に存在する組織血流の領域を黒色にした分離マップを作成する。分離マップは例えば、大きな血管を白色にし、背景血流を黒にした二値化したマップであるので、血流値を数値として比較できる情報を持っていない。合成マップ作成部１４で作成される合成マップは、マップを概観し、血管の位置を特定するために必要であり、血流値を知るためにも必要である。さらに、分離マップを作成することにより、血流動態の異なる血管系統領域を分離マップで識別し、表示をより見易くすることができる。

このように、心拍検出が成功した場合は、合成マップ・付加情報マップ重畳部１９では、合成マップ作成部１４で作成した合成マップを下層にし、波形マップ作成部１６で作成した波形マップ、及び分離マップ作成部１８で作成した分離マップを上層に重畳する。図１５の表示部９は、合成マップ・付加情報マップ重畳部１９で作成した合成マップの血流の分布に一致した、関心領域の分離マップ、及び波形マップを重畳したマップを表示する。図２１は、心拍検出が成功した場合に、重畳表示されるマップを例示している。図２１の円ｐ内では、分離マップを重畳して表示している。円ｐ内の白い領域が表層の血流領域を示しており、黒い領域が背景の血流領域を示している。図中のバーグラフｑ，ｒは、波形数値として歪度（以下Skewと呼ぶ）をグラフ表示すると共に、数値でも示している。ｑは表層血流の動態指標であるSkewを、ｒは背景血流の動態指標Skewを表している。Skewは歪度（Skewness）に定数を乗した数値である。

上述のように、心拍解析部１２では、血流解析部１１からの血流マップを解析して、血流値が最低となる血流マップを探索して、各心拍の最低血流マップを決定するが、このとき、最低血流マップを決定することができなかった場合、時系列に並んだ血流マップの全てが、心拍検出できない場合として、合成マップ作成部１４に送られる。

図２２は、心拍検知がうまく行えなかった場合を説明する図である。図２２において、測定時間4secの間の血流値を、血流の時間変化ｓとして例示している。図示の血流値は、大きくは脈動していない場合を例示している。合成マップ作成部１４では、血流解析部１１で作成した時系列に並んだ複数枚の血流マップｔを、先頭のマップｕから最後のマップｖまでを平均化し、図２３に示す一枚の血流マップを作成する。これを合成マップと呼ぶ。任意の時間に開始した最初の測定による血流マップが、先頭マップとなり、かつ、所定時間後に測定を終了した最後の測定が、最後のマップとなる。最初のマップが拍動の途中でも構わず、最初のマップから最後のマップまで平均化されることになる。

表層・背景分離部１７は、心拍検知がうまく行えなかった場合に波形数値を利用することができないが、合成マップ上の関心領域内の血流値を基にヒストグラムを得て、この分布から分離することが可能なので、合成マップのみでも分離マップを作成可能である。図２４は、心拍検知がうまく行えなかった場合に重畳表示されるマップを例示している。図中に、波形数値等の何らの表示も無い。

心拍検知が失敗した場合には、合成マップ及び分離マップのみを作成して、合成マップ・付加情報マップ重畳部１９で重畳する。表示部９に提供できる情報は、血流の速い領域と遅い領域を視認できるマップのみになる。心拍検出できないデータには、血流動態を指標化したマップを作成できず、動態指標をマップに表現できない。拍動に合わせて変化する複数の血流マップを観察できるが、客観的な指標を持たないため、個体間の比較をこのマップのままでは行えない。皮膚のような末梢では、拍動が顕著にみられないため、心拍検知に失敗しやすく、波形マップを作成できないため、血流動態を比較できないことがあった。波形数値は、拍動の周波数の強度に注目した動態指標ではないので、拍動による血流変化が少ないまったりとした血流変動では拍動の強度比較が行えなかった。

合成マップ自体は、時間変化の情報を含んでおらず、合成マップのみを基に分離マップを作成しても、拍動の強さ等の情報は表現できない。また時間情報を持っている波形マップは、心拍解析部１２による心拍検知を必要としており、心拍検知に失敗してしまうデータでは、血流動態指標を表現している波形マップを得ることができなかった。このように、従来技術では心拍検出に失敗すると、血流動態に関する情報は平均血流値のみになり、血流の質のような動態にかかわる情報は失われるので、血流を測定したデータすべてに対して、十分な機能提供ができているというわけではなかった。

特公平５−２８１３３号公報 特公平５−２８１３４号公報 特開平４−２４２６２８号公報 特開平８−１１２２６２号公報 特開２００３−１６４４３１号公報 特開２００３−１８０６４１号公報 国際公開第２０１４／１７５１５４パンフレット 国際公開第２００８／６９０６２パンフレット

糖尿病の患者では、病態が深刻になると足病になり足先が壊死してくる。この場合、足先の特定の部位は部分的に拍動が一定せず、血流動態が健常部位とは異なってくる。足先のまだ健常な部位と壊死しかかっている拍動のない領域の鑑別は、健常な領域を切断せずに残したいとする患者およびドクターにとって重要な情報となりうるので、非接触で簡便に二次元測定できる血流動態画像化診断装置が求められていた。

図１５及び図１６を参照して説明した従来の診断機器は、血流の周期的な拍動性を利用し、時系列の連続する血流マップを基に、心拍解析と呼ぶ処理を必要としている。心拍検出できる場合は、心拍マップから血流波形を抽出し、関心領域内の歪度などの血流動態を波形数値として数値化し、マップ化して臨床上診断に有益な情報を提供することができる。

しかし心拍検出できない場合には、時系列の血流マップから拍動を無視した平均化した血流マップを一枚作成し、マップに視覚化するのみである。血流の速い領域、及び遅い領域を視覚化し、関心領域内の平均血流値を比較することができるのみであった。
炎症が発生した個所では、血流が通常よりも多くあることは観察されながらも、観察領域の一部分において、どの程度の血管抵抗性があり、窮屈な流れを有しているかということまでは、合成マップを観察するだけではわからなかった。これまでの波形数値は、心拍抽出が成功したデータしか出力できず、すべての血流データに適用できるというわけではなかった。

特に皮膚血流の場合は、血流量はある程度あるものの、末梢側の血流であるため心拍に同期した規則的な血流変化が減少する。この様態は末梢側に向かうほど顕著で、拍動と思しき周期的な血流変化は非常に確認しづらくなる。また末梢血流には拍動に同期しないもっとゆっくりとした血流変動もあるため、心拍抽出処理が失敗し、波形数値の算出および心拍マップの作製ができないケースが多かった。このように従来の心拍解析の手法は、測定によって得られるすべての血流データについて適用できるというわけではなく、皮膚血流動態の指標化には限界があった。

以上、本発明が解決しようとする課題は、従来の波形数値とは異なる新たな動態指標として、拍動（ビート）の強度を表す心拍強度をマップ上に画像化し、個体間はもとよりマップ内の領域間で血流動態の違いを比較可能にすることである。これによって、これまで心拍検知に失敗していたデータを含むすべてのデータについて、血流動態の客観的数値指標となる心拍強度を得ることができる。

本発明の血流動態画像化診断装置は、血球を有する生体組織の観察領域にレーザー光を照射するレーザー光照射系と、生体組織の観察領域からの反射光を検出する複数の画素からなる受光部と、受光部からの信号に基づき連続的に複数の画像を取り込む画像取込部と、複数の画像を記憶する画像記憶部と、該記憶された複数画像の対応する各画素の出力信号の時間的変化から生体組織内の血流速度を演算する演算部と、該演算結果の二次元分布を血流マップとして表示する表示部とを備える。

演算部は、血流速度の演算により得られた血流マップの時間変化から、所定領域内における時系列血流データが有する基本周波数の信号強度を求め、該信号強度を基に心拍の強度を表す心拍強度を計算する心拍強度演算部を備える。表示部は、心拍強度演算部で計算した心拍強度を画像化して表示する。

また、表示部は、時系列に並んだ複数枚の血流マップを平均化することにより作成した合成マップに重畳して、心拍強度演算部で計算した心拍強度を表示する。また、表示部は、心拍強度と共に、その心拍強度に対応する周波数を画像化して表示する。この周波数は、心拍数に変換した数値として表示することができる。

また、表示部は、血流変化から生じる血流動態を数値化した波形数値を表す波形マップを表示したり、表層・背景分離情報に基づき作成した分離マップを表示することができる。また、表示部は、所定領域の表示の傍に、心拍強度或いは周波数の大きさを視覚的に表現することができる。

心拍強度演算部は、所定領域の時系列血流データを使用して、血流の時間変化をフーリエ級数展開して周波数空間へ変換し、該周波数空間に変換した血流データを基にパワースペクトルを予測し、この予測したパワースペクトルの分布内から最大パワーをもつ周波数を基本周波数として、それに対応する信号強度を求めて、心拍強度を計算する血流パワースペクトルモデル解析部と、二次元に分布する心拍強度を画像化する心拍強度マップ作成部を備える。また、血流パワースペクトルモデル解析部は、予測したパワースペクトルの分布内から予め設定した一つまたは二つ以上の周波数の範囲から、この周波数範囲に対応する信号強度を求めて、心拍強度を計算する。合成マップは、時系列に並んだ複数枚の血流マップを全て平均化することにより作成するか、若しくは、一つまたは二つ以上の拍動について決定した先頭マップを基準として、血流マップ同士を平均化した心拍マップを作成し、さらにこの心拍マップを平均化することにより作成することができる。

時系列血流データを取得する所定領域は、血流マップ上に観察者にとって興味ある領域を任意の形で設定した領域、又は、血流マップ全体、若しくは、血流マップを分割した領域である。
また、本発明の血流動態画像化診断方法は、レーザー光照射系によって血球を有する生体組織の観察領域にレーザー光を照射すると共に、複数の画素からなる受光部によって生体組織の観察領域からの反射光を検出し、画像取込部によって受光部からの信号に基づき連続的に複数の画像を取り込むと共に、画像記憶部によって複数の画像を記憶し、該記憶された複数画像の対応する各画素の出力信号の時間的変化から生体組織内の血流速度を演算し、該演算結果の二次元分布を血流マップとして表示するものであり、血流速度の演算により得られた血流マップの時間変化から、所定領域内における時系列血流データが有する基本周波数の信号強度を求め、該信号強度を基に心拍の強度を表す心拍強度を計算し、かつ、この計算した心拍強度を画像化して表示する。

従来、拍動による血流変化が少ないまったりとした血流変動では心拍検知をすることができないために、血流動態の客観的数値指標を得ることができないが、本発明によれば、心拍検知に失敗するような血流変動の場合であっても、すべてのデータについて血流動態の客観的数値指標となる心拍強度を得ることができる。

また、本発明によれば、心拍抽出処理なしに、すべての血流データについて血流動態指標を数値化して表示することにより、個体間はもとよりマップ内の各領域間における血流動態の違いを、定量的、視覚的に比較し、診断できる効果がある。

従来、心拍検知をすることができた場合には、歪度などの血流動態を数値化した波形数値を得ることができる。歪度は、動脈性の鋭い立ち上がり波形を有する部位と、静脈性の緩やかに上下する波形を有する部位を区別する数値である。しかし、このような波形数値は、末梢側では心拍検知に失敗することもあって、成否が安定しない。これに対して、心拍の強度を表す心拍強度は、確実にどこの場所でも血流動態を表現できる。

本発明に基づき構成した血流動態画像化診断装置の概略全体構成図の例を示す。 本発明が特徴とする演算部の構成図の例を示す。 血流パワースペクトルモデル解析部の演算例の流れ図を示す。 血流の時間変化を例示する図である。 作成した心拍強度マップの例である。 合成マップと心拍強度マップを重畳した表示例である。 ＬＡＳＳＯ法の最小化問題を解くために実現可能なコード例である。 図６の領域１から抽出した時系列の血流データをグラフ表示した例である。 図６の領域１から抽出した時系列の血流データから作成した予測パワースペクトルの例である。 図６の領域２から抽出した時系列の血流データをグラフ表示した例である。 図６の領域２から抽出した時系列の血流データから作成した予測パワースペクトルの例である。 図６の痛風患者が寛解した合成マップと心拍強度マップを重畳したマップの例である。 図１２の領域１から抽出した時系列の血流データをグラフ表示した例である。 図１２の領域１から抽出した時系列の血流データから作成した予測パワースペクトルの例である。 従来技術による血流動態画像化診断装置の概略全体構成図の例を示す。 従来技術による血流動態画像化診断装置の演算部の例を示す。 心拍検知に成功した場合の血流の時間変化を示す。 各拍動に対応する複数の血流マップから合成マップを作成する方法を説明する図である。 一心拍における血流の時間変化を示すグラフである。 心拍検知が行えた場合に作成した合成マップを例示している。 合成マップ上に分離マップ、および波形マップを重畳した例を示す。 心拍検知がうまく行えなかった場合を説明する図である。 心拍検知がうまく行えなかった場合に作成した合成マップを例示している。 心拍検知がうまく行えなかった場合に重畳表示されるマップを例示している。

以下、例示に基づき本発明を説明する。図１は、本発明に基づき構成した血流動態画像化診断装置の概略全体構成図である。レーザー光照射系２１は、皮膚血流や内臓・臓器血流、眼底血流等のような血流を有する生体組織２０（例えば、足背）にハーフミラー２２を介してレーザー光を照射する。受光部２３は、受光面上に多数の画素を有するＣＣＤ、ＣＭＯＳ等の固体撮像装置、レーザー反射光を固体撮像装置に結像する受光レンズ、及び撮像素子の出力を増幅する増幅回路などを有している。固体撮像装置は二次元に配置した撮像素子だけでなく一次元に配置した撮像素子を用いてもよい。タイミングパルス２５に基づき駆動される撮像素子は受光レンズにより結像された生体組織像を、タイミングパルス２５に基づいて電気信号に変換し、フレーム蓄積方式で信号電荷を読み出して映像信号として増幅して出力する。

画像取込部２４は、出力された映像信号に対して利得制御等のアナログ処理した出力をさらにＡ／Ｄ変換したデジタル信号、及びタイミングパルスに基づき、一心拍または二心拍以上の所定時間間隔（例えば１／３０秒間隔）で連続的に複数の画像を取り込む。画像記憶部２６は、この取り込まれた画像を記憶する。演算部２７は、記憶された複数画像の対応する各画素の出力信号の時間変化から生体組織内の血流速度を演算する。以上に説明した血流動態画像化診断装置の構成は、演算部２７を除いて図１を参照して前述した従来構成（特許文献８参照）と同一にすることができる。

本発明の演算部２７は、演算により得られた血流マップの時間変化から、血流マップ上の一つまたは二つ以上の関心領域内における心拍の強度を表す心拍強度を計算する。また、心拍強度と共に、その心拍強度に対応する周波数を出力することができる。表示部２８では、血流マップ情報に加え、心拍強度情報を重畳表示する。加えて、従来技術同様に、波形数値を表す波形マップ及び表層・背景分離情報に基づき作成した分離マップ等の診断に有用な情報を血流マップに重畳させ、二次元分布情報としてマップ表示できる。本発明の演算部２７の構成についての詳細は、図２を使って説明する。

図２は、本発明が特徴とする演算部２７を例示する構成図である。この演算部２７には、血流パワースペクトルモデル解析部３９及び心拍強度マップ作成部４０からなる心拍強度演算部３８が設けられており、この血流パワースペクトルモデル解析部３９には、関心領域情報と共に、血流解析部２９から血流マップ情報が入力される点、血流パワースペクトルモデル解析部３９の出力信号が心拍解析部３０に入力される点、及び心拍強度マップ作成部４０の出力信号が合成マップ・付加情報マップ重畳部３７に入力される点以外では、図１６を参照して前述した従来技術構成と同様に動作する。それ故に、重複する部分についての詳細な説明は省略する。

血流解析部２９は、画像記憶部２６から取り込まれた画像を解析して、生体組織での血流分布を二次元画像表示した血流マップを作成する。血流解析部２９は、複数画像の対応する各画素の出力信号の時間変化から、生体内の血流速度を演算する。この演算した血流値を、図４において、血流の時間変化として例示している。血流値は、通常脈動するが、図中上方の矢印４−２（４個）は、各々の拍動における最低血流マップの時間的な位置を示している。

心拍解析部３０では、上述の従来技術と同様に、各拍動ごとに最低血流マップを探索する手段から得られた複数の一心拍当たりのマップ数を求める。加えて、血流パワースペクトルモデル解析部３９で得られた心拍強度に対応する周波数をもとに、一心拍当たりのマップ数を求め、これら二つの手法で得た一心拍当たりのマップ数を比較し、最適な一心拍当たりのマップ数を決定する。この手法により、従来技術に比べ、より心拍検知が成功する確率が向上する。次に、図４のように、血流解析部２９で得た時系列血流マップの先頭から、決定した一心拍当たりのマップ数で順次複数マップを区切り（４−１）、各拍動における最低血流マップを探索する。図４の場合、４−３、４−４、４−５のように拍動ごとの血流マップの集団が区切られる。

これによって、心拍解析部３０は、得られた各心拍の最低血流マップを先頭マップとして決定し、心拍検出成功として心拍マップ作成部３１に進む。心拍解析部３０以降の処理は従来技術と同様である。簡単に言えば、心拍マップ作成部３１は、決定した先頭マップを基準として、二つ以上の拍動の場合は血流マップ同士を平均化した血流マップを作成する。一つの拍動でも心拍マップを作成することができるが、平均化処理は必要としない。時間軸に並べたこの平均化した複数の血流マップを、心拍マップと呼ぶ。合成マップ作成部３２では、心拍マップ作成部３１で作成された心拍マップを平均化して一枚の合成マップを作成する。

心拍マップ上の平均血流値を、時系列順に並べると、拍動による血流変化が確認でき、この血流変化を波形数値計算の基とする。波形数値作成部３３は、従来技術同様に、この拍動による血流変化から、関心領域内の歪度などの血流動態を数値化した波形数値を作成することができる。波形マップ作成部３４は、波形数値作成部３３で作成した波形数値を画像上に配置してマップ化した波形マップを作成することができる。

表層・背景分離部３５では、作成した合成マップを基に、基準とした閾値を境に、表層と背景を分離する。分離マップ作成部３６は、従来技術同様に、分離した表層と背景をそれぞれ二値化した分離マップを作成することにより表示を見易くすることができる。また、分離した表層、背景の分離情報は、それぞれ表層の血流領域の関心領域、背景血流の関心領域として各々再定義し、血流パワースペクトルモデル解析部３９の入力情報として使用しても良い。

心拍解析部３０で最低血流マップを決定することができなかった場合、即ち、心拍検知がうまく行えなかった場合、時系列に並んだ血流マップの全てが、合成マップ作成部３２に送られる。合成マップ作成部３２では、血流解析部２９で作成した時系列に並んだ複数枚の血流マップを、先頭のマップから最後のマップまでを平均化し、一枚の血流マップを作成する。これを合成マップとする。

表層・背景分離部３５は、心拍検知がうまく行えなかった場合、波形数値作成部３３で作成される波形数値を利用することができないが、合成マップ上の関心領域内の血流値を基にヒストグラムを得て、この分布から表層と背景を分離することが可能なので、合成マップのみでも分離マップを作成可能である。このように、心拍検知が失敗した場合には、合成マップ及び分離マップのみを作成して、合成マップ・付加情報マップ重畳部３７で重畳する。

なお、合成マップ自体は、上述のように、心拍マップが作成される場合は心拍マップを平均化することにより、或いは、心拍マップが作成されなくとも、時系列に並んだ複数枚の血流マップを全て平均化することにより作成することが可能である。但し、波形マップで表示すべき波形数値を作成するためには、心拍マップを作成することが必要となる。また、表層・背景分離部３５で、より精度の高い分離をするためには、心拍マップを作成することが望ましい。

次に、本発明の特徴とする心拍強度演算部３８について説明する。血流解析部２９は、画像記憶部２６に記憶された複数画像の対応する各画素の出力信号の時間変化を取得して解析し、その解析結果として、時系列に並んだ複数枚の血流マップは、心拍解析部３０だけでなく、血流パワースペクトルモデル解析部３９にも出力される。

血流パワースペクトルモデル解析部３９は、関心領域毎に時系列に並んだ血流マップ上の血流値を平均化して、時系列血流データを作成する。関心領域は、血流マップ上に一つあるいは二つ以上の観察者にとって興味ある領域を矩形、楕円、ポリゴン、グリッドなど任意の形で設定した領域である。血流マップ全体の場合や、血流解析部２９で解析した血流マップを画像処理して得られるマップ上に適切に任意の形状で分割した領域である場合もありうる。

血流パワースペクトルモデル解析部３９では、関心領域毎に時系列血流データを入力として心拍の強度を表す心拍強度と心拍強度に対応する周波数を出力することができる。時系列血流データとしては、関心領域内に存在する血流値からある程度の数、例えば１０００程のサンプルを抽出し、抽出した血流値を基に平均化などの算術演算をしてもよい。また時系列血流データは、各心拍の血流値が最大値になる経過時間を採用する方法を用いてもよい。

心拍強度は、時系列血流データを基に周波数強度の分布であるパワースペクトルを求めた後に、得られたパワースペクトルの特徴、つまり起伏のあるパワーの強弱の形状を基にしている。例えばパワースペクトルの分布内から最大パワーをもつ周波数に対応するパワーに定数を乗した数値を心拍強度とし、この心拍強度に対応する周波数と合わせて比較できる血流動態指標とする。

時系列血流データをフーリエ展開すれば、第一基本波の振幅が心拍強度に対応して、第一基本波の周波数（基本周波数）が、心拍強度をもつ周波数に対応する。本発明では心拍強度の計算に必要なフーリエ展開したときに得られるフーリエ級数の係数を予測している。ここで予測とは、非常に長い時系列血流データを取得しなければ得られないパワースペクトルを、限られた短い時系列血流データを用いて前もって推量することである。

心拍強度マップ作成部４０では、関心領域別に拍動の強度をマップ化した心拍強度マップを作成する。図５は、心拍強度マップを単独で例示する図である。複数（２個の円として例示）の関心領域１，２の傍に、バーグラフで心拍強度の大きさを視覚的に表現している例である。

心拍強度に対応する周波数は、心拍数に変換した数値を表示することもできる。変換した数値が心拍数に近い数値である場合、拍動によるビートであるか、また末梢血流によく観察される交感神経、副交感神経系の血流変動なのか観測者が判断できるので、関心領域の血流動態の把握に有用であり、臨床上診断に有用である。

合成マップ・付加情報マップ重畳部３７では、合成マップ作成部３２で作成した合成マップを下層にし、心拍強度マップ作成部４０で作成した関心領域別に拍動の強度をマップ化した心拍強度マップ（図５）を上層にし、これらマップを図６のように重畳する。

図６に示すように、関心領域１，２の傍にバーグラフで心拍強度の大きさを視覚的に表現している。また、図中に６−１、６−２で示すように、上の行に心拍強度を、下の行に周波数を数値で表示している。バーグラフでは、表示する数値の程度を色分けして、高い数値ほど赤などの暖色にし、低い値では青などの寒色にしてもよい。心拍強度または心拍強度に対応する周波数は、他の有効なグラフや数値表現などを使って描画してもよい。前回測定したときの心拍強度などの数値を装置に記憶しておき、前回値と比較するグラフも描画できる。

これによって、合成マップ・付加情報マップ重畳部３７は、合成マップの関心領域に一致して、関心領域の心拍強度による拍動の強度分布を配置したマップを一枚作成し、図１の表示部２８で表示する。表示部２８では、合成マップ、心拍強度マップを、重畳することなく、それぞれ別のマップとして表示してもよい。また追加の画像処理及び追加情報をマップ上に描画し、これらのマップを強調し表示してもよい。

さらに、合成マップ・付加情報マップ重畳部３７では、心拍検知に成功した場合には、血流動態を波形数値として数値化した波形マップとか、領域情報を付加した分離マップ（従来技術の図２１参照）など複数のマップを合成マップの上層に重畳しても良い。

次に、心拍強度を計算する血流パワースペクトルモデル解析部３９内部の演算の流れについて、図３の流れ図を使ってさらに詳細に説明する。

まずS1で、関心領域に対する時系列血流データを取得する。取得方法については先に述べたとおりである。時系列血流データは途中マップが欠損したデータでもよい。また時間軸に等間隔に並んでいる必要もない。ｎ個の時系列血流データY（ｙ１、ｙ２、・・・、ｙｎ）にそれぞれ先頭マップからの時間を計算した時間情報T（ｔ１、ｔ２、・・・、ｔｎ）を付加する。

一般的に時系列データをフーリエ変換すると、パワースペクトルが得られる。心拍強度はパワースペクトルを基に算出しており、パワースペクトルを作成することが必要になる。S2では、S1で作成した時間情報を基に、次の周波数空間にフーリエ級数展開する行列を数１のように計算する。

ここでMはパワースペクトルの周波数分割数であり、例えば１００を設定できる。またｆ_{1 、}ｆ_2、・・・、ｆ_Mは、パワースペクトルの周波数要素であり、例えば、０．１、０．２、_、・・・、１０[Hz]を設定できる。

時系列血流データは、S２で作成したフーリエ級数展開する行列を使って以下のように表される。

ここでXはパワースペクトルの周波数要素ごとの係数のベクトルであり、次のように２M個の要素からなる。

S3では、周波数空間に変換した時間情報と時系列血流データから、パワースペクトルを作成する。関心領域の時系列血流データのパワースペクトルは、ある特定の周波数、例えば足の表面の血流測定の場合、拍動がよく表れる領域では心拍数に同期した周波数や、交感神経、副交感神経支配下の０．５[Hz]以下の周波数のパワーが強く表れる。これらの生体にかかわる周波数の信号はスパースであると考えられる。スパースとはこの場合、パワースペクトル上の大部分の周波数成分がほぼ０であり、拍動信号や交感神経、副交感神経系の信号などの意味のある信号成分ぐらいしかない状態である。パワースペクトルがスパースであるということは、限られた時系列血流データから、パワースペクトルを予測できる。

スパースである場合、作成には例えば、LASSO（ｌｅａｓｔ ａｂｓｏｌｕｔｅ ｓｈｒｉｎｋａｇｅ ａｎｄ ｓｅｌｅｃｔｉｏｎ ｏｐｅｒａｔｏｒ）法、またはＲｅｄｇｅ回帰等、別の回帰アルゴリズムを使用して、最適な予測パワースペクトルを作成できる。

パワースペクトルの計算は、関心領域内の平均化した血流値ではなく、測定点毎に時系列データを取得してフーリエ変換を行うのが通常である。このため、長時間の測定した時系列データを必要としていた。例えば心電図のパワースペクトルを得るための測定では、１００から６００拍ほどの心電のデータを取得し、得られたデータからフーリエ変換を行って、周波数解析が行われている。予測パワースペクトルを作成することにより、これまでのような長い時間をかけた測定を必要とせず、４秒程度の短い時系列血流データを基に、０．１[Hz]程度の低周波数から例えば１００[Hz]までのパワースペクトルを予測して、周波数分布を得ることができる。パワースペクトルを予測し周波数分布を得るために、例えば、LASSOなどの手法（実施例１参照）を使うことができる。

求めた最適な予測パワースペクトルの係数列を以下のように定義すると、

iの総数は周波数分割数であり、例えばi=１,2,・・・１００である。

次にＳ４では、Ｓ３で求めた予測パワースペクトルを解析し、心拍強度を一つ計算する。心拍強度は予測パワースペクトルを基に計算した数値であり、例えば以下の式で計算する。

ここでＢＳを満たすf_iを心拍強度に対応する周波数とする。Cは心拍強度のスケールを調整する比例定数である。

心拍強度は、予め設定した一つまたは二つ以上の周波数の範囲内から求めることができる。例えば０．０５〜０．１５[Hz]や、０．１５〜０．４０[Hz]などである。心拍強度の他の演算としては、予測パワースペクトルの平均値を求める方法、モードをもつパワースペクトルの分布の山になる局所的な面積を求める方法、予め設定した周波数範囲を積分する方法などがあり、これらを採用してもよい。

Ｓ５では、Ｓ４で解析した結果である心拍強度と、心拍強度に対応する周波数を関心領域毎に記憶部に保存できる。一度計算した心拍強度と心拍強度に対応する周波数は、記憶域に保存しているため関心領域が変更される都度、計算しなくてもよいというメリットがある。必要に応じて、血流パワースペクトルモデル解析部３９から、関心領域毎の心拍強度情報を、心拍強度マップ作成部４０に引き渡すことができる。

Ｓ５で記憶した心拍強度に対応する周波数は、血流波形が拍動に同期した形状である場合、６０を乗ずると心拍数に一致する。例えば関心領域が皮膚表面の動脈性の拍動を含む領域であれば、一般成人の場合、心拍数は５０〜９０であり、心拍強度に対応する周波数は０．８〜１．５[Hz]になる。

関心領域が皮下の小動脈などをあまり含まず、末梢血流が支配的な領域では拍動が検出されず、心拍強度に対応する周波数は０．５[Hz]未満になることが多い。この低周波の領域のパワーは、交感神経系、副交感神経系の活動の程度を表しているという報告もあり、この心拍強度を使って、人、動物など測定対象に対して、末梢領域での神経系の活動状況を判定できる可能性が高い。

本発明者らの検討から、２秒以上のデータ、つまり毎秒３０フレームの固定撮像装置を使用する場合、６０フレーム程度からパワースペクトルを予測可能である。本発明者らの開発した装置を使用し、心拍強度を計算するために最適なフレーム数は、４秒ほどの時系列血流データであり、１２０フレームのデータから末梢血流の血流動態を把握できる心拍強度が得られる。時系列血流変動が６秒程度になってくると、心拍に同期した血流変化よりも交感神経、副交感神経系の信号強度が強くなり、末梢血流の動態様態を識別する検出力が弱まる傾向がある。

血流パワースペクトルモデル解析部３９で作成した関心領域毎の心拍強度は、関心領域の数だけ存在し、二次元のマップ上に対応する。心拍強度マップ作成部４０では、空の仮想マップを作成し、合成マップに対応した心拍強度の数値を、該当する仮想マップ上の位置に描画する。数値の視認性を上げるため、関心領域の近くまたは関心領域内にバーグラフなど数値を視認性良く確認できるグラフを描画してもよい。

一例としてLASSO法を用いた予測パワースペクトルを作成する方法について説明する。時系列血流データおよび最適な予測パワースペクトルの係数列を以下のようなモデル関数に表すことができる。

という意味である。

である。

数６及び数７を満たす最小化問題を解くことによって、最適な予測パワースペクトルの係数列を解くことができる。ＬＡＳＳＯ法の最小化問題を解くために実現可能なコード例（例えばＲ言語で記述したもの）を図７に記載する。

本発明による効果的な実施例を提示するため、実際に本書で説明した血流動態画像化診断装置を構成し、実際に足に疾患のある一例として痛風患者の足の血流を測定した結果を図６に示す。左足が炎症部で患者は痛みを伴っており、色が白いところほど血流値が高く、黒色のところほど血流が低い。関心領域を円で２か所設置しており、関心領域１は炎症がひどい部位で、関心領域２は炎症がさほど見られない部位である。

図８は関心領域１の血流値の経時変化をプロットしたもので、4秒間の血流値の変化を表している。ひどく炎症を起こしている結果、患部は腫れ上がり血流も炎症作用の影響から血流値が高く、また血流の変動も拍動に同期して大きく変化していることが確認できる。血流波形は収縮期に急峻に立ち上がり、いったん下がった後、バウンスするように再度跳ね上がる二峰性の動態変化を有している。

図８の信号を基に本発明の機能を用いてパワースペクトルを描画した結果が図９である。パワースペクトル内でパワーの最大値を持つ周波数（心拍強度に対応する周波数）は１．２[Ｈｚ]であることがグラフからわかり、血流パワースペクトルモデル解析部の出力と一致していた。心拍数を計算すると７２回であるので、拍動に同期している信号であることが推測される。合わせて心拍強度を算出した結果は、BS２０．３であった。BS（任意単位）は、パワースペクトルをもとに計算したパワーの最大値を意味している。グラフ縦軸はパワースペクトルのパワーを表している。

図８の波形は二峰性の動態変化であり、二峰性の信号では周波数が高いところにもピークがあったが、パワースペクトルを観察すると拍動に同期した信号のパワーのほうが、よりパワーが強い。本例は、得られた心拍強度が拍動による動態指標として、より的確に計算できた一例である。

一方で、関心領域２について、血流値の経時変化のプロットを図１０に示す。経時変化の様子を見る限り、関心領域１の血流値の経時変化である図８に比べ、拍動の立ち上がりがなだらかであり、またピーク後の立ち下がりも図８ほど急峻でなく、なだらかに減少していることが確認できる。

関心領域２のパワースペクトルを図１１に示す。心拍強度に対応する周波数は、１．２[Ｈｚ]であり、心拍強度は５．６と計算された。拍動に同期したパワーは、関心領域１に比べて小さかった。体のどこかの部位が腫れている経験があれば、拍動に同期して大きく脈（＝ビート）を打つことを感覚的に経験しているが、この症例では、腫れている部位（関心領域１）と、そうでない部位（関心領域２）の脈を打つ程度が、同一の足でありながら、心拍強度により領域間の差異を客観的な指標として比較できている一例である。

図６の痛風患者が、後日、寛解した時の合成マップの様子及び心拍強度を測定した結果を図１２に示す。痛風時に設定した関心領域１を、図６に同様に関心領域１として設定した。血流の経時変化を図１３に、パワースペクトルを図１４に示す。

関心領域１の心拍強度は、１．４であり、図６の時に見られている炎症が治まっているおかげで、拍動はなだらかで静かな血流動態になっていることが分かる。図１３の経時変化からも拍動に同期した大きな血流変動は見られない。

本来、末梢側の血流は、末梢血管網に対して適度な流れを維持しており、はっきりとした血流変動は見出しづらい。これは小川から大きなため池に流れ込む水の流れを想定すると、理解しやすい。流れこむ大きなため池が大きければ大きいほど、ため池の表面の波は穏やかで静かな変動になる。一方で、小川が激流となり池に流れ込んでくると、ため池の表面は大きく波面を乱すことになり、表面の変動は大きくなり変動のパワーは強くなる。

図１２の測定時には炎症が治まり、血流変動は落ち着いているので、足に流れ込む血流は阻害されておらず、適度な流れに調節されていることが心拍強度の数値からも分かる。さらに、十分なサンプル数が得られれば、心拍強度による末梢血流動態の把握および鑑別も十分可能である。

臨床分野では、形成外科の分野で皮弁移植や植皮など皮膚の移植を行っているが、皮膚の定着の程度を判定する際に、本発明に基づき構成した血流動態画像化装置の果たす役割は大きい。皮膚が定着することは、いわば血流の再建にほかならず、血流動態の把握は今後重要な検討項目になってくるに違いない。非侵襲的に測定できることから、被験者に苦痛を与えず、またレーザーを照明するといっても対象に大きなエネルギーを与えないため、被験者に重篤な被害をもたらさないことも大きな利点である。

これまで本発明者らは、皮膚血流をはじめ内臓・臓器血流、或いは眼底血流の血流動態を測定できる血流画像化装置を開発してきた。これらの開発品は、時間分解能において同等のスペックを有しており、同等の時間分解能があれば、皮膚血流と同様に血流が観察できる対象であれば、動態指標の数値化が可能である。本発明を適用する範囲は、皮膚血流のみならず、これらのこれまで測定が行えた対象について、心拍強度を得られることは容易に想像がつく。

以上、本開示にて幾つかの実施の形態のみを単に一例として詳細に説明したが、本発明の新規な教示及び有利な効果から実質的に逸脱せずに、その実施の形態には多くの改変例が可能である。

２０ 生体組織
２１ レーザー光照射系
２２ ハーフミラー
２３ 受光部
２４ 画像取込部
２５ タイミングパルス
２６ 画像記憶部
２７ 演算部
２８ 表示部
２９ 血流解析部
３０ 心拍解析部
３１ 心拍マップ作成部
３２ 合成マップ作成部
３３ 波形数値作成部
３４ 波形マップ作成部
３５ 表層・背景分離部
３６ 分離マップ作成部
３７ 合成マップ・付加情報マップ重畳部
３８ 心拍強度演算部
３９ 血流パワースペクトルモデル解析部
４０ 心拍強度マップ作成部

Claims (12)

1. 血球を有する生体組織の観察領域にレーザー光を照射するレーザー光照射系と、前記生体組織の観察領域からの反射光を検出する複数の画素からなる受光部と、前記受光部からの信号に基づき連続的に複数の画像を取り込む画像取込部と、前記複数の画像を記憶する画像記憶部と、該記憶された複数画像の対応する各画素の出力信号の時間的変化から生体組織内の血流速度を演算する演算部と、該演算結果の二次元分布を血流マップとして表示する表示部からなる血流動態画像化診断装置において、
   前記演算部は、血流速度の演算により得られた血流マップの時間変化から、所定領域内における時系列血流データが有する基本周波数の信号強度を求め、該信号強度を基に心拍の強度を表す心拍強度を計算する心拍強度演算部を備え、
   前記表示部は、前記心拍強度演算部で計算した心拍強度を画像化して表示することから成る血流動態画像化診断装置。
2. 前記表示部は、時系列に並んだ複数枚の血流マップを平均化することにより作成した合成マップに重畳して、前記心拍強度演算部で計算した心拍強度を表示する請求項１に記載の血流動態画像化診断装置。
3. 前記表示部は、前記心拍強度と共に、その心拍強度に対応する周波数を画像化して表示する請求項２に記載の血流動態画像化診断装置。
4. 前記周波数は、心拍数に変換した数値として表示する請求項３に記載の血流動態画像化診断装置。
5. 前記表示部は、血流変化から生じる血流動態を数値化した波形数値を表す波形マップを表示する請求項２に記載の血流動態画像化診断装置。
6. 前記表示部は、表層・背景分離情報に基づき作成した分離マップを表示する請求項２に記載の血流動態画像化診断装置。
7. 前記表示部は、所定領域の表示の傍に、心拍強度或いは周波数の大きさを視覚的に表現する請求項２に記載の血流動態画像化診断装置。
8. 前記心拍強度演算部は、前記所定領域の時系列血流データを使用して、血流の時間変化をフーリエ級数展開して周波数空間へ変換し、該周波数空間に変換した血流データを基にパワースペクトルを予測し、この予測したパワースペクトルの分布内から最大パワーをもつ周波数を前記基本周波数として、それに対応する信号強度を求めて、前記心拍強度を計算する血流パワースペクトルモデル解析部と、二次元に分布する心拍強度を画像化する心拍強度マップ作成部を備えた請求項１に記載の血流動態画像化診断装置。
9. 前記心拍強度演算部は、前記所定領域の時系列血流データを使用して、血流の時間変化をフーリエ級数展開して周波数空間へ変換し、該周波数空間に変換した血流データを基にパワースペクトルを予測し、この予測したパワースペクトルの分布内から予め設定した一つまたは二つ以上の周波数の範囲から、この周波数範囲に対応する信号強度を求めて、前記心拍強度を計算する血流パワースペクトルモデル解析部と、二次元に分布する心拍強度を画像化する心拍強度マップ作成部を備えた請求項１に記載の血流動態画像化診断装置。
10. 前記合成マップは、時系列に並んだ複数枚の血流マップを全て平均化することにより作成するか、若しくは、一つまたは二つ以上の拍動について決定した先頭マップを基準として、血流マップ同士を平均化した心拍マップを作成し、さらにこの心拍マップを平均化することにより作成する請求項２に記載の血流動態画像化診断装置。
11. 時系列血流データを取得する前記所定領域は、血流マップ上に観察者にとって興味ある領域を任意の形で設定した領域、又は、血流マップ全体、若しくは、血流マップを分割した領域である請求項１に記載の血流動態画像化診断装置。
12. レーザー光照射系によって血球を有する生体組織の観察領域にレーザー光を照射すると共に、複数の画素からなる受光部によって前記生体組織の観察領域からの反射光を検出し、画像取込部によって前記受光部からの信号に基づき連続的に複数の画像を取り込むと共に、画像記憶部によって前記複数の画像を記憶し、該記憶された複数画像の対応する各画素の出力信号の時間的変化から生体組織内の血流速度を演算し、該演算結果の二次元分布を血流マップとして表示する血流動態画像化診断方法において、
    血流速度の演算により得られた血流マップの時間変化から、所定領域内における時系列血流データが有する基本周波数の信号強度を求め、該信号強度を基に心拍の強度を表す心拍強度を計算し、かつ、この計算した心拍強度を画像化して表示することから成る血流動態画像化診断方法。
    

Images