JP2015097638A - 生体情報測定方法 - Google Patents

Abstract

【課題】不要な振動ノイズが多い状況下においても、心拍信号等の生体情報を正確に測定することができる生体情報測定方法を提供する。
【解決手段】Nを3〜5程度の自然数とし、時刻tにおいてドップラーレーダモジュールが出力するIQ直交２信号をVI(t)とVQ(t)とすると、これらを全時刻tにおいてN+1次多項式を用いて最小二乗法近似してから、N階時間微分を求め、高階微分結果によって得られた振幅ピーク位置として心拍信号等の生体情報を抽出する。
【選択図】図１

Description

本発明は、ドップラーレーダ波を人体の胸部に当てることによって胸部表面位置の変化を測定し、その情報の中から微細な心拍信号や呼吸信号等の生体情報を抽出する生体情報測定方法に関するものである。

近年の家電業界や医療業界では、「ヘルスケア」が流行語となり、大きな注目を集めている。「健康管理」とでも訳すべき「ヘルスケア」という概念の根底には、対処療法中心の医療から予防医療へのシフトが必要であるという共通認識がある。その予防医療を普及させる前提として、家庭においてもオフィスにおいても通勤途中の車の中でも、２４時間常に人間の心拍・呼吸・血圧・体温等の生体情報をモニターしてビッグデータ化する時代の到来が予想されている。特に長時間運転者を車中で連続モニターする場合は、疲労度や居眠り予兆を検知することによって、単なる「健康管理」の域を超えて事故防止にまで繋がるメリットも期待される。

ところで、生体情報を２４時間モニターする場合には、人体にケーブル等の測定器具を装着して人体を拘束することは許されないため、非接触や非侵襲の測定技術が必須となる。特に、車の中で運転者の生体情報をモニターする場合は、人体を拘束することは事故に直結するために尚更である。

そのような非接触の測定技術の１つとして、ミリ波帯電波を用いたドップラーレーダ方式が注目されている。ここで、ミリ波帯電波としては、ARIB規格STD-T73で規定された24.15GHz帯が有望視されている。ドップラーレーダを用いた非接触心拍測定方法については、例えば特許文献１や非特許文献１，２等の多数の文献において開示されている。

実際の測定系の構成を図８の模式図に示す。

図８に示すドップラーレーダモジュール１０２においては、発振回路１１１が発生した高周波正弦波は、点１１６において２経路に分割される。一方は送信アンテナ１１３から送信され、人体１０１の胸部に当たって反射し、受信アンテナ１１４から取り込まれる。受信信号と元の信号はミキサ回路１１２にて位相比較され、その位相差に相当する低周波IF出力が、出力端子１１５から信号処理部１０３へと出力される。この位相差は、アンテナ１１３，１１４と人体胸部との間の距離に略比例した情報が含まれており、これを信号処理部１０３で詳細解析することによって距離の変動を逆算することができ、そこから胸の表面の運動を推定することができる。尚、ミキサ回路１１２として安価なシングルミキサを用いた場合には、出力端子１１５は１本しかないが、バランス型直交ミキサを用いた製品では、IQ直交2信号が出力されるために出力端子１１５は2本になる。

ドップラーレーダモジュール１０２は、現在では安価な市販品を用いるのが一般的であり、多くの場合、性能を決定するのは主に信号処理部１０３となる。出力端子１１５から得られた実際の信号のスペクトルの一例を図９に示す。

図９に示すように、出力端子１１５から得られる信号には、通常は複数種類の生体信号やノイズが混在してしまい、どれが所望の生体信号なのか判別できない場合が多い。24.15GHzの電波の波長は1.2cmもあり、数mm程度の大きな振幅を持つ呼吸運動は検出し易いが、特に心拍運動の振幅は僅かμmオーダであり、もともと単純計算ではノイズに埋もれてしまう。そのため、従来技術においては、微細な心拍信号を抽出する方法として、周波数領域で線形バンドパスフィルタを掛けて1〜1.5Hz付近の信号のみ取り出す方法を採用している。このようなアプローチの代表例として、特許文献１にはフィルタ性能向上の方法の一例が詳細に開示されている。

特開２０１３−１５３７８２号公報 松井岳巳、「高性能レーダーを用いた自動車運転時の非接触バイタルサインモニタリング」、科学研究費補助金研究成果報告書、課題番号19519174、平成21年5月20日 川本貴志、「心拍由来体表面振動を用いた低筋力作業時筋疲労解析」」、慶應義塾大学大学院理工学研究科2012博士学位論文

しかしながら、従来のドップラーレーダによる生体情報測定方法は、振動ノイズに極めて弱いという問題があった。ここで言う振動とは、身体全体の体動や、心拍運動に対する呼吸運動や、車載の場合はエンジンや路面からの振動等、様々なものが考えられる。

従来技術が振動ノイズに極めて弱い原因は、特に微細な心拍信号を測定する場合には、測定したい心拍信号を単純に周波数領域で1〜1.5Hz付近に存在する信号としてしか識別できていないためである。体動や車の振動は、ランダムな周波数成分を含むために頻繁1〜1.5Hz付近に重なってくる。又、呼吸運動についても、呼吸運動自体は一般的に0.25〜0.5Hz付近が多いと言われているが、意識的・無意識的に幾らでも変化するものであり、その高調波成分まで含めれば心拍信号の1〜1.5Hz付近に重なる事態は避けられない。呼吸運動の高調波成分による妨害問題については、特許文献１で詳しく議論されている。従来技術のように周波数領域の線形バンドパスフィルタに依存した技術では、いくらフィルタ性能を強化しようとも、根本的に周波数成分が重なってしまった場合の妨害は原理的に回避不能である。

本発明は上記問題に鑑みてなされたもので、その目的とする処は、不要な振動ノイズが多い状況下においても、生体情報、特に微細な心拍を正確に測定することができる生体情報測定方法を提供することにある。

上記目的を達成するため、請求項１記載の発明は、ドップラーレーダ波を人体の胸部に当てることによって胸部表面位置の変化を測定し、その測定情報の中から微細な心拍信号等の生体情報を抽出する生体情報測定方法において、
コンピューターが、
放射波と反射波の位相差変化として得られるIQ直交２信号に対して、Nを3〜5程度の自然数とし、時刻tにおいてドップラーレーダモジュールが出力するIQ直交２信号をVI(t)、VQ(t)とすると、これらのVI(t)、VQ(t)を全時刻tそれぞれにおいてN+1次多項式：
f(t)=a_(N+1) * t^(N+1) + a_N * t^N + ・・・ + a_1 * t + ａ_0
を用いて最小二乗法近似してからN階時間微分による、N=3〜5程度の高階微分を求め、当該高階微分結果によって得られた振幅ピーク位置として心拍信号等の生体情報を抽出することを特徴とする。

請求項２記載の発明は、請求項１記載の発明において、前記コンピューターが、前記高階微分結果を、更に円周座標系へ変換することによって、円周方向Ｎ階微分Vφ_N(t)と半径方向Ｎ階微分Vr_N(t)を求め、求められたVφ_N(t)とVr_N(t)と、周波数fが測定したい生体情報信号よりも遥かに高く設定されたMorlet関数或いは類似の局所的周期関数との時間領域相関係数を計算し、得られた複素相関係数のピーク位置として心拍信号等の生体情報信号を得ることを特徴とする。

請求項３記載の発明は、請求項２記載の発明において、前記Morlet関数の次数mが2〜10、周波数fが7〜25Hzであることを特徴とする。

請求項４記載の発明は、請求項１記載の発明において、コンピューターが、Mを自然数として、M階微分結果が成すベクトル（ VI_M(t), VQ_M(t) ）とM±1階微分結果が成すベクトル（ VI_(M±1)(t),
VQ_(M±1)(t) ）からはベクトル外積により、M階微分結果が成すベクトル（ VI_M(t), VQ_M(t) ）とM±2階微分結果が成すベクトル（ VI_(M±2)(t),
VQ_(M±2)(t) ）からはベクトル内積により、時間領域のピーク位置として心拍信号等の生体情報信号を得ることを特徴とする。

本発明によれば、周波数領域の線形処理ではなく、時間領域において非線形処理を行うようにし、具体的には、先ず時間波形に対して3〜5階の高階微分を行うことによって特徴的な信号成分のみを浮き上がらせ、次にその信号成分と形が良く似た局所的周期関数との相関積分を求めることによって更に信号成分を浮き上がらせるようにしたため、不要なノイズ成分については時間領域の高階微分と相関成分の２つの作用によって減衰し、たとえ周波数領域で見ると所望信号と重なっているように見えても、不要なノイズ成分を除去することができ、不要な振動ノイズが多い状況下においても、微細な心拍信号等の生体情報を確実に抽出して正確に測定することができる。

発明の実施の形態１に係る生体情報測定方法によって心拍信号を抽出する処理手順を示すフローチャートである。 HSV波形の一例を示す図である。 Morlet関数の一例を示す図である。 本発明の実施の形態１に係る生体情報測定方法において得られる心拍信号の時間波形を示す図である。 本発明の実施の形態２に係る生体情報測定方法によって心拍信号を抽出する処理手順を示すフローチャートである。 本発明の実施の形態２に係る生体情報測定方法において得られる心拍信号の時間波形を示す図である。 接触測定と被接触測定によって得られた心拍数の相関を示す図である。 従来技術及び本発明共通の測定系の構成を示す模式図である。 ドップラーレーダモジュールの出力の一般的な周波数スペクトルの模式図である。

以下に本発明の実施の形態を添付図面に基づいて説明する。
［実施の形態１］
ホン発明の実施の形態１に係る生体情報測定方法によって心拍測定を実施するための測定系の構成は、図８に示したものと同じであるため、これについての再度の説明は省略し、以下の説明では図８に示した符号をそのまま使用する。尚、本発明は、特に微細な心拍振動を抽出することにおいて顕著な効果を発揮するため、以下の実施の形態においては、心拍測定方法について説明する。

本発明方法においては、ドップラーレーダモジュール１０２の出力１１５としてIQ直交2信号が得られるものを用いるが、このようなドップラーレーダモジュール１０２は例えば独InnoSent社のIPS-265等が簡単且つ安価に入手できる。

本発明の要点は、コンピューター等のハードウェア内に設けられる信号処理部１０３にあり、この信号処理部１０３において下記の時間領域処理を行うことを特徴とする。以下、この処理を図１に示すフローチャートに従って説明する。

本実施の形態に係る心拍測定方法は、ドップラーレーダ波を人体の胸部に当てることによって胸部表面位置の変化を測定し、その測定情報の中から微細な心拍信号を抽出する方法であって、放射波と反射波の位相差変化として得られるIQ直交２信号に対して、Nを3〜5程度の自然数とし、時刻tにおいてドップラーレーダモジュール１０２が出力するIQ直交２信号をVI(t)、VQ(t)とする（図１のステップＳ１）。そして、これらのVI(t)、VQ(t)を全時刻tにおいてN+1次多項式：
f(t)=a_(N+1) * t^(N+1) + a_N * t^N + ・・・ + a_1 * t + ａ_0 …（１）
を用いて最小二乗法近似する（ステップＳ２，Ｓ３）。

次に、上記（１）式によって求められるf(t)をN階時間微分する。具体的には、N=3〜5程度の高階微分を行うことによって、VI(t)、VQ(t)のＮ回微分の近似値VI(t)_N(t)、VQ(t)_N(t)をそれぞれ求め（ステップＳ４，Ｓ５）、その結果が成す直交座標ベクトル（I(t)_N(t)、VQ(t)_N(t)を円周方向φ、半径方向ｒから成る円周座標ベクトル（Vφ_N(t)、Vr_N(t)）に変換する（ステップＳ６）。

そして、求められたVφ_N(t)とVr_N(t)と、周波数fが測定したい心拍信号よりも遥かに高く設定されたMorlet関数或いは類似の局所的周期関数との時間領域相関係数を計算する。本実施の形態では、時間領域で、Vφ_N(t)、Vr_N(t)とMorlet関数を相関積分することによって時間領域相関係数を計算する（ステップＳ７，Ｓ８）。この場合、Morlet関数は、次数mが2〜10、周波数fが7〜25Hzであるとする。

上述のようにして時間領域相関係数が計算によって求められると、相関係数のピーク位置を心拍信号として抽出し、この抽出された心拍信号によって心拍を測定する（ステップＳ９）。

ここで、本発明に係る心拍測定方法の理論的根拠について説明する。

図２は前述の非特許文献２に記載された「心拍起因の体表面振動（HSV）」と呼ばれるものの波形を示す図である。このように体内深部の心臓が振動することによって、体表面でも微細な振動が起こる。このHSV波形の特徴は、図２の細かいギザギザ（I音及II音）で示されている通り、高周波成分を極めて多く含んでいることである。これは、心臓自体が高周波振動する訳ではないが、複数ある心臓弁が少しずつタイミングをずらして開閉する結果として、あたかも高周波振動が発生したかのように見えるためである。その結果として、心拍運動自体（1〜1.5Hz）は人間の可聴域（20Hz〜20kHz）よりも低い成分であるにもかかわらず、医者が聴診器を胸に当てると、HSVの高周波成分という形で可聴域の心音を聴くことができる。本発明は、ドップラーレーダにおいてこの医者の聴診器と同じ原理の応用を目指したものである。

本発明では、図２の高周波成分（I音及II音）が、ウェーブレット変換で良く使われるような局所的周期関数に形が似ていることに注目した。通常のフーリエ変換は、時間領域で無限に広がる周期関数cos(2πft)+i・sin(2πft) との相関積分を行うものである。これに対して、ウェーブレット変換における局所的周期関数は、基本的に前記周期関数に対して時間領域の窓関数を掛けて、時間領域における広がりを局所的に制限したものである。ここで、iは虚数単位、fは周波数、tは時間を指す。

局所的周期関数には何種類もあるが、代表的なものとしてMorlet関数を図３に示す。Morlet関数は、窓関数をexp{-(2πft/m)^2}
としたものであり、自然数mは次数と呼ばれるパラメータであって、図３中の突起の数を増減させる。図３（ａ）はm=4、f=10Hzとした場合の例である。このMorlet関数には重要な性質があり、時間tに対して微分を何回繰り返しても、周期関数部分の位相が毎回90度回転するだけで、波形そのものは変化しない。図３（ｂ）は図３（ａ)を１階微分した場合の例である。

本発明に係る心拍測定方法の原理は、上記理解に基づいている。先ず、時間領域で高階微分を行うことによって、微分の度に周期性が無いノイズ成分は減衰して行き、周期性があるMorlet関数成分だけが残って浮き上がる。そして、最後に、残ったMorlet関数類似成分に対してMorlet関数そのものとの相関積分を行うことによって、残ったノイズ成分も排除することができる。その際、Morlet関数の周波数fは、図２の高周波成分（I音及II音）を見れば分かる通り、測定したい心拍信号そのもの（1〜1.5Hz）よりも遥かに高い周波数に設定しておく。実験にて試行錯誤した結果、f=7〜25Hz、m=2〜10の範囲が適当なようであった。最終的に心拍情報は、相関係数のピーク位置として得ることができる。

尚、一般的には、ノイズを含む測定データに対して３〜５階もの高階微分を行うことは、ノイズが著しく増幅発散するため、行われていなかった。しかし、本発明においては、上記特性に着目して、全時刻ｔそれぞれにおいて測定データではなく、多項式で最小二乗近似することによってノイズを除去し、測定データ自体ではなく、この近似多項式を微分することによって、高階微分にもかかわらず高精度を確保することに成功した。

次に、本発明方法における信号処理の結果として得られた各段階の時間波形の一例を図４（ａ）〜（ｅ）に示す。

図４（ａ）は図１のA点で得られた時間波形を示し、この時間波形は、ドップラーモジュールのIQ直交2信号であるVI(t)とVQ(t)そのものである。この波形を眺めても、心拍信号は余りにも微細過ぎてどこに存在するのか肉眼では識別できない。

４（ｂ）は図１のB点で得られた時間波形を示し、VI(t)とVQ(t)を1階微分した結果であるVI_1(t)とVQ_1(t)である。VI_1(t)とVQ_1(t)そのままでは分かりにくいため、これの円周方向φ成分を計算して表示してある。この段階では、約4秒周期の大きな呼吸運動は確認できるが、約1秒周期の微細な心拍運動はやはり識別困難である。

図４（ｃ）は図１のC点で得られた時間波形を示し、VI(t)とVQ(t)を3階微分した結果であるVI_3(t)とVQ_3(t)である。VI_3(t)とVQ_3(t)そのままでは分かりにくいため、これの円周方向φ成分を計算して表示してある。本発明者の経験によれば、このような3階微分の頃からノイズが消えて心拍運動が明瞭に見えてくる。

図４（ｄ）は図１のD点で得られた時間波形を示し、VI(t)とVQ(t)を5階微分した結果であるVI_5(t)とVQ_5(t)である。VI_5(t)とVQ_5(t)そのままでは分かりにくいため、これの円周方向φ成分を計算して表示してある。本発明者の経験によれば、3階微分を超えた頃からMorlet関数に類似した成分が主体になり、それ以上何回微分を繰り返しても波形は余り変化しなくなる。

図４（ｅ）は図４（ｄ）とMorlet関数との時間領域相関積分を取った結果であり、縦軸は相関係数を表している。前述のように、Morlet関数の周波数fとしては、測定したい心拍運動の周波数（1〜1.5Hz）ではなく、敢えて遥かに高い周波数（7〜25Hz）に設定している。これにより、図５(ｅ）に示すように、矢印↓で示した位置においてピークが発生し、１秒弱の周期の心拍運動を抽出できていることが確認できる。
[実施の形態２]
次に、本発明の実施の形態２を図５に示すフローチャートに従って説明する。

本実施の形態も前記実施の形態１と同様に、ドップラーレーダ波を人体の胸部に当てることによって胸部表面位置の変化を測定し、その測定情報の中から微細な心拍信号を抽出する心拍数測定方法であって、前記実施の形態１において説明した、VI(t)、VQ(t)のＮ回微分の近似値VI(t)_N(t)、VQ(t)_N(t)を求める処理（図１のステップＳ１〜Ｓ５）までは同じである（図５のステップＳ１１〜Ｓ１５）。

而して、本発明に係る心拍測定方法においては、VI(t)、VQ(t)のＮ回微分の近似値VI(t)_N(t)、VQ(t)_N(t)が求められると、Mを自然数として、M階微分結果が成すベクトル（ VI_M(t), VQ_M(t) ）とM±1階微分結果が成すベクトル（ VI_(M±1)(t),
VQ_(M±1)(t) ）からはベクトル外積を取り、M階微分結果が成すベクトル（ VI_M(t), VQ_M(t) ）とM±2階微分結果が成すベクトル（ VI_(M±2)(t),
VQ_(M±2)(t) ）からはベクトル内積を取る（ステップＳ１６）。そして、外積・内積の時間領域ピーク位置を心拍信号として抽出し、この抽出された心拍信号によって心拍を測定する（ステップＳ１７）。

本実施の形態に係る心拍測定方法においても、N階微分が終わるまでの処理は前記実施の形態１における心拍測定方法と同じであるため、その処理の結果得られる時間波形は図４（ａ）〜（ｄ）に示すものと同じである。

而して、本実施の形態に係る心拍測定方法においては、前述のように高階微分結果が成すベクトル同士の外積や内積を計算し、そのピークの形で心拍信号を抽出するようにしたため、最終的には図６に示す時間波形が得られる。図６に示す時間波形は、外積については5階微分結果ベクトル(
VI_5(t)、 VQ_5(t) )と4階微分結果ベクトル( VI_4(t)、 VQ_4(t)
) を掛けたもの、内積については5階微分結果ベクトル(
VI_5(t)、 VQ_5(t) )と3階微分結果ベクトル( VI_3(t)、 VQ_3(t)
) を掛けたものある。この時間波形においても、図６において矢印↓で示した位置においてピークが発生し、１秒弱の周期の心拍運動を抽出できることが確認できる。

ここで、アイドリング状態の車の中で心拍を測定した場合の一例を図７に示す。

図７（ａ）は従来技術である非特許文献１の図４（Ａ）を掲載したものであり、横軸は正確な接触型心電図計による心拍測定結果であり、縦軸は本発明と同様にドップラーレーダによる非接触測によって得られた心拍測定結果である。図７（ａ）に示す結果によれば、車がアイドリング状態で振動ノイズが存在する場合、非接触測定結果が不正確にばらつくことが分かる。

それに対して、図７（ｂ）は第１発明に係る測定方法を用いた場合の結果である。実験は安全確保のため助手席で行い、ドップラーレーダモジュールは胸の真正面のダッシュボード上に取り付けた。図７（ｂ）に示す結果によれば、車がアイドリング状態で軽い振動ノイズが存在しても、測定精度が劣化しないことが確認できる。

以上のように、本発明によれば、周波数領域の線形処理ではなく、時間領域において非線形処理を行うようにし、具体的には、先ず時間波形に対して3〜5階の高階微分を行うことによって特徴的な信号成分のみを浮き上がらせ、次にその信号成分と形が良く似た局所的周期関数との相関係数を求めることによって更に信号成分を浮き上がらせるようにしたため、不要なノイズ成分については時間領域の高階微分と相関成分の２つの作用によって減衰し、たとえ周波数領域で見ると所望信号と重なっているように見えても、不要なノイズ成分を除去することができ、不要な振動ノイズが多い状況下においても、心拍信号を確実に抽出して心拍を正確に測定することができる。

尚、上記実施の形態１，２においては、心拍測定方法について説明したが、心拍のような所定間隔でピークが現れる所謂心電図の鋭いピークであるＲ波を利用した、Ｒ波の発生時刻と、１つ前のＲ波の発生時刻の時間差、間隔を計測する｛ＲＲ間隔（RR Interval）の特性を有するものの測定についても本発明を有効に適用することができる。

又、上記ＲＲ間隔特性を有した心拍信号以外でも、図２に示したＨＳＶ波形と類似した細かい振動成分を持つその他の生体情報の測定に対しても本発明を適用することができる。

１０１ 人体
１０２ ドップラーレーダモジュール
１０３ 信号処理部
１１１ 発振回路
１１２ ミキサ回路
１１３，１１４ アンテナ
１１５ 出力端子

Claims (4)

1. ドップラーレーダ波を人体の胸部に当てることによって胸部表面位置の変化を測定し、その測定情報の中から微細な心拍信号等の生体情報を抽出する生体情報測定方法において、
   コンピューターが、
   放射波と反射波の位相差変化として得られるIQ直交２信号に対して、Nを3〜5程度の自然数とし、時刻tにおいてドップラーレーダモジュールが出力するIQ直交２信号をVI(t)、VQ(t)とすると、これらのVI(t)、VQ(t)を全時刻tそれぞれにおいてN+1次多項式：
   f(t)=a_(N+1) * t^(N+1) + a_N * t^N + ・・・ + a_1 * t + ａ_0
   を用いて最小二乗法近似してからN階時間微分による、N=3〜5程度の高階微分を求め、当該高階微分結果によって得られた振幅ピーク位置として心拍信号等の生体情報を抽出することを特徴とする生体情報測定方法。
2. 前記コンピューターが、前記高階微分結果を、更に円周座標系へ変換することによって、円周方向Ｎ階微分Vφ_N(t)と半径方向Ｎ階微分Vr_N(t)を求め、求められたVφ_N(t)とVr_N(t)と、周波数fが測定したい生体情報信号よりも遥かに高く設定されたMorlet関数或いは類似の局所的周期関数との時間領域相関係数を計算し、得られた複素相関係数のピーク位置として心拍信号等の生体情報信号を得ることを特徴とする請求項１記載の生体情報測定方法。
3. 前記Morlet関数の次数mが2〜10、周波数fが7〜25Hzであることを特徴とする請求項２記載の生体情報測定方法。
4. コンピューターが、Mを自然数として、M階微分結果が成すベクトル（ VI_M(t), VQ_M(t) ）とM±1階微分結果が成すベクトル（ VI_(M±1)(t),
   VQ_(M±1)(t) ）からはベクトル外積により、M階微分結果が成すベクトル（ VI_M(t), VQ_M(t) ）とM±2階微分結果が成すベクトル（ VI_(M±2)(t),
   VQ_(M±2)(t) ）からはベクトル内積により、時間領域のピーク位置として心拍信号等の生体情報信号を得ることを特徴とする請求項１記載の生体情報測定方法。
   

Images