JP6784403B2

JP6784403B2 - 心拍変動の推定方法、心拍変動の推定プログラム及び心拍変動推定システム

Info

Publication number: JP6784403B2
Application number: JP2017168657A
Authority: JP
Inventors: 津村　徳道; 徳道津村; 亮太三橋; 弦樹岡田
Original assignee: Chiba University NUC
Current assignee: Chiba University NUC
Priority date: 2017-09-01
Filing date: 2017-09-01
Publication date: 2020-11-11
Anticipated expiration: 2037-09-01
Also published as: JP2019042145A

Description

本発明は、近赤外光を用いた心拍変動の推定方法、心拍変動の推定プログラム及び心拍変動推定システムに関する。

非特許文献１には、以下の技術が記載されている。防犯カメラ等で広く用いられる１バンド近赤外カメラで撮影した動画像に対して、設定した顔画像の関心領域の平均画素値を時間軸に解析することで脈波信号の原信号を取得する。その後、経験的モード分解による信号の傾きを除去した後、周波数解析により心拍数を推定する技術が記載されている。

W.Zeng, Q. Zhang, Y. Zhou, G. Xu, and G. Liang, "Infrared Video based Non-invasive Heart Rate Measurement,"in Proceedings of the 2015 IEEE Conference on Robotics and Biomimetics. IEEE, 2015, pp. 1041‐1046.

上記の従来手法では照明環境が一定の場所を想定して測定が行われていた。しかし、実環境においてこの技術を実装することを考えた際、多くの場合で照明環境は変動することが予想され、現状では照明環境が変動すると心拍数の推定精度が大きく低下する。従って、照明環境が変動するような実環境においても安定した心拍数の推定、という課題が挙げられる。

本発明の一つの観点によれば、上記課題を解決するために、心拍変動推定方法を、対象者の少なくとも２種類の波長域の動画像を撮影し、撮影した２種類の波長域の動画像に基づいて照明変動成分及びヘモグロビン成分を算出し、算出したヘモグロビン成分に基づいて心拍変動の推定を行うものとした。

また、本発明の他の観点によれば、心拍変動推定方法を、対象者の少なくとも２種類の波長域の動画像を撮影し、撮影した動画フレーム中の各画素の対数をとることで変換される対数空間上において、照明変動成分と、ヘモグロビン成分の少なくとも２成分に基底変換を行い、基底変換した情報に基づいて心拍変動の推定を行うものとした。さらに、少なくとも２種類の波長域のうち、１種類の波長域が770nm以上790nm以下であり、他の１種類の波長域が890nm以上910nm以下とすると望ましい。

また、本発明の他の観点によれば、心拍変動推定プログラムを、コンピュータに、対象者の少なくとも２種類の波長域の動画像を撮影する手順と、前記撮影した２種類の波長域の動画像に基づいて照明変動成分及びヘモグロビン成分を算出する手順と、前記算出したヘモグロビン成分に基づいて心拍変動の推定を行う手順とを実行させるためのプログラムとした。

また、本発明の他の観点によれば、心拍変動推定システムを、対象者の少なくとも２種類の波長域の動画像を撮影するカメラと、前記カメラにより撮影した２種類の波長域の動画像に基づいて照明変動成分及びヘモグロビン成分を算出し、前記算出したヘモグロビン成分に基づいて心拍変動の推定を行う心拍変動推定部と、を備えるものとした。

本発明によれば、照明環境が変動するような環境においても、安定して心拍変動を推定することができる心拍変動の推定方法、心拍変動の推定プログラム及び心拍変動推定システムを提供することができる。

本発明の実施形態の概要を示す図である。対数空間上における画素値の分布を示す図である。１バンドで推定した心拍数と本手法である２バンドで心拍数を推定した結果を示す図である。肌の構造モデルを示す図である。仮想肌の分光反射率曲線を示す図である。シミュレーションを用いた分光曲線データを示す図である。ヘモグロビン成分の推定手法を示す図である。ノイズの量と評価値の関係を示す図である。撮影環境を示す図である。近赤外２バンドカメラによる撮影画像を示す図である。濃度空間における２バンド近赤外画像の画素地の分布を示す図である。中心波長780nmの１バンド近赤外から得られた脈波信号を示す図である。傾き除去後の脈波信号（中心波長780nm）を示す図である。傾き除去後の中心波長780nmの脈波信号のパワースペクトル密度を示す図である。バンドパスフィルタ後の中心波長780nmの脈波信号を示す図である。脈波信号のピーク検出を示す図である。近赤外１バンドカメラによる撮影画像を示す図である。照明変動がない時の近赤外１バンドカメラによる撮影画像を示す図である。照明変動がない時のバンドパスフィルタ後の脈波信号を示す図である。照明変動がある時の近赤外１バンドカメラによる撮影画像を示す図である。照明変動がある時のバンドパスフィルタ後の脈波信号を示す図である。ヘモグロビン・陰影成分分離を示す図である。照明変動がない時の信号処理後の脈波信号を示す図である。照明変動がある時の信号処理後の脈波信号を示す図である。推定した心拍数と相対誤差率を示す図である。１バンド測定、２バンド測定における推定心拍数との相対誤差率を示す図である。濃度空間における画素値の分布を示す図である。

本発明の実施形態の概要を、図１を参照しながら以下説明する。防犯カメラ等で広く用いられている１バンド近赤外カメラを２バンド近赤外カメラ（780nm±10nmと900nm±10nm）に拡張することで照明変動成分の影響を低減する。撮影した動画フレーム中の各画素の対数を取ることで変換される対数空間上において、照明変動成分と、血液の主成分であるヘモグロビン成分の２成分に基底変換を行う。この変換により、照明変動成分に影響を受けることなく心拍数の推定を行う。

成分分離を行ったうえで所望の成分のみに対してノイズ除去等の信号処理を行うことで被写体の動きや照明変動などに起因した陰影変化の影響を除外しつつ、生体情報を抽出する。

実際に、対数空間における撮影動画像の画素値の分布を図２に示す。図から分かるように照明変動成分は(1, 1)ベクトル方向に、もう一方向にヘモグロビン成分の変動が観測されていることが確認できる。

照明環境が変動している測定環境において、従来手法の１バンドで推定した心拍数と本手法である２バンドで心拍数を推定した結果を図３に示す。これより、従来手法では脈波の原信号に照明変動が大きく影響し心拍数の推定に大きな影響を与えているのに対して、提案手法では照明変動に影響を受けることなく、ヘモグロビン成分の変動をおよそ96.27%の精度で推定されていることが確認できる。

本発明の実施例を以下詳細に説明する。

第1章序論
1.1 背景
近年、夜間における精神状態のモニタリング技術が必要とされている。例えば、接触型の心拍計などの取り付けが困難な患者や乳幼児を、夜間でも連続的にモニタリングすることができれば、いち早く異常を検知することができる。また、夜間に運転中、もしくはトンネルの中などの照明環境が変化するような環境におけるドライバーを連続的にモニタリングすることができれば、ドライバーの眠気を検知するドライブモニタリングといった技術にも応用することができる。

また、近年カメラを使った非接触による生体情報の検出に関する研究が報告されている。しかし、これらの従来手法は照明環境の変化がない場合の計測であり、照明環境の変動が生じることにより心拍数の推定密度が低下するという問題点が存在する。

そこで、本検討では先行研究の1バンド近赤外画像に対して、2バンド近赤外画像を用いることで照明環境の変動に依存しない、暗所でも利用可能な生体情報の検出手法の提案を目的とした。

1.2 本実施例の構成
第2章では、近赤外領域での肌光学モデルについて述べる。第3章では、肌光学モデルを基に、2つの近赤外波長を用いて脈波測定を行う手法について述べる。第4章では、脈波測定を行うシミュレータの作成と、その評価方法について述べる。第5章では、実際に被験者を対象とした脈波測定について述べる。最後に第6章で、本実施例のまとめと今度の課題を述べる。

第2章近赤外領域での肌光学モデル
2.1 近赤外光
私たちが生活するために欠かせない光は、紫外線、可視光線、赤外線の三つに大別される。紫外線はおよそ100[nm]~380[nm]の波長域を持ち、高エネルギーを持つことから皮膚や目、免疫系へダメージを与える可能性がある。可視光線はおよそ380[nm]~780[nm]の短波長の波長域を持ち、普段私たちが目で見ている物体は、物体が反射した可視光が私たちの目に入ることでその物体を知覚している。赤外線は長波長の波長域を持ち、波長域により近赤外、中赤外、遠赤外の三つに大別され、近赤外線は目で見ることができないため、防犯カメラ等の犯罪防止を目的とした用途として一般的に普及している。そこで、本実施例では、およそ780nm~1,000nmの波長域を持つ近赤外線を用いて脈波測定を行う。

2.2 皮膚のモデルと観測信号
人間の皮膚は多層構造をとっており、表皮、真皮、そして皮下組織の三つに大別される。皮膚にはメラニン、ヘモグロビン、ビリルビンなど様々な色素が含まれ、皮膚の色調はこれらの色素により左右される。これらの色素の中でもメラニン、ヘモグロビンによる変化が皮膚の色調に大きく影響している。メラニン色素は表皮に、ヘモグロビン色素は真皮に毛細血管が張り巡らされているため真皮に多く存在している。従って、表皮をメラニン層、真皮をヘモグロビン層と仮定し、メラニン色素とヘモグロビン色素を空間的に独立とみなせる。このことから、可視光源の照明環境下においては、人間の肌を図４ (a)に示すような、表皮のメラニン層、真皮のヘモグロビン層からなる肌の二層構造として扱うことができる。先行研究ではこの仮定の下で解析を行い、濃度空間上で、陰影成分を除去した観測信号に独立成分分析を行うことで、ヘモグロビン成分とメラニン成分を推定した。

先行研究は通常のRGBカメラで行った研究であり、RGBカメラから得られる画像は、可視光から得られた情報のみで構成される。近赤外光は可視光よりも長波長であることから生体内部への透過深度が深く、毛細血管よりさらに奥の細動脈と呼ばれる、生体深部の血管を測定できると言われている。従って、近赤外領域の光が肌に入射した際の振る舞いは、図４(b)に示すような皮膚の深部のみでの反射が生じると仮定できる。すなわち、表皮における反射は考えず、細動脈付近でのヘモグロビンによる反射が生じると考えられる。本実施例では図４に示したような肌の一層構造モデルを想定し、この仮定をもとに検討を行った。

吸光度分布より、仮想肌における肌の反射率Rは以下の数式１（式(2.1)）で表される。

ここで、aは酸化ヘモグロビンの色素濃度、は波長の際の酸化ヘモグロビンの吸光度を表す。また、近赤外光を肌の一層構造モデルに照射した際のモノクロカメラのセンサ応答Iは数式２（式（2.2））で表される。

ここで、Rは肌の反射率、Cはカメラ感度、Bは近赤外透過フィルタに使用したバンドパスフィルタの感度、Sは照明強度、λは波長を表す。照明強度Sは本来、照明のちらつきにより変動するが、本実施例では、照明強度は一定とした。また、皮また、肌への入射光を常に単色光と仮定し、周辺の光源情報を除去する必要があるため、皮膚画像の撮影は暗室化で行われる。ここで、式(2.1)の対数を取ることにより観測信号を画像空間から濃度空間に変換する。本実施例では、この空間において新たな基底ベクトルを決定することによりヘモグロビンベクトルの推定を行った。

第3章シミュレータによる近赤外透過フィルタの最適化
3.1 はじめに
本章では、第2章で述べた近赤外領域における肌光学モデルに基づき、脈波計測に最適な2バンド近赤外波長の組み合わせを決定するためのシミュレータを作成した。本実施例では、酸化ヘモグロビンの吸光度曲線を用いることで実際の肌を想定した仮想肌を作成し、使用するモノクロカメラ、近赤外フィルタの分光曲線を用いることで、2バンド近赤外画像を撮影するカメラシステムを構築した。

3.2節では仮想肌におけるヘモグロビンの分布について述べる。3.3節では、ヘモグロビン分布の脈動による変動について述べる。3.4節では、仮想肌を撮影した際の画像化方法について述べる。 3.5節では、算出した画素値からヘモグロビン成分を推定する手法を述べる。最後に3.6節では、推定したヘモグロビン成分の耐ノイズ性を評価する関数について述べる。

3.2 ヘモグロビンの分布
人間の血液は絶えず循環しているため、肌の各部では、それぞれ異なるヘモグロビンの色素が分布している。そのため、本シミュレータでは10×10の小領域画像を作成し、その各画素位置にそれぞれ異なるヘモグロビン成分を分布させた。本実施例では、ヘモグロビンの分布を実際の肌の反射率と近づけるために、ヘモグロビンの色素濃度を0.20~0.213の範囲で100段階変化させて小領域画像に分布させることで仮想肌でのヘモグロビンの分布の様子を想定した。図５に仮想肌の分光反射率を示す。

3.3 脈動によるヘモグロビンの変動
3.2節では、仮想肌におけるヘモグロビンの分布について述べた。しかし実際の肌におけるヘモグロビン成分は、心臓から全身に送られる血液の脈動によりその値が変動し、時間変化によって異なる分布を示すと考えられる。この脈動による変化が生じる前のヘモグロビン成分の値と、脈動によって変化が生じた後のヘモグロビン成分の値を本シミュレータにより推定することで、ヘモグロビン成分を検出することができる。今回は、実際の肌における脈動を再現するために、0から0.00068の範囲でヘモグロビンの色素濃度を変動させた。

3.2節と本節より、ヘモグロビンの分布と、時間変化によるその変動を想定した仮想肌モデルを作成した。本実施例では、この脈動による酸化ヘモグロビンの時間変化を、脈波信号として扱う。以降のシミュレーションでは、この仮想肌モデルを用いてシミュレーションを行う。

3.4 画像化手法
3.2節、3.3節で作成した仮想肌モデルをカメラで撮影したときの画素値vは以下の数式３（式(3.1)）で表される。

ここで、Rは肌の反射率、Cはカメラ感度、Bは近赤外透過フィルタに使用したバンドパスフィルタの感度、Sは照明強度、は波長を表す。照明強度Sは本来、照明のちらつきにより変動するが、本実施例では、照明強度は一定としてシミュレーションを行った。カメラ感度Cは、図６ (a)に示すモノクロカメラ（DMK0234UV、アルゴ社）の分光曲線を使用した。バンドパスフィルタの感度Bは、特定の周波数だけ透過する矩形波を用いた。中心波長はそれぞれ780[nm]、820[nm]、860[nm]、900[nm]、940[nm]の5種類であり、半値全幅は10[nm]である。それに加えて、5つの干渉型バンドパスフィルタ（Edmund Optics Japan社）の分光曲線を同様に、近赤外光を透過するパンドパスフィルタの感度Bとして使用した。中心波長は矩形波と同じくそれぞれ780[nm]、820[nm]、860[nm]、900[nm]、940[nm]の5種類であり、半値全幅は10[nm]である。5つのバンドパスフィルタの分光分布曲線を重ね合わせた波形を図６(b)に示す。

シミュレーションは、矩形波の中から2種類選択し、各フィルタ感度のときの画素値を算出した。干渉型バンドパスフィルタの際も同様に2種類のバンドパスフィルタを選択し各フィルタ感度の際の画素値を、一般的な8ビットカメラを想定し、量子化ビット数を8ビットとして算出した。これにより、任意の2バンド近赤外画像が取得される。バンドパスフィルタの組み合わせは、波長が離れていくにつれ、酸化ヘモグロビンが高い吸光度特性を持つことから、780[nm]のバンドパスフィルタを基準に他の4つのバンドパスフィルタとの組み合わせを検証した。

3.5 ヘモグロビン成分の推定
3.4節で得られた2バンドの小領域画像から、ヘモグロビン成分を抽出する。まず、画像として観測された2バンドの信号は式(3.1)で示したように、様々な要素の積算で表されており、ここからヘモグロビンの成分だけを抽出することは容易ではない。そこで、この小領域画像の各画素値に対して対数を取ることで画像を濃度空間へ変換し、積算によって表されていた画素値を構成する要素の線形和で表現する。

次に、画像から濃度空間へ変換した点群に対し、新たな基底ベクトルを求める。照明の変動は濃度空間において強度のみの情報を持つことから、照明変動に関する成分は全て(1, 1)上に変換される。そのため、新たな基底ベクトルの第1ベクトルを照明変動の成分である(1, 1)に決定し、照明変動による成分を表しているため以降、陰影成分と呼ぶ。図７に示す濃度空間上にある点Aを、ヘモグロビン成分と陰影成分の新たな2成分に分離すると、点Aは以下の数式４（式(3.2)）のように表現することができる。

ここで、e_x=(1,0)、e_y=(0,1)でありx、 yは各波長の画素値の対数を取った値である。e_hとe_sはヘモグロビン成分と陰影成分を軸としたときの基底ベクトルでありx’、y’はe_hとe_sを軸としたときの点Aのヘモグロビン成分，陰影成分を表す。式(3.2)を行列で表すと数式５（式(3.3)）となる。

今、求めたい値はヘモグロビン成分と陰影成分の新たな2成分を軸とした基底変換後におけるx’、y’の値であるから、ヘモグロビン成分と陰影成分を要素とする基底変換行列の逆行列を計算することでヘモグロビンの色素濃度を推定することができる。また、本実施例では10×10の小領域画像を想定しているため、合計で100個の画素値が算出される。近赤外画像におけるヘモグロビンの変動値は1画素未満で現れるため、ヘモグロビンの成分を推定する際には整数値である画素値では変動を観測することはできない。従って、推定した100個のヘモグロビン成分に対してその平均値を取ることでヘモグロビン成分の変動を観測する。推定した100個のヘモグロビン成分を平均化した値は数式６（式（3.4））で表される。

ここで、iはi番目の画素値、Nは全体の画素数を表す。

今回のシミュレーションでは、ヘモグロビンが脈動によって変動する前の平均画素値と、脈動によってヘモグロビンの色素濃度が変化したあとの平均画素値の2つの値を算出し、3.3節で設定したヘモグロビン成分の変動が検出できるか判定する。

本シミュレータでは、ヘモグロビンの成分は脈動による変動の前と後の2段階を想定しているため、濃度空間でのヘモグロビン成分の差分をとることでヘモグロビンベクトルを取得した。陰影ベクトルは、濃度空間に変換した際に強度のみの情報を持つことになるため、常に(1,1)へと固定される。

3.6 ノイズを付加した際のヘモグロビン成分の推定
3.4節において、式(3.1)で示したようにカメラから出力される画素値は肌の反射率とカメラの感度、そしてバンドパスフィルタの感度の積算で表されると述べた。しかし実際の計測では、シャッターを閉じている状態でも微小な電流が流れる暗電流ノイズ等、ノイズが画素値に与える影響を想定することは非常に重要である。本節では、小領域画像にノイズを加えていったときに、どの程度のノイズまで加えてもヘモグロビン成分の抽出が可能か検証する。

ノイズを考慮した画素値Iは数式７で表される。

ここで、Rは肌の反射率、Cはカメラ感度、Bはバンドパスフィルタの感度、Sは照明強度、Nはノイズ、λは波長を表す。今回は、微小な電流と熱によって固定パターンノイズの原因となる暗電流ノイズを想定してNの値を増加する。ノイズNには正規分布から無作為に値を取り出す関数を使用し、増加していく定数と掛け合わせることで固定パターンノイズを作成した。また、極端に大きな値の乱数により、ノイズが画素値に対して支配的になることを防ぐため、ノイズを考慮した画素値の計算には、画素値の算出を100回行い、その平均を取ることで、各画素値に対する影響が一様になるようにノイズを調整してから1つの画素値を出力した。

第4章シミュレーション結果
4.1 はじめに
フィルタ感度とヘモグロビンの吸光度分布の関係から、2バンドのバンドパスフィルタを用いるとき、フィルタの波長が離れるにつれて酸化ヘモグロビンの吸光度特性に差が生じる。この各波長から取得できる吸光度特性の差を利用することで、1バンドでは取得できなかった酸化ヘモグロビンの特性を補正していく。本実施例で使用した酸化ヘモグロビンの吸光度分布は780nm~940nmの範囲を使用しており、この範囲においては波長が離れるにつれて吸光度特性の差が大きくなっていくことが図５からも確認できる。従って、波長が離れるにつれて多くの吸光度特性の差が取得できるため、ヘモグロビンの変動を推定しやすくなるという性質がある。一方、図５よりフィルタの波長が離れるにつれカメラの感度は下がり、ノイズによる画素値への影響が大きくなる、という性質がある。この2つの性質から、ヘモグロビンの吸光度分布とカメラの感度はトレードオフの関係にあるといえる。従って、ヘモグロビンの感度が十分異なり、ノイズから極端な影響を受けない程度にカメラの感度を保てる、最適なフィルタの組み合わせを選択することが必要である。

本章では、第3章で述べたシミュレータを用いて各フィルタの組み合わせを検証した結果について述べる。4.2節では、ノイズ毎のヘモグロビン成分の推定結果が、正しく分離されているかどうかを評価する分離評価関数について述べる。4.3節では、各フィルタにおけるシミュレーション結果について述べる。最後に4.4節で本章の考察を行う。

4.2 分離評価関数
この節では、3.5節、3.6節で述べたヘモグロビン成分の推定手法がどれだけノイズに対して耐性を持つか評価する関数について述べる。理想は、ヘモグロビン成分の分布から、脈動により変動したあとのヘモグロビンの分布を確認できることである。ノイズがない場合、この変動は容易に確認できるが3.6節で述べたようにノイズが画素値に加算されるにつれ、この変動が徐々にノイズによって小さくなっていく。また、求めたヘモグロビン成分がノイズの影響により大きな分散を持った場合も、ヘモグロビンの変動がノイズによって小さくなっていく、という2つの性質から、特定のノイズ量に対して推定したヘモグロビン成分がヘモグロビンの変動を保持しているか評価する分離評価関数Eは以下の数式８（式(4.1)）で表される。

ここで、I_h+Δhは脈動によるヘモグロビン変動後の推定値、I_hは変動前のヘモグロビンの推定値V_h+Δhはヘモグロビン変動後の分散、V_hは変動前のヘモグロビンの推定値の分散を表す。今回のシミュレーションはノイズを単調増加させていった時の評価値を求めたため、ノイズがないときの各フィルタの評価値は大きな値を取り、ノイズが増えるごとに評価値はゼロに近づいていく。このことから、ノイズが増えても評価値が最も高い値を保つバンドパスフィルタの組み合わせは、ノイズが増えても他のバンドパスフィルタの組み合わせより変化を検出することができることから、評価値が最も高いバンドパスフィルタの組み合わせを最適なフィルタの組み合わせとして最適化を行った。

4.3 シミュレーション結果
第3章で述べたシミュレーション環境で近赤外透過フィルタの最適化を行った。今回シミュレーションを行ったバンドパスフィルタの組み合わせは780[nm]と820[nm]、780[nm]と860[nm]、780[nm]と900[nm]、780[nm]と940[nm]の4通りである。各フィルタの組み合わせで2バンド近赤外画像を取得し、脈動によるヘモグロビン成分の変化前と変化後のヘモグロビン成分を推定した。これにノイズを徐々に加えていった際の分離評価関数による評価結果を図８に示す。ノイズの量を増加させた時の、推定したヘモグロビン成分は次第にノイズが支配的になるため、評価値は図８から分かるようにゼロに近づいていく。このことから、ノイズの増加に対して評価値が最も高い値を保っているのが780[nm]と900[nm]の組み合わせであることが分かる。

4.4 考察
4.3節で得られたシミュレーション結果より、どのフィルタに関しても4.2節で述べたようにノイズが増加するにつれ、ヘモグロビンの成分がノイズに埋もれていくため評価値はゼロに近づいていくことが分かる。しかし、各フィルタの組み合わせによっては、ノイズの大きさが変化することで評価値を大幅に減少しているフィルタの組み合わせや、評価値を保ったままのフィルタの組み合わせも存在する。これはシミュレータで用いた酸化ヘモグロビンの吸光度分布曲線、カメラ感度の分光曲線、バンドパスフィルタの分光曲線と、カメラの画像化システムの特性から決定されるためである。このことについて、各フィルタの組み合わせ毎の評価値について以下に述べる。

780[nm]と820[nm]のフィルタの組み合わせでは、評価値が最も早くゼロに近づいていることがわかる。これは、図６のヘモグロビンの吸光度分布曲線より、780[nm]と820[nm]ではヘモグロビンの吸光度分布に差がほとんどないため、3.5節で述べた濃度空間における画素値の分布がほぼ一直線になってしまっていると考えられる。そのため、新たに決定した陰影成分、ヘモグロビン成分の基底ベクトルが重なり、うまく分離が行えていないと考えられる。

780[nm]と860[nm]のフィルタの組み合わせでは、780[nm]と820[nm]の評価値と比較して、ノイズが増加してもヘモグロビンの変動が推定できていることが分かる。これは、780[nm]と820[nm]のフィルタの組み合わせの時と比べ、ヘモグロビンの吸光度分布がお互いに離れた感度を持っているためである。そのため、濃度空間における画素値の分布が広がりを持ち、陰影成分、ヘモグロビン成分に分離することができたと考えられる。

780[nm]と900[nm]のフィルタの組み合わせが、今回シミュレーションを行った中で最も高い評価値を持ったフィルタの組み合わせとなった。こちらも上と同様にヘモグロビンの吸光度分布が離れた感度を持っているため、濃度空間での画素値の分布が780nmと860nmの時よりさらに広がり、陰影成分、ヘモグロビン成分に分離することができたと考えられる。

780[nm]と940[nm]のフィルタの組み合わせでは、ノイズの増加に対して評価値が大きな影響を受け、780[nm]と860[nm]のフィルタの組み合わせよりも小さな評価値となっている。これは、図６から酸化ヘモグロビンの吸光度特性は僅かに変化しているが、その影響よりも図６(b)のシミュレーションに使用したモノクロカメラの感度が大きく影響している。シミュレーションに用いたモノクロカメラは近赤外光源に対しても感度を持っているが、波長が1000[nm]に向かうに従いカメラ感度がゼロに近づいていることが分かる。つまり、波長が1000[nm]に向かうにつれて出力される画素値は小さくなり、その分ノイズによる影響を受けやすくなる。ノイズに影響を受けやすくなるということは、僅かなノイズの増加によって、ヘモグロビン成分の変動が含まれた画素値はノイズが支配的となるため、僅かなノイズでも評価値が小さくなってしまう。

これらの結果より、本実施例でシミュレーションした結果では、酸化ヘモグロビンの吸光度特性が十分に離れており、カメラの感度もノイズに大きく影響されなかった780[nm]と900[nm]のフィルタの組み合わせが最適なフィルタの組み合わせであると考えられる。

第5章実測
5.1 はじめに
本章では、第4章で行った2バンドでの脈波計測シミュレーションの結果をもとに、実際の被験者に対して脈波測定を行った結果を述べる。5.2節では撮影環境を述べる。5.3節ではヘモグロビン成分推定手法について述べる。5.4節では脈波検出手法として、5.4.1節ではパワースペクトル密度推定、5.4.2節ではバンドパスフィルタ、5.4.3節ではピーク検出手法、5.4.4節では心拍数の推定手法について述べる。5.5節では推定値の精度検証手法について述べる。5.6節では被験者に対して測定を行った結果を述べる。最後に5.7節で本測定の考察を行う。

5.2 撮影環境
撮影環境を図９に示す。顎台で顔を固定した被験者に人工太陽灯を照射し、30[fps]に設定した2バンド近赤外カメラを用いて120[s]間、顔画像の撮影を行った。人工太陽灯は可視光源として用いられることが一般的であり、通常であれば紫外領域と近赤外領域の光はフィルタにより遮断されている。本実施例では近赤外領域の光源が必要なため、このフィルタを外して測定を行った。また、人工太陽灯は点灯させてから30~40分経過することで安定した照明状態となる。まず、十分に照明環の変動が起こらない場合での撮影を行った。次に、点灯直後の人工太陽灯の照明状態は不安定であり、照明にちらつきが発生することから、本実施例では、このちらつきを照明環境の変動であると仮定して測定を行った。このようにして、照明環境の変動がない場合とある場合の2種類の環境で測定を行った。被験者から人口太陽灯までの距離は0.30[m]、被験者からカメラまでの距離は0.50[m]である。被験者には、リラックスした状態でなるべく動かないよう指示した。また、人工太陽灯の光が直接目に入らないよう、安全面を考慮し測定中は目をつぶるよう指示した。カメラでの撮影と同時に、正解値として腕に電極を取り付け、心電図を測定した。心電図は5.7節で、推定した脈波信号との精度検証の際に用いる。撮影は暗室下で行われ、それぞれのモノクロカメラに入射する光は干渉型バンドパスフィルタを用いたことから、780[nm]を中心波長とした光と900[nm]を中心波長とした光のみとなるため、暗所での計測と同等の実験環境であるとして測定を行った。

5.3 ヘモグロビン成分の推定
本節では、撮影された2バンドの近赤外画像からヘモグロビン成分を推定し、陰影成分とヘモグロビン成分の2成分に分離する手法について述べる。

5.2節で述べた環境で撮影された、照明変動が含まれた2バンド近赤外画像を図１０に示す。図１０(a)は中心波長が780[nm]のバンドパスフィルタを用いた近赤外画像、図１０(b)は中心波長が900[nm]のバンドパスフィルタを用いた近赤外画像である。3.5節でも述べたように、撮影した画素値は肌におけるヘモグロビンの反射率、撮影に使用したモノクロカメラの感度、バンドパスフィルタの感度が積算されることにより表されているため、RGB空間からヘモグロビン成分だけ推定することは容易ではない。そこで、2バンド近赤外画像の各画素に対して対数をとることで、2バンド近赤外画像は濃度空間へ変換され、積算で表されていた画素値を構成する各要素はそれぞれの線形和で表される。このとき、照明の変動に関しては強度のみの情報を持つため濃度空間上では(1，1)ベクトルに固定される。実際の計測画像を濃度空間に変換した際の各画素値の分布を図１１に示す。図１１より(1,1)ベクトル方向での変動は照明変動による画素値の変化を表しており、これを以降陰影ベクトルと呼ぶ。一方、陰影ベクトルとは異なる他辺の変動がヘモグロビン成分による変動であり、この変動の向きに新たな基底ベクトルとなるヘモグロビンベクトルを取得することで、2バンドの近赤外画像からヘモグロビン成分と陰影成分を推定することができる。陰影成分は(1,1)で固定されていることから、ヘモグロビン成分のベクトルを決定することができれば、近赤外画像からヘモグロビン成分と陰影成分の推定が行える。

ヘモグロビンベクトルの推定手法は、まず濃度空間において1度ずつ基底ベクトルを回転させることで、360個のヘモグロビン成分の脈波信号を取得する。その中から、全ての脈波信号に対して後述する心拍数の推定を行い、その値が心電図で測定した正解値と誤差が最小になるベクトルをヘモグロビンベクトルとして決定した。

5.4 脈波検出手法
ヘモグロビン画像の額部分に関心領域を設定し、横軸を時間変動、縦軸を平均画素値として解析を行うことで脈波信号を取得することができる。今回は例として中心波長780nmの1バンド近赤外画像を解析することで得られた脈波信号を図１２に示す。この脈波信号には、被験者の僅かな体動やカメラのノイズ等の生体情報以外の様々な要素も含まれている。そのため、本節ではそれらのノイズを除去し、ヘモグロビン情報のみを抽出した脈波信号を推定する信号処理手法について述べる。

5.4.1 パワースペクトル密度推定
取得した脈波信号に対し傾き除去を行い、図１３に示す傾き除去後の脈波信号が得られた後、高速フーリエ変換を用いてパワースペクトル密度を推定することで、入力信号に含まれている周波数成分を解析した。脈波信号は周期的な信号であるため、そのパワースペクトル密度はある一定の周波数に強いピークを持つ。傾き除去後の脈波信号にパワースペクトル密度推定を行った結果が図１４であり、信号中に含まれる最も強い周波数が1.6Hzであることを示している。

5.4.2 バンドパスフィルタ
先行研究では、脈波信号からノイズを取得する際に成人の心拍数を45bpm~180bpm、即ち0.75Hz~3.0Hzを透過するバンドパスフィルタを適用することでノイズ除去を行った。一方、本手法では5.4.2節で行ったパワースペクトル密度推定より、0.75Hz~3.0Hzの間で最も強いパワースペクトルを持つ周波数f_maxに対して、[ f_max−0.2, f_max+0.2 ] Hzのみを透過させるバンドパスフィルタを適用することで、被験者の脈波信号に特化したバンドパスフィルタを設計した。さらに、バンドパスフィルタの透過帯を狭めることにより生体情報以外のノイズ除去精度が向上すると考えられる。図１５に、パワースペクトル密度推定から求めた最大ピーク周波数1.6Hzから、その近傍0.2Hzを透過するバンドパスフィルタを適用した結果を示す。

5.4.3 ピーク検出
バンドパスフィルタを適用して得られた脈波信号に対して近傍領域を設定し、その近傍領域内での最大値をピークとして検出する。例として、30fpsで撮影した、60bpmの心拍数を持つ被験者を例にすると、1秒間に1回のペースで脈動によるピークが発生するため、ピーク検出の際に用いる近傍領域は15フレームが適切となる。本手法では、パワースペクトル密度推定より算出された周波数と撮影環境が30fpsであることから、近傍領域の決定を適宜手動により調節しピーク検出を行った。この手法によりピーク検出を行った結果を図１６に示す。

5.4.4 心拍数の推定
ピーク検出を行った結果から、ピークとピークのフレーム差を算出することにより、脈波間隔と呼ばれる特徴量を取得することができる。先行研究では、この特徴量を周波数解析することにより精神状態の推定を行った。ピーク検出から求められる脈波間隔から、RR_{int erval} を脈波間隔としたときに、心拍数HR_estは以下の数式９（式(5.1)）で求められる。

5.5 精度検証手法
本実施例では、2バンド近赤外カメラによる撮影と同時に、接触型の心拍計を用いて脈波信号を計測した。心拍計はポリグラフ（日本光電：RMT-1000）を用いて計測を行い、得られた脈波信号から非接触式で利用した手法と同じ手法でピーク検出を行った。また、誤検出や検出漏れがないように、目視で波形とピークを確認し、閾値などの調整を行い、補正を行った。

また、カメラで計測される脈波信号は心臓が拍動してから顔に血流が届くまでにタイムラグがあるため、波形のピークが記録される時間の誤差を精度検証に用いることができない。そのため、(5.2)式で取得できる心拍数を用いて検証を行う。

カメラによる推定値と、心拍計から得られた脈波信号から、算出した正解値の心拍数を評価する手法として、絶対誤差率を用いた。絶対誤差率は数式１０（式(5.2)）で示されるように推定値と理論値の差を理論値で除算した百分率で表される。

ここで、HR_estは測定画像から推定した心拍数HR_gtは心電図から算出した心拍数を表す。絶対誤差率は推定値と理論値が近いほど値が小さくなり、その値がゼロに近づくほど推定精度が高いことを表す。

5.6 実験結果
照明変動のない環境で撮影した2バンド近赤外画像に対し、図１７のように額に関心領域を設定することで、関心領域中での平均画素値の時間変動を解析する。横軸に時間変動、縦軸に平均画素値として各フレームの平均画素値の変動として解析を行った、各近赤外波長における1バンド脈波信号を図１８に示す。

次に、照明変動のある時の環境で撮影された近赤外画像に対して、各1バンドの近赤外画像から脈波信号を取得した結果を図２０(a)、(b)に示す。その後、傾き除去、バンドパスフィルタ、ピーク検出を行った後に得られる心拍数推定に用いる脈波信号を図２１(a)、(b)に示す。

一方、撮影した2バンド近赤外画像から新たにヘモグロビンと陰影の基底ベクトルを取得し、ヘモグロビン画像と陰影画像に分離した結果を図２２に示す。本実施例ではヘモグロビン画像に着目して、被験者の額部分に関心領域を設定し、横軸を時間変動、縦軸を平均画素値として解析を行うことで、図２３(a)に示す、照明変動がない時のヘモグロビン成分における脈波信号を取得した。この脈波信号に対し先ほどと同様に傾き除去を行い、パワースペクトル密度推定により、0.75Hz~3.0Hzの間で最も強いパワースペクトルを持つ周波数に対して、[ f_max−0.2, f_max +0.2 ]Hzのみを透過させるバンドパスフィルタを適用した。このバンドパスフィルタを適用することで得られた脈波信号に対してピーク検出を行った結果を図２３(b)に示す。同様に、照明環境を変化させて撮影した画像からヘモグロビン成分の脈波信号と、その信号処理後の結果をそれぞれ図２４に示す。さらに、検出したピーク数から、脈波間隔を算出することで推定した心拍数と、心拍計より算出した正解値の心拍数を図２５に示す。このとき、780nm、900nmそれぞれの1バンド近赤外画像のみから取得した脈波信号より推定した心拍数も図２５に併記した。推定結果は照明環境の変動がある場合とない場合の2通りの結果を示した。また、この表から1バンドで測定した際の絶対誤差率と2バンドで測定した際の絶対誤差率を比較した結果を図２６に示す。

5.7 考察
照明環境が変化しない環境における心拍数の推定結果より、1バンドでの測定結果に関しては絶対誤差率が10%未満の精度で心拍数を推定することができたが、先行研究は全て1%台の精度で心拍数を推定しているため、本実施例の推定精度は低い。その原因として、顔の動きによる影響が考えられる。本実施例では測定の際に顎台を用いて顔を固定し測定を行っているが、実際には2分間の計測の中で被験者は意図せず顔の位置を動かしてしまう。先行研究では、顔の動き追跡は心拍数の推定に関して大きな影響を与える、と述べていることから被験者の顔の動きが推定精度の低下に影響を与えてしまったものと考えられる。一方、2バンド近赤外画像による心拍数の推定結果は、絶対誤差率が0.45%と、先行研究の1バンド近赤外画像から推定した心拍数よりも高精度で推定を行うことができた。これは、図２７に示した濃度空間における画素値が、陰影成分とヘモグロビン成分の2方向に分布していたため、適切に分離が行えたからだと考えられる。

照明変動のある環境で測定した1バンド近赤外画像から解析した脈波信号は、傾き除去とバンドパスフィルタを適用することで図２１(a)、(b)のような波形を取得した。しかし、心拍計から算出した心拍数との絶対誤差率はどちらも15%以上と、推定精度は低かった。これは、照明変動が周期的な信号ではないことから、パワースペクトル密度推定から決定した周波数帯にも照明変動に当たる周波数が含まれており、バンドパスフィルタだけでは照明変動によるノイズを除去することができなかったためだと考えられる。この結果からも分かるように、1バンドの近赤外画像では、照明の変動により心拍数の推定精度に大きな影響を与えている。

一方、照明変動がない時のヘモグロビン画像から推定した心拍数と正解値との絶対誤差率は3.73%と、1バンドによる測定に比べ、より高精度で心拍数の推定が行えていることが分かる。これは、5.3節でも述べたように濃度空間において照明変動はひとつの陰影ベクトルとして表され、ヘモグロビンベクトルは正解値と絶対誤差率の最も小さい値をとるベクトルに最適化したためである。しかし、先行研究と比較するとまだ十分な精度であるとは言えない。その原因の一つとして、2バンド近赤外画像からヘモグロビン画像と陰影画像に分離した図１８の画像より、ヘモグロビン画像の鼻の部分にまだ陰影が残っていることが確認できる。この結果より、濃度空間における陰影ベクトルとヘモグロビンベクトルの決定が最適でなかった可能性が挙げられる。図２７(b)より、撮影画像が濃度空間において陰影成分、ヘモグロビン成分方向に分布していることは確認できる。しかし、分離画像では正しく分離が行えていないことから、ベクトルの最適化手法に問題があると考えられる。本実施例で測定した被験者においては、濃度空間における2バンド画像の画素値の分布は第3象限において平行四辺形のように分布しており、それに対して原点からの基底ベクトルの最適化を行っていたため、分布に沿った適切なベクトルに決定することができなかった。これはバイアスベクトルを与えることにより原点を移動させてから基底ベクトルを決定することで改善すると考えられる。また、ヘモグロビンベクトルが(1, 1)と(-1, -1)の時の際は基底変換行列の逆行列が計算できず、今回のシミュレーションで想定していなかったため、今後このベクトルも含めて最適化を行う。また、照明変動がない時の画像から求めたヘモグロビンベクトルの値を使って分離を行うことで、現時点よりもより高精度に心拍数を推定できる可能性があると考えられる。

第6章結論
6.1 本実施例のまとめ
本実施例では、2バンドの近赤外画像を用いることで照明変動があるような実環境を想定した暗所での脈波信号の検出手法を提案した。ヘモグロビンベクトルに関しては、決定した陰影ベクトルを基に、1度ずつ基底ベクトルを回転させて心拍数まで推定することにより、正解値と比較して、最も小さい誤差率となるベクトルを、ヘモグロビンの変動成分として決定した。こうして濃度空間上で新たな基底ベクトルを取得することで、2バンド近赤外画像をヘモグロビン画像と陰影画像に分離した。これにより得られたヘモグロビン画像の時間変動を解析することによりヘモグロビン成分の脈波信号を取得した。これを、先行研究でも挙げられているように傾き除去、バンドパスフィルタ、ピーク検出による信号処理を行うことで脈波間隔という特徴量を取得し、心拍数の推定を行った。1バンドでの脈波信号に比べ、2バンドでの脈波信号から推定した心拍数がより正解値との推定誤差率が小さいことから、2バンド近赤外画像を用いることで照明の変動による心拍数の推定精度への影響を大幅に軽減して心拍数の推定を行うことができた。

本発明は、心拍変動の推定方法、心拍変動の推定プログラム及び心拍変動推定システムとして産業上の利用可能である。

Claims

対象者の少なくとも２種類の波長域の動画像を撮影し、前記撮影した動画フレーム中の各画素の対数をとることで変換される少なくとも２つの基底ベクトルで定められる対数空間上において、第１ベクトルとして新たに決定された照明変動成分と、前記第１ベクトルとは異なる角度のベクトルに新たに決定されたヘモグロビン成分の少なくとも２成分に行列の基底変換を行い、前記対数変換された画素データを前記基底変換して得られたヘモグロビン成分の情報に基づいて心拍変動の推定を行う心拍数推定方法。
前記動画像は波長が1000nmに向かうに従い感度がゼロに近づくカメラで撮像されており、
前記少なくとも２種類の波長域のうち、１種類の波長域が770nm以上790nm以下であり、
他の１種類の波長域が890nm以上910nm以下である請求項１記載の心拍変動推定方法。
コンピュータに、対象者の少なくとも２種類の波長域の動画像を撮影する手順と、前記撮影した動画フレーム中の各画素の対数をとることで変換される少なくとも２つの基底ベクトルで定められる対数空間上において、第１ベクトルとして新たに決定された照明変動成分と、前記第１ベクトルとは異なる角度のベクトルに新たに決定されたヘモグロビン成分の少なくとも２成分に行列の基底変換を行い、前記対数変換された画素データを前記基底変換して得られたヘモグロビン成分の情報に基づいて心拍変動の推定を行う手順とを実行させるための心拍数変動推定プログラム。
対象者の少なくとも２種類の波長域の動画像を撮影するカメラと、前記カメラにより撮影した動画フレーム中の各画素の対数をとることで変換される少なくとも２つの基底ベクトルで定められる対数空間上において、第１ベクトルとして新たに決定された照明変動成分と、前記第１ベクトルとは異なる角度のベクトルに新たに決定されたヘモグロビン成分の少なくとも２成分に行列の基底変換を行い、前記対数変換された画素データを前記基底変換して得られたヘモグロビン成分の情報に基づいて心拍変動の推定を行う心拍変動推定部と、を備える心拍数変動推定システム。