(57)【特許請求の範囲】
【請求項１】
圧力変化を測定する２個以上のセンサ群を備え、１個以上の加速度計に接続された支持部材上に存在する対象の生体パラメータを推定するための推定方法であって、
前記センサ群により測定された圧力変化に対応する第１信号に基づいて、前記センサ群に含まれるセンサそれぞれにセンサ固有モデルを提供するステップと、
前記加速度計からの第２信号に基づいて、前記対象の生体パラメータ推定に使用する１つの選択センサを、タイムフレームごとに前記センサ群から選択する選択ステップと、
前記選択センサに対応する前記センサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定ステップと、
を含み、
前記選択ステップにおいて、
帰無仮説パラメータ、確率密度関数パラメータ、カルバック・ライブラー情報量パラメータ、歪度パラメータおよび尖度パラメータの少なくともいずれか１つを含む統計パラメータおよびｑ−ハーストパラメータの少なくともいずれか一方を含むパラメータであって、前記第１信号および前記第２信号のいずれか一方からタイムフレームごとに抽出されたパラメータを備えた特徴ベクトルを算出し、
前記特徴ベクトルを非線形モデルに入力して、前記特徴ベクトルのパラメータを、タイムフレームごとに、前記選択センサの番号に割り当てる生体パラメータ推定方法。
【請求項２】
前記非線形モデルは、非線形関数と確率に基づいた決定とを組み合わせることにより、前記特徴ベクトルのパラメータを、タイムフレームごとに、前記選択センサの番号に割り当てる請求項１に記載の生体パラメータ推定方法。
【請求項３】
圧力変化を測定する２個以上のセンサ群を備え、１個以上の加速度計に接続された支持部材上に存在する対象の生体パラメータを推定するための推定方法であって、
前記センサ群により測定された圧力変化に対応する第１信号に基づいて、前記センサ群に含まれるセンサそれぞれにセンサ固有モデルを提供するステップと、
前記加速度計からの第２信号に基づいて、前記対象の生体パラメータ推定に使用する１つの選択センサを、タイムフレームごとに前記センサ群から選択する選択ステップと、
前記選択センサに対応する前記センサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定ステップと、
を含み、
前記推定ステップは、前記センサ群のいずれに対しても、前記第１信号のそれぞれに基づいて、状態ベクトルの、初期周波数を含む状態パラメータを推定する初期化ステップを含み、
前記初期化ステップにおいて、前記センサ群のそれぞれに対して、前記状態パラメータごとに、共分散値をプロセスノイズ値にマッピングする非線形マッピングモデルを利用して、プロセスノイズ値を取得し、取得したプロセスノイズ値を前記推定ステップに利用し、
非線形マッピングモデルの中間層のバイアスベクトルと前記共分散値より計算した変換値を乗じた前記中間層の係数ベクトルとのシグモイド関数を計算する第１計算ステップと、
前記非線形マッピングモデルの出力層のバイアスベクトルと、前記第１計算ステップの結果を乗じた前記出力層の係数ベクトルとのマッピング関数を計算する第２計算ステップと、
を実行することにより、前記マッピングモデルが、共分散値をプロセスノイズ値に割り当てる生体パラメータ推定方法。
【請求項４】
前記選択ステップにおいて、前記第１信号および前記第２信号を用いて、前記選択センサを選択する請求項１、２または３に記載の生体パラメータ推定方法。
【請求項５】
前記推定ステップは、前記センサ群のいずれに対しても、前記第１信号のそれぞれに基づいて、状態ベクトルの、初期周波数を含む状態パラメータを推定する初期化ステップを含む請求項１または２に記載の生体パラメータ推定方法。
【請求項６】
前記初期化ステップにおいて、前記センサ群のそれぞれに対して、前記状態パラメータごとに、共分散値をプロセスノイズ値にマッピングする非線形マッピングモデルを利用して、プロセスノイズ値を取得し、取得したプロセスノイズ値を前記推定ステップに利用する請求項５に記載の生体パラメータ推定方法。
【請求項７】
非線形マッピングモデルの中間層のバイアスベクトルと前記共分散値より計算した変換値を乗じた前記中間層の係数ベクトルとのシグモイド関数を計算する第１計算ステップと、
前記非線形マッピングモデルの出力層のバイアスベクトルと、前記第１計算ステップの結果を乗じた前記出力層の係数ベクトルとのマッピング関数を計算する第２計算ステップと、
を実行することにより、前記マッピングモデルが、共分散値をプロセスノイズ値に割り当てる請求項６に記載の生体パラメータ推定方法。
【請求項８】
前記推定ステップは、前記センサ群のいずれに対しても、前記選択ステップにおいて計算された、センサごとのセンサ固有モデルを用いて、前記信号の前サンプル（ｋ−１）の状態ベクトルに基づいて、前記信号の現サンプル（ｋ）の状態ベクトルを算出する状態ベクトル算出ステップを含む請求項５、６または７に記載の生体パラメータ推定方法。
【請求項９】
前記状態ベクトル算出ステップにおいて、前記センサ群の状態パラメータごとに得られたプロセスノイズ値を用いて、状態ベクトルを算出する請求項８に記載の生体パラメータ推定方法。
【請求項１０】
前記状態ベクトル算出ステップは、
前記前サンプルの前状態ベクトルとセンサごとのセンサ固有モデルの線形モデルとに基づいて前記現サンプルの現状態ベクトルを予測する予測ステップと、
予測された前記現状態ベクトルを、前記センサごとのセンサ固有モデルの非線形モデルを用いて更新する更新ステップと、
前記選択ステップにおいて計算された１個のセンサに基づいて、生体パラメータ推定のために用いるべきセンサ固有モデルを、前記選択センサに提供されるセンサ固有モデルに切り替える切替ステップと、
を含む請求項８または９に記載の生体パラメータ推定方法。
【請求項１１】
前記センサ群に含まれる各センサからの信号を前処理し、前処理済み信号が、前記第１信号として、前記初期化ステップおよび前記推定ステップに出力される前処理ステップをさらに含む請求項５に記載の生体パラメータ推定方法。
【請求項１２】
前記前処理ステップは、
前記第１信号に帯域フィルタ処理を施して第１帯域フィルタ処理済み信号とする第１帯域フィルタ処理ステップと、
前記第１帯域フィルタ処理済み信号を第１規格化信号に規格化する第１規格化ステップと、
前記第１規格化信号に非線形変換を施して変換済信号とする変換ステップと、
前記変換済信号に帯域フィルタ処理を施して第２帯域フィルタ処理済み信号に変換する第２帯域フィルタ処理ステップと、
前記第２帯域フィルタ処理済み信号を中心化して中心化信号を生成する中心化ステップと、
前記中心化信号を規格化して前記前処理済み信号とする第２規格化ステップと、
を含む請求項１１に記載の生体パラメータ推定方法。
【請求項１３】
処理装置用のコンピュータプログラムであって、
前記処理装置に、請求項１から１２までのいずれかの生体パラメータ推定方法を実行させるコンピュータプログラム。
【請求項１４】
請求項１３に記載のコンピュータプログラムを格納したコンピュータ可読媒体。
【請求項１５】
前記処理装置の内部メモリに直接ロード可能な請求項１３に記載のコンピュータプログラム。
【請求項１６】
圧力変化を測定する２個以上のセンサ群を備え加速度計に接続された支持部材上に存在する対象の、生体パラメータを推定する生体パラメータ推定装置であって、
前記対象の生体パラメータ推定に用いる１つの選択センサを、前記加速度計からの第２信号に基づいて、タイムフレームごとに、前記センサ群から選択する選択手段と
前記センサ群により測定された圧力変化に対応する第１信号に基づいて前記センサ群に含まれるセンサごとに提供されたセンサ固有モデルのうち、前記選択センサに対応するセンサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定手段と、
を備え、
前記選択手段は、非線形モデルを備え、
帰無仮説パラメータ、確率密度関数パラメータ、カルバック・ライブラー情報量パラメータ、歪度パラメータおよび尖度パラメータの少なくともいずれか１つを含む統計パラメータおよびｑ−ハーストパラメータの少なくともいずれか一方を含むパラメータであって、前記第１信号および前記第２信号のいずれか一方からタイムフレームごとに抽出されたパラメータを備えた特徴ベクトルを算出し、
前記特徴ベクトルを前記非線形モデルに入力して、前記特徴ベクトルのパラメータを、タイムフレームごとに、前記選択センサの番号に割り当てる生体パラメータ推定装置。
【請求項１７】
前記非線形モデルは、非線形関数と確率に基づいた決定とを組み合わせることにより、前記特徴ベクトルのパラメータを、タイムフレームごとに、前記選択センサの番号に割り当てる請求項１６に記載の生体パラメータ推定装置。
【請求項１８】
圧力変化を測定する２個以上のセンサ群を備え加速度計に接続された支持部材上に存在する対象の、生体パラメータを推定する生体パラメータ推定装置であって、
前記対象の生体パラメータ推定に用いる１つの選択センサを、前記加速度計からの第２信号に基づいて、タイムフレームごとに、前記センサ群から選択する選択手段と
前記センサ群により測定された圧力変化に対応する第１信号に基づいて前記センサ群に含まれるセンサごとに提供されたセンサ固有モデルのうち、前記選択センサに対応するセンサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定手段と、
を備え、
前記推定手段は、前記センサ群に含まれる各センサに対し、前記第１信号のそれぞれに基づいて、状態ベクトルの、初期周波数を含む状態パラメータを推定する初期化処理を行ない、
前記推定手段は、前記初期化処理において、前記センサ群のそれぞれに対して、前記状態パラメータごとに、共分散値をプロセスノイズ値にマッピングする非線形マッピングモデルを利用してプロセスノイズ値を取得し、取得したプロセスノイズ値を用いて前記対象の生体パラメータを推定し、
前記推定手段は、前記非線形マッピングモデルを含み、
前記非線形マッピングモデルは、
非線形マッピングモデルの中間層のバイアスベクトルと前記共分散値より計算した変換値を乗じた前記中間層の係数ベクトルとのシグモイド関数を計算する第１計算ステップと、
前記非線形マッピングモデルの出力層のバイアスベクトルと、前記第１計算ステップの結果を乗じた前記出力層の係数ベクトルとのマッピング関数を計算する第２計算ステップと、
を実行することにより、共分散値をプロセスノイズ値に割り当てる生体パラメータ推定装置。
【請求項１９】
前記選択手段は、
前記第１信号および前記第２信号に基づいて、前記選択センサを選択する請求項１６または１７に記載の生体パラメータ推定装置。
【請求項２０】
前記推定手段は、前記センサ群に含まれる各センサに対し、前記第１信号のそれぞれに基づいて、状態ベクトルの、初期周波数を含む状態パラメータを推定する初期化処理を行なう請求項１６または１７に記載の生体パラメータ推定装置。
【請求項２１】
前記推定手段は、前記初期化処理において、前記センサ群のそれぞれに対して、前記状態パラメータごとに、共分散値をプロセスノイズ値にマッピングする非線形マッピングモデルを利用してプロセスノイズ値を取得し、取得したプロセスノイズ値を用いて前記対象の生体パラメータを推定する請求項２０に記載の生体パラメータ推定装置。
【請求項２２】
前記推定手段は、前記非線形マッピングモデルを含み、
前記非線形マッピングモデルは、
非線形マッピングモデルの中間層のバイアスベクトルと前記共分散値より計算した変換値を乗じた前記中間層の係数ベクトルとのシグモイド関数を計算する第１計算ステップと、
前記非線形マッピングモデルの出力層のバイアスベクトルと、前記第１計算ステップの結果を乗じた前記出力層の係数ベクトルとのマッピング関数を計算する第２計算ステップと、
を実行することにより、共分散値をプロセスノイズ値に割り当てる請求項２１に記載の生体パラメータ推定装置。
【請求項２３】
前記推定手段は、前記選択手段において計算されたセンサごとのセンサ固有モデルを用いて、前記センサ群に含まれる各センサについて、前記第１信号の前サンプル（ｋ−１）の状態ベクトルに基づいて、前記第１信号の現サンプル（ｋ）の状態ベクトルを算出する請求項２０、２１または２２に記載の生体パラメータ推定装置。
【請求項２４】
前記推定手段は、前記センサ群の状態パラメータごとに得られたプロセスノイズ値を用いて、前記状態ベクトルを算出する請求項２３に記載の生体パラメータ推定装置。
【請求項２５】
前記推定手段は、前記状態ベクトルを算出するため、
前記前サンプルの前状態ベクトルとセンサごとのセンサ固有モデルの線形モデルとに基づいて前記現サンプルの現状態ベクトルを予測し、
予測された前記現状態ベクトルを、前記センサごとのセンサ固有モデルの非線形モデルを用いて更新し、
前記選択手段において計算された１個のセンサに基づいて、生体パラメータ推定のために用いるべきセンサ固有モデルを、前記選択センサに提供されるセンサ固有モデルに切り替える請求項２３または２４に記載の生体パラメータ推定装置。
【請求項２６】
前記推定手段は、前記センサ群に含まれる各センサからの信号を前処理し、前処理された信号を、前記第１信号として用いる請求項２０に記載の生体パラメータ推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、支持部材の上に置かれた対象の生体パラメータを推定する技術に関するものである。言い換えれば、本発明は人体の生理的パラメータ推定技術、特に、例えば車両の運転手および乗客などの脈拍や呼吸数のような生理的パラメータを推定する技術に関するものである。
【背景技術】
【０００２】
生体パラメータのモニタリングに関しては、特許文献１乃至３に記載がある。例えば脈拍信号および呼吸信号は、圧電センサを利用してモニタリング可能である。圧電センサは、圧力変化を測定する。例えば家庭や乗用車の座席に取り付けられた圧電センサにより、血流圧が引き起こす圧力変化を測定することができる。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】ＦＲ２９４３２３３Ａ１
【特許文献２】ＦＲ２９４３２３４Ａ１
【特許文献３】ＦＲ２９４３２３６Ａ１
【発明の概要】
【発明が解決しようとする課題】
【０００４】
本発明は、特許文献１乃至３に記載されている生体パラメータモニタリングを改善し、例えば座席やベッドなどの支持部材上に、圧力変化測定用の２個以上のセンサを埋め込むことにより生体パラメータをロバストに推定する方法を提供することを目的としている。
【課題を解決するための手段】
【０００５】
上記目的を達成するため、本発明に係る方法は、
圧力変化を測定する２個以上のセンサ群を備え、１個以上の加速度計に接続された支持部材上に存在する対象の生体パラメータを推定するための推定方法であって、
前記センサ群により測定された圧力変化に対応する第１信号に基づいて、前記センサ群に含まれるセンサそれぞれにセンサ固有モデルを提供するステップと、
前記加速度計からの第２信号に基づいて、前記対象の生体パラメータ推定に使用する１つの選択センサを、タイムフレームごとに前記センサ群から選択する選択ステップと、
前記選択センサに対応する前記センサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定ステップと、
を含み、
前記選択ステップにおいて、
帰無仮説パラメータ、確率密度関数パラメータ、カルバック・ライブラー情報量パラメータ、歪度パラメータおよび尖度パラメータの少なくともいずれか１つを含む統計パラメータおよびｑ−ハーストパラメータの少なくともいずれか一方を含むパラメータであって、前記第１信号および前記第２信号のいずれか一方からタイムフレームごとに抽出されたパラメータを備えた特徴ベクトルを算出し、
前記特徴ベクトルを非線形モデルに入力して、前記特徴ベクトルのパラメータを、タイムフレームごとに、前記選択センサの番号に割り当てる。
上記目的を達成するため、本発明に係る他の方法は、
圧力変化を測定する２個以上のセンサ群を備え、１個以上の加速度計に接続された支持部材上に存在する対象の生体パラメータを推定するための推定方法であって、
前記センサ群により測定された圧力変化に対応する第１信号に基づいて、前記センサ群に含まれるセンサそれぞれにセンサ固有モデルを提供するステップと、
前記加速度計からの第２信号に基づいて、前記対象の生体パラメータ推定に使用する１つの選択センサを、タイムフレームごとに前記センサ群から選択する選択ステップと、
前記選択センサに対応する前記センサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定ステップと、
を含み、
前記推定ステップは、前記センサ群のいずれに対しても、前記第１信号のそれぞれに基づいて、状態ベクトルの、初期周波数を含む状態パラメータを推定する初期化ステップを含み、
前記初期化ステップにおいて、前記センサ群のそれぞれに対して、前記状態パラメータごとに、共分散値をプロセスノイズ値にマッピングする非線形マッピングモデルを利用して、プロセスノイズ値を取得し、取得したプロセスノイズ値を前記推定ステップに利用し、
非線形マッピングモデルの中間層のバイアスベクトルと前記共分散値より計算した変換値を乗じた前記中間層の係数ベクトルとのシグモイド関数を計算する第１計算ステップと、
前記非線形マッピングモデルの出力層のバイアスベクトルと、前記第１計算ステップの結果を乗じた前記出力層の係数ベクトルとのマッピング関数を計算する第２計算ステップと、
を実行することにより、前記マッピングモデルが、共分散値をプロセスノイズ値に割り当てる。
【０００６】
上記目的を達成するため、本発明に係る装置は、
圧力変化を測定する２個以上のセンサ群を備え加速度計に接続された支持部材上に存在する対象の、生体パラメータを推定する生体パラメータ推定装置であって、
前記対象の生体パラメータ推定に用いる１つの選択センサを、前記加速度計からの第２信号に基づいて、タイムフレームごとに、前記センサ群から選択する選択手段と
前記センサ群により測定された圧力変化に対応する第１信号に基づいて前記センサ群に含まれるセンサごとに提供されたセンサ固有モデルのうち、前記選択センサに対応するセンサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定手段と、
を備え、
前記選択手段は、非線形モデルを備え、
帰無仮説パラメータ、確率密度関数パラメータ、カルバック・ライブラー情報量パラメータ、歪度パラメータおよび尖度パラメータの少なくともいずれか１つを含む統計パラメータおよびｑ−ハーストパラメータの少なくともいずれか一方を含むパラメータであって、前記第１信号および前記第２信号のいずれか一方からタイムフレームごとに抽出されたパラメータを備えた特徴ベクトルを算出し、
前記特徴ベクトルを前記非線形モデルに入力して、前記特徴ベクトルのパラメータを、タイムフレームごとに、前記選択センサの番号に割り当てる。
上記目的を達成するため、本発明に係る他の装置は、
圧力変化を測定する２個以上のセンサ群を備え加速度計に接続された支持部材上に存在する対象の、生体パラメータを推定する生体パラメータ推定装置であって、
前記対象の生体パラメータ推定に用いる１つの選択センサを、前記加速度計からの第２信号に基づいて、タイムフレームごとに、前記センサ群から選択する選択手段と
前記センサ群により測定された圧力変化に対応する第１信号に基づいて前記センサ群に含まれるセンサごとに提供されたセンサ固有モデルのうち、前記選択センサに対応するセンサ固有モデルを用いて、タイムフレームごとに、前記対象の生体パラメータを推定する推定手段と、
を備え、
前記推定手段は、前記センサ群に含まれる各センサに対し、前記第１信号のそれぞれに基づいて、状態ベクトルの、初期周波数を含む状態パラメータを推定する初期化処理を行ない、
前記推定手段は、前記初期化処理において、前記センサ群のそれぞれに対して、前記状態パラメータごとに、共分散値をプロセスノイズ値にマッピングする非線形マッピングモデルを利用してプロセスノイズ値を取得し、取得したプロセスノイズ値を用いて前記対象の生体パラメータを推定し、
前記推定手段は、前記非線形マッピングモデルを含み、
前記非線形マッピングモデルは、
非線形マッピングモデルの中間層のバイアスベクトルと前記共分散値より計算した変換値を乗じた前記中間層の係数ベクトルとのシグモイド関数を計算する第１計算ステップと、
前記非線形マッピングモデルの出力層のバイアスベクトルと、前記第１計算ステップの結果を乗じた前記出力層の係数ベクトルとのマッピング関数を計算する第２計算ステップと、
を実行することにより、共分散値をプロセスノイズ値に割り当てる生体パラメータ推定装置。
【発明の効果】
【０００７】
本発明によれば、必要に応じて推定に用いるセンサを変えることにより、ロバストに生体パラメータを推定できる。
【図面の簡単な説明】
【０００８】
【図１】本発明の実施形態に係る装置の構成を示す概略ブロック図である。
【図２】本発明の実施形態に係るセンサの自動的変更原理を説明するための概略ブロック図である。
【図３】相異なる２つの運転ケースに対して、正規分布からの変動を示す図である。
【図４】１４個のセンサを使用した本発明の実施形態の場合の、全てのセンサの推定確率密度関数を示す図である。
【図５】運転中に、センサを特定せずに全てのセンサに対して２０秒間の脈拍推定を実施した結果を示す図である。
【図６】本発明の実施形態に係る、ｑ−ハーストパラメータをセンサ番号に割り付けるモデルを説明するための概略ブロック部である。
【図７】本発明の実施形態に係る、センサ前処理工程を示す概略ブロック図である。
【図８】前処理に使用されるバンドパスフィルタの振幅応答および位相応答を示す図である。
【図９】前処理によって非線形的に変換されるセンサ信号を示す図である。
【図１０】本発明の実施形態に係る初期周波数推定原理を説明するための概略ブロック図である。
【図１１】ＥＳＰＲＩＴ周波数推定原理を示す図である。
【図１２】本発明の実施形態に係る初期周波数決定用のクラスタリング原理を示す図である。
【図１３】ニューラルネットワークを用いた非線形フィッティング方法を示す概略ブロック図である。
【図１４】ＩＭＭ−ＥＫＦ処理を説明するための図である。
【図１５】本発明の実施形態に係る自動的センサ変更を利用した脈拍推定結果を示す図である。
【発明を実施するための形態】
【０００９】
本発明によれば、支持部材上に存在する対象の、心拍数や呼吸リズムのような生体パラメータを推定する。支持部材は、例えば家にある椅子や車両の座席、またはベッドなどである。支持部材は、例えば圧電センサなどの圧力変化測定用センサを２個以上備えている。
【００１０】
支持部材はさらに少なくとも１つの加速度計を備えている。支持部材が車両座席の場合には加速度計は座席台に取り付けられ、特に車両からの振動ノイズなどの周辺ノイズを検出する参照センサとして、３軸方向のノイズを検出する。
【００１１】
運転中、ドライバーの身体は頻繁に動き、その動きは時に大きな振幅を発生させる。ここで考慮されている身体の動きとは、例えば個別または同時に動く、脚、手、身体全体、頭などの、全ての部分の動きである。身体の動きは、運転に関係する動き、および、生理学的または心理学的状態に関係する動き、という２つのカテゴリーに分類できる。
【００１２】
ここでは、例えば圧電センサなどの少なくとも２つのセンサが、座席などの支持部材に取り付けられており、ノイズは無いかまたは既に除去されているものと仮定する。
【００１３】
（装置構成）
図１は本実施形態に係る装置の構成を示す概略ブロック図である。本装置は、インタラクティブなマルチモデル拡張カルマンフィルタによりセンサを自動的に変更する。
【００１４】
本実施形態で説明する信号と関数は、ディジタル化されている。
【００１５】
図１に示す通り、本装置は、指定された時刻における脈拍を測定するのに使用するのに最も適していると思われるセンサを予測し選択するための自動センサ変化推定部（ＡＳＣＥブロック）３０と、新しいセンサが選択されたときに線形状態空間モデル２１と非線形状態空間モデル２２とをを切り替える、インタラクティブなマルチモデル拡張カルマンフィルタ（ＩＭＭ−ＥＫＦブロック）２０とを含む。
【００１６】
例えば車両の座席などの支持部材（図には示されていない）のセンサ群（例えば座席上部の圧電センサおよび座席底部の圧電センサ）からの信号６１が、装置１によって処理される。信号６１は、ノイズが除去され、ディジタル化されたセンサ出力である、例えば圧電センサからの出力である。装置１は、例えば車両座席などの支持部材に取り付けられた加速度計からの信号１６１をも処理する。信号１６１は加速度計からディジタル化して出力される。
【００１７】
本発明によれば、ＩＭＭ−ＥＫＦブロック２０では、各タイムフレームで１個のセンサだけを使用する。したがって、Ｎ個のセンサがあると仮定すると、ＡＳＣＥブロック３０は、所定のタイムフレームＴごとに、Ｎ個のセンサから１つのセンサを選択する。タイムフレームＴは、例えば５００ｍｓ、１ｓ、または２ｓのオーダーである。
【００１８】
装置１は本発明の一実施形態を示し、圧力変化を測定するＮ個（Ｎ≧２）のセンサを備えた支持部材上の対象の生体パラメータ４０を推定する。支持部材にはＭ個（Ｍ≧１）の加速度計が接続される。装置１には、Ｎ個のセンサの一つ一つについて、Ｎ個のセンサからの信号６１に基づく、線形／非線形状態空間モデル（センサ固有モデル）２１、２２が含まれる。
【００１９】
信号６１は、Ｎ個のセンサのそれぞれにより測定された圧力変化に対応している。選択処理においては、各タイムフレームに対してＡＳＣＥブロック３０がＮ個のセンサから１個のセンサを選択し、Ｍ個の加速度計からの信号１６１、または信号１６１と信号６１の両方に基づき、ＩＭＭ−ＥＫＦブロック２０において支持部材の生体パラメータを推定するために使用する。推定処理においては、ＩＭＭ−ＥＫＦブロック２０が、タイムフレームＴごとに、生体パラメータ推定および追跡ブロック２４において、選択した１個のセンサに設定されモデル切替ブロック２３により切り替えられたセンサ固有モデル２１、２２を利用して対象の生体パラメータ４０を推定する。
【００２０】
ＡＳＣＥブロック３０におけるブロックＨ0は、ノイズが正規分布（ガウス分布）を示す場合にはテストを含む。ＡＳＣＥブロック３０におけるブロックＰＤＦは、ランダム変数（ここではノイズ）の相対確率を表現する確率密度関数で構成される。ＡＳＣＥブロック３０におけるブロックＫＬはカルバック・ライブラー情報量により構成され、さらにＡＳＣＥブロック３０におけるｑ−ハーストブロックは、ランダムマルチフラクタル理論によるｑ−ハーストパラメータにより構成される。
【００２１】
１．センサ自動変更
本実施形態によるセンサ自動変更の原理を図２に示す。
【００２２】
タイムフレームＴごとに、各加速度計信号１６１、または加速度計信号１６１およびセンサ信号６１に対して特徴パラメータを計算する。加速度計は振動ノイズをも測定することになり、なんらかの身体の動きにも影響を受ける。
【００２３】
ＡＳＣＥブロック３０のセンサ自動変更モデル３１には特徴ベクトルが含まれており、最終的には、ＩＭＭ−ＥＫＦブロック２０での推定処理において使用するセンサ番号（センサ＃）を決定する。
【００２４】
１．１  特徴ベクトル
特徴ベクトルとは、加速度計およびセンサから抽出されたパラメータの集合のことである。パラメータの種類は制限されていない。例えば以下のパラメータが使用可能である。
【００２５】
ａ）ｑ−ハーストパラメータ
これらのパラメータは、ランダムマルチフラクタル理論に基づいている。背後にある考え方は、ノイズおよびセンサ信号に対して不変である特定の構造を見出そうとするものである。ｑ−ハースト指数は、加速度計信号１６１およびセンサ信号６１のマルチフラクタル構造をパラメータ化するために利用される。ｑ−ハーストパラメータは、マルチフラクタルのトレンド除去ゆらぎ解析を利用して計算する。
【００２６】
最初に下記の累積和によって信号を積分する。
【数１】
0006738177.tif 000002
ここでｘintは加速度計信号またはセンサ信号の積分値、ｘは加速度計信号またはセンサ信号を示す。
【００２７】
平均振幅変動、すなわち二乗平均平方根（ＲＭＳ）は下記の式で計算される。
【数２】
0006738177.tif 000003
ここでＦはスケールｓおよびセグメントインデックスｖのＲＭＳ値である。ｘfitはｘintの２次トレンド除去を示し、下記の式で計算される。
【数３】
0006738177.tif 000004
ここで、ａはフィッティング係数である。
【００２８】
ｑパラメータにより一般化すると、マルチフラクタル振幅変動に対して下記の式が得られる。
【数４】
0006738177.tif 000005
したがって、ｌｏｇ2Ｆq（ｓ）の傾きを評価することにより、ｑハーストパラメータが得られる。
【００２９】
ハーストパラメータにおいて「ｑ」は追加の尺度であり、ｑが負の値の場合には非常に小さい変動を有するセグメントが増幅され、ｑが正の値の場合には非常に大きい変動を有するセグメントが増幅される。
【００３０】
したがって、いくつかのスケールのｓ、例えばｓ＝４，８，１６，３２，６４，１２８，２５６，５１２の各値に対して計算を実行する。
【００３１】
上記の計算を、所用の各ｑの値、例えばｑ＝−５，−１，５の値に対して繰り返す。
【００３２】
ｑ−ハーストパラメータは勾配であり、もしあるｑ値に対して３個の値が対応する場合、１つの特徴ベクトルに対し、センサの数に応じてＮ＊３の長さを有する。例えば座席に１４個のセンサーがある場合にセンサに対して上記計算を実施すると、センサ自動変更モデル３１の入力となる特徴ベクトルの長さは４２となる。
【００３３】
上記の特徴ベクトルを利用することにより、最良のセンサ、すなわち生体パラメータ推定に利用できるセンサを決定することが可能である。
【００３４】
ｂ）統計パラメータ
もう一つの可能性は、統計的記述子を組み合わせて特徴ベクトルを作成することである。かかる計算は各タイムフレームＴに対しても実行される。
【００３５】
分布が正規分布に従うという帰無仮説（Ｈ0）検定
この計算は、全てのセンサまたは加速度計に対して実行され、正規性のリリーフォース検定から構成される。
【００３６】
データ分布が正規分布であればＨ0＝０であり、さもなければＨ0＝１である。
【００３７】
最初のステップは、センサ信号６１と加速度計信号１６１との両方またはいずれかの平均値と標準偏差値の計算であり、次の式で与えられる。
【数５】
0006738177.tif 000006
ここでｘkは現在のタイムフレームのｋ番目のデータサンプルである。このタイムフレームにはＬ個のサンプルが存在する。ｘバーおよびσは、それぞれサンプルの平均値およびサンプルの標準偏差である。
【００３８】
第２ステップでは、かかるタイムフレームの全てのサンプルに対して規格化した値を計算する。
【数６】
0006738177.tif 000007
従ってリリーフォース検定、ＬＦは次式で計算できる。
【数７】
0006738177.tif 000008
ここでＬＦは、平均が０で分散が１の正規分布（Ｆ（ｘ））と、Ｚk値実験分布との差の絶対値上限である。
【００３９】
ＬＦが検定臨界値より大きい場合には、検定は実施しない（臨界値は表によって与えられる）。
【００４０】
したがって、正規確率分布と比較した場合のノイズ確率分布の変動性の尺度が導出される。各タイムフレームに対して前記比較が実施される。
【００４１】
図３は、センサのノイズ確率密度が正規分布である場合とない場合の差異を示す。ｘ軸はフレーム番号としての時刻であり、ｙ軸はセンサ番号である。白色で正規分布の場合を、黒色で正規分布でない場合を示す。
【００４２】
確率密度関数（ＰＤＦ）
上記はノイズ分布がガウス分布にならない場合も多いことを示しており、従って全ての加速度計やセンサに対する分布の形状を詳細に示す説明を追加するのも興味深い。
【００４３】
ＰＤＦは、加速度計信号またはセンサ信号のデータのヒストグラムを取ることによって計算することも可能であるが、真のＰＤＦに収束するカーネル密度推定方法を用いるのが好ましく、次式により与えられる。
【数８】
0006738177.tif 000009
ここで、ｈは帯域幅、Ｋはカーネルである。例えばカーネルとしてはガウス関数、一様、Epanechnikovなどが考えられる。
【００４４】
図４は、所定の時刻でのノイズ確率密度関数を推定したものを示す。１４個全てのセンサは水平軸上におかれ、各センサの対応範囲はグラフ上の矢印で示している。
【００４５】
図４からわかるように、各センサにより分布は異なっている。いくつかはフラット（例えばセンサ＃１２）であり、いくつかは大きなテールをひき（例えばセンサ＃７）、またいくつかは非対称である（例えばセンサ＃１、＃４、＃８）。これらの分布は各タイムステップによっても変化し、従ってセンサまたは加速度計に何が発生しているかに関して重要な情報を与えてくれる。
【００４６】
カルバック・ライブラー（ＫＬ）情報量交差ＰＤＦ
各時刻におけるＰＤＦ推定値が得られている。他のセンサと比較した場合の特定のセンサ分布の変化の尺度を知っておくことも重要である。
【００４７】
したがって所定の時刻ｔi-1においては全てのセンサに対して、時刻ｔiにおいては交差する全てのセンサに対して、ＫＬ情報量を計算する。
【００４８】
ＫＬ情報量は次式で表される。
【数９】
0006738177.tif 000010
ここでｆおよびｈは、情報量を計算するための２つのＰＤＦである。
【００４９】
従って１４個のセンサがある場合には、所定の時刻において１４×１４＝１９６個の値を計算する。
【００５０】
追加の統計パラメータ
さらに３次および４次の統計モーメント、すなわち歪度（skewness）と尖度（kurtosis）がそれぞれ計算可能であり、特徴ベクトルに追加される。これらは次式で与えられる。
【数１０】
0006738177.tif 000011
【数１１】
0006738177.tif 000012
特徴ベクトルは、これら全てのパラメータにより構成される。
【００５１】
１．２  ターゲット
ターゲットは、センサ数である。各タイムフレームに対して、例えば最初に全てのセンサ上での脈拍推測地をチェックするなどして、最良と考えられるセンサを決定しなければならない。
【００５２】
図５は、所定の対象とランに対して最良のセンサを見つける方法の一例を示す。拡張カルマンフィルタを用いて２０秒間脈拍を推定したものを、全てのセンサに対して拡大して示す。時間は１秒区切りに分割されており（ここでは最初の１０秒だけを示す）、各区切りにおいて最良のセンサが示されている（＃５はセンサ番号５を示す、等々）。図５には、レファランスとしてＥＣＧ信号も示している。
【００５３】
従って図２に示すセンサ自動変更モデル３１は、入力の特徴ベクトルと出力のターゲットとの間で非線形なマッピングを実行することになる。
【００５４】
１．３  モデル
非線形モデルはニューラルネットワーク（ＮＮ）に基づいて使用されている。ＮＮによりあらゆるタイプの非線形性がモデル化可能である。
【００５５】
分類を実行する（すなわち最良だと考えられるセンサを決定する）ために用いるニューラルネットワークの構造を図６に示す。
【００５６】
図６には、入力特徴ベクトルがセンサごとに３つのｑ−ハースト・パラメータに基づいている場合の、分類器（すなわち決定）モデル３１の構造の一例を示す。
【００５７】
入力特徴ベクトルは３つのｑ−ハースト・パラメータに限定されるものではなく、３つ以上のパラメータを使用してもよい。
【００５８】
図６に示すように、各センサあたり３つのｑ−ハースト・パラメータに対応して、入力では４２個のパラメータが存在し、出力では１４個のパラメータが存在する。出力の各要素はセンサ番号に対応している。すなわち、要素１はセンサ＃１に、要素２はセンサ＃２に対応している。
【００５９】
センサ番号は、最大値を持つ、すなわち最も高い確率を持つ出力要素に決定される。全てのセンサ番号に対して秒単位で確率が得られるため、かかる構造が他の種類の決定時にも利用される。
【００６０】
かかる構造は、勿論入力特徴ベクトルの種類にも適用する必要がある。
【００６１】
かかる構造は、中間層（隠れ層hidden layer）および出力層の存在を示している。中間層のニューロン数は、構造と入力数（特徴ベクトル）に依存する。ニューロン数は例えば２５０個というような数であるが、信頼できるセンサ番号を予測するために必要であればいくつであってもよい。中間層には重量（ｗ）、バイアス（ｂ）、および非線形関数が含まれ、この場合（ニューロンが１個）にはシグモイドである。出力層には重量（ｗ）、バイアス（ｂ）、および非線形関数が含まれ、ここではソフトマックス関数である。
【００６２】
使用されるシグモイドおよびソフトマックス関数は次式で表される。
【数１２】
0006738177.tif 000013
ソフトマックス関数により、確実性の尺度（すなわち事後確率）が得られる。
【００６３】
最初にパラメータ計算用にデータを準備する。このことを以下に説明する。
【数１３】
0006738177.tif 000014
入力行列には、１４個全てのセンサのｑ−ハーストパラメータに対応して４２行、さらに１個のパラメータあたり与えられる秒数データに対応してＣ個の列が必要である。前述のように、新規の列はそれぞれ相異なる時刻に対応する。ターゲット行列は入力行列と同じ数の列を有するが、行数はそれぞれ各センサに対応して１４行のみである。従って選択された最良センサに対応したものだけが１に設定され、他の値は０でなければならない。
【００６４】
以上よりモデル３１のパラメータが計算できる。この処理はトレーニングとも呼ばれるが、反復計算を利用し、さらに同じターゲット値に対して相異なる入力値を持ち得るからである。この計算は、スケール共役勾配法のバックプロゲーションアプローチを用いて実行される。
【００６５】
統計的ディスクリプタ特徴ベクトルなどの異なる特徴ベクトルに対しても、同じ計算方法が利用可能である。
【００６６】
非線形モデル３１のパラメータが一度計算されれば、脈拍推定に利用されるセンサを決定する際に利用できる。
【００６７】
引き続き各タイムフレームにおいて、センサまたは加速度計用にｑ−ハーストパラメータまたは統計的ディスクリプタが計算され、パラメータの特徴ベクトル、Ｘrearureが作成される。
【００６８】
前記ベクトルは以下の計算に使用される。まず入力を次式を使用して変換する。
【数１４】
0006738177.tif 000015
ここでｘoffsetは特徴ベクトルから取り除くべきオフセット、Ｇはゲインを示す。これらのオフセットとゲインは、前に示したオフライン手順によって計算済である。
【００６９】
次に中間層の出力を次式にて計算する。
【数１５】
0006738177.tif 000016
ここでＢhiddenは中間層のバイアスベクトルであり、中間層のニューロン数、Ｑに等しい長さを有する。
【００７０】
Ｗhiddenは中間層の係数行列であり、「Ｑ×Ｐ」の大きさを有する。Ｐは特徴ベクトルの長さである。ｙsigmoidは非線形のシグモイド関数であり、その式は既に示した。中間層出力、ｙhiddenは長さＱのベクトルである。
【００７１】
最後のステップは出力層の出力を計算することであり、次式で与えられる。
【数１６】
0006738177.tif 000017
ここでＢoutは出力層のバイアスベクトルであり、センサ数Ｎと等しい長さを有する。
【００７２】
Ｗoutは出力層の係数行列であり、「Ｎ×Ｑ」の大きさを有する。ｙsoftmaxは非線形なソフトマックス関数であり、その式はすでに示した。出力ｙoutは、センサ数Ｎと等しい長さを有するベクトルである。
【００７３】
以上より最良と思われるセンサは、ｙoutが最大となる要素番号のものに対応する。
【００７４】
ニューラルネットワークは必ずしも唯一の可能性を示すものではなく、隠れマルコフモデルを利用してもよい。
【００７５】
センサ数はＩＭＭ−ＥＫＦブロック２０に伝えられる。
【００７６】
２．ＩＭＭ−ＥＫＦによる頑健な脈拍推定
特許文献１−３には、拡張カルマンフィルタの説明がある。
【００７７】
本アプリケーションにおいては、脈拍推定の際にセンサが１個しか保持されてなく、しかも時間とともに変化するので、センサ自動変更に対して前述のアプローチを取る必要がある場合には、スイッチ観測モデルと状態空間モデルとの切り替え用に、必要な場合にはカルマンフィルタを適用する。
【００７８】
ＩＭＭ−ＥＫＦブロック２０の役割は、振動ノイズにもよらず、たとえ身体がセンサと接触していても身体の動きにもよらず、脈拍を頑健に推定することである。ＩＭＭ−ＥＫＦブロック２０は、モデル切替が可能なＥＫＦ拡張版である。基本的にはミキシング、フィルタリング、および組み合わせという３つの主要ステップを含む。
【００７９】
モデル変更は主として確率のミキシングに基づいて決定される。全てのモデルに対してＥＫＦ評価を実行する。
【００８０】
本発明によれば、推定に利用するセンサ番号の情報がわかっており、オリジナルのＩＭＭ−ＫＴＦアプローチは、この情報が利用できるように修正されている。Ｎ個のセンサ間で切り替える可能性があるので、Ｎ個の状態空間モデルとＮ個の観測モデル（線形状態空間モデル２１および非線形状態空間モデル２２）とが提供されている。
【００８１】
２．１  圧電センサの前処理
低減されたとはいえまだノイズが残っているＩＭＭ−ＥＫＦブロック２０のセンサ信号６１を利用する前に、前処理を実行して推定性能を最大化する。図７はセンサ信号６１に対する前処理を示す（例えば座席底部の圧電センサおよび座席上部の圧電センサとから構成される圧電信号）。
【００８２】
ａ）  バンドパスフィルタ
まずセンサ信号６１に帯域フィルタ処理を施す。この目的のために使用されるバンドパスフィルタは、脈拍推定に関係する周波数に中心がおかれている。無限インパルス応答（ＩＩＲ）フィルタは、次の特性を有するように設計されている。
〇周波数中心値：  ｆ0＝１．３Ｈｚ
〇通過帯域：  ｂ＝２．５Ｈｚ
〇オーダー：  ３
このようなフィルタをデザインするために、次の計算を実行する。
【数１７】
0006738177.tif 000018
さらに次のベクトルを定義する。
【数１８】
0006738177.tif 000019
フィルタのオーダーが３であるので、これらの値を次のように２回畳み込む。
【数１９】
0006738177.tif 000020
Ｂ＝Ｂ0；Ａ＝Ａ0
Ｂ＝Ｂ＊Ｂ0（２回実施）
Ａ＝Ａ＊Ａ0（２回実施）
ここで、「＊」は畳み込みを示す（乗算ではない）。
【００８３】
ＩＩＲフィルタのデザインが終了したたので、双一次変換を利用してデザインしたものをＺドメインに移す。
【００８４】
双一次変換を実行するために、次式で示すワープ周波数を使用する。
【数２０】
0006738177.tif 000021
フィルタの特性を図８に示す。図８において、実線がバンドパスフィルタの振幅応答を、一点鎖線が位相応答を示す。
【００８５】
ｂ）  規格化
引き続きバンドパスフィルタの信号を規格化する。最初の規格化において、センサ信号を、最初のフレーム（すなわち５秒間）の標準偏差で割り算する。
【数２１】
0006738177.tif 000022
ここでＳn,iはセンサｉの規格化した信号、Ｓiはセンサ信号、Ｓi|o<t<Tframeはセンサiのフレームサンプル、stdは標準偏差σである。
【００８６】
ｃ）  非線形変換
振動や身体の動きにより、状況に応じてセンサ信号に大きな振幅変動がもたらされる。この大きな振幅変動が、脈拍推定に大きな影響を及ぼす可能性がある。そこで非線形変換のステップは、一定の振幅値ではあるが振動が同じ状態で維持されている場合の信号を与えることを目的としている。
【００８７】
非線型変換関数として双曲線正接を利用し、次式にて非線形変換後の信号を計算する。
【数２２】
0006738177.tif 000023
図９は非線形変換結果（実線）を変換前のセンサ信号（一点鎖線）と比較したものを示す。同じ図面上で変換後の信号と比較できるように、一点鎖線は振幅で規格化している。実際のセンサ信号値は１０，０００倍である。
【００８８】
図９から分かるように、入力にいかなる変動があっても振幅変動は一定であり、振動周期も同じ状態に維持されている。
【００８９】
ｄ）  状態パラメータの初期化
非線形変換の後、ａ）項で説明したものと同じバンドパスフィルタを適用する。これにより、信号が正弦波により近くなる。
【００９０】
ｅ）  中心化
引き続き再度帯域フィルタ処理を実行し、以下に示す中心化処理を実施することにより、信号からサンプルの平均値を取り除くことができる。
【数２３】
0006738177.tif 000024
ここでＹc,iは中心化後の信号、Ｙbp,iはｄ）項で説明した帯域フィルタ処理後の信号、Ｙbp,iバーはサンプルの平均値である。
【００９１】
ｆ）  規格化
この最後のステップでは、前処理済みのセンサ信号に対して最終的な規格化を実施する。現在の時間フレームＴでの信号継続時間を標準偏差により規格化する。
【数２４】
0006738177.tif 000025
【００９２】
２．２  状態パラメータの初期化
推定手順を容易にするため、ＩＭＭ−ＥＫＦブロック２０の状態パラメータはベクトル形式に設定されている。Ｎ個のセンサがあるのでＮ個のモデル２１、２２が存在し、さらに推定に利用するためＮ個の状態ベクトルも存在する。前述のようにＮは例えば１４に等しくてもよいし、これに限定される必要もない。
【００９３】
推定すべきパラメータは、周波数ｆ、振幅Ａ、および位相φである。従って状態ベクトルは次式のようになる。
【数２５】
0006738177.tif 000026
Ｎ個の状態ベクトルｘiバーには、  ＩＭＭ−ＥＫＦブロック２０において初期化すべきパラメータが３個存在する。これら３個のパラメータは、例えば最初の１．５秒間センサ信号を利用して推定可能である。これらのパラメータ値をランダムに設定することもできるが、選択する値が良くない場合には収束時間が長くなる。
【００９４】
図１０では、周波数の推定に２段階のステップが必要なことを示す。すなわち全てのセンサに対する周波数推定そのもの、および使用する値の決定である。
【００９５】
推定には、部分空間アプローチ、すなわち信号の固有値および固有ベクトルによる分解を使用している。ＥＳＰＲＩＴ（Estimation of Signal Parameters via Rotational Invariance Techniques、回転不変技術による信号パラメータ推定）を用いた。原理を図１１に示す。
【００９６】
ＥＳＰＲＩＴによる周波数推定では、部分空間の間での決定論的な関係を利用している。最初にデータ行列Ｘを構築する必要があり、次のように実施される。
【数２６】
0006738177.tif 000027
ここでｘはセンサ信号、Ｄはウィンドウ・サイズである。例えばＤ＝８である。Ｏは使用したサンプルの数である。
【００９７】
次に行列Ｘに対して特異値分解（ＳＶＤ）を適用し、Ｘ＝ＬＳＵというように書き換えることが可能である。ここでＬは左特異ベクトルの［Ｏ×Ｏ］行列、Ｕは右特異ベクトルの［Ｄ×Ｄ］行列である。Ｓは主対角上に大きさが降順に並んだ特異値の［Ｏ×Ｄ］行列である。
【００９８】
Ｕは、多次元ベクトル空間の正規直交基底を形成している。この部分空間は、信号部分空間（Ｕｓ）とノイズ部分空間（Ｕｎ）とに分割可能である。部分空間の間の閾値は、Ｐ＝５に設定されている。このことはＵｓが、右側に特異値があり、Ｐの値が最大となる行列であることを意味している。
【００９９】
次のステップでは、部分空間を分離して、Ｕ1とＵ2とに分離する。Ｕ1には１からＤ−１までの要素を含まれ、Ｕ2には２からＤまでの要素が含まれる。分離した部分空間との間の回転特性により、周波数が推定できる。
【０１００】
続いて次式を用いてΨを計算する。
【数２７】
0006738177.tif 000028
周波数推定方法を図１１に示す。ここでΦpにはΨの固有値が含まれている。
【０１０１】
全てのセンサ信号に対して上記を実行した後、各センサに対して得られた周波数推定値をクラスタリングして決定を下す。
【０１０２】
図１２に示すように、前述の周波数推定手順により得られた全てのセンサ（この場合１４個のセンサ）周波数推定値に対してクラスタリングを繰り返し実行する。基本的に、距離ｄがクラスタの重心（クラスタ内の全ての周波数推定値の平均値）に近い場合、かかるクラスタに新しい周波数推定値を追加する。そうでない場合には、新しいクラスタを作成する。
【０１０３】
最後に要素の数が最も多いクラスタを選択し、初期周波数は本クラスタ中の周波数推定値の平均値となる。
【０１０４】
振幅と位相は、脈拍推定に影響を与えないように０（ゼロ）に設定してもよい。詳細な推定値が必要であれば、最小二乗法によるパラメータ推定が利用可能である。
【０１０５】
これらの値をＮ＝１４の全ての状態ベクトルに対して設定する。
【０１０６】
２．３  プロセス共分散およびノイズ共分散の初期化
プロセスノイズは、状態パラメータに関係したノイズである。プロセスノイズを定義するためには２つのパラメータ、すなわち周波数推定ノイズｄfおよび振幅推定ノイズｄAである。
【０１０７】
ＩＭＭ−ＥＫＦ脈拍推定に対しては、ｄAはそれほど重大な影響を及ぼさない。単にｄA＝５×１０-11（実験的に見出した）に設定してもよい。
【０１０８】
ＩＭＭ−ＥＫＦ脈拍推定に対してｄfは大きな影響力を持ち、脈拍の良好な追跡を可能にする一方、間違えば追跡が困難になる。ｄfの値は、例えば各センサの初期１．５秒の信号を使用して、オンラインで計算される。この値は各センサごとに異なる可能性がある。
【０１０９】
センサ共分散値とクラメル・ラオの下限計算に基づいたｄfの最適値との間で、非線形モデルが見出された。これらの最適値は、脈拍を推定する際にｄfに対して相異なる値を設定してＥＫＦにより得られる結果を確認することにより、経験的にも容易に見出すことが可能である。最良の脈拍推定が得られる場合が、最良のｄf値に対応する。
【０１１０】
換言すれば、入力（観測共分散）および既知のｄf最適値をフィッティングするモデルが求められる。
【０１１１】
多項式によるフィッティングを採用すれば、誤差は大きくなる。ニューラルネットワークにより、複雑な非線形性がモデル化可能である。図１３に、ニューラルネットワークを用いた非線形フィッティング手順を示す。このニューラルネットワークは、図６に示したセンサ自動変更用のものと若干異なる。ここで出力層は単なるマッピングである。出力層には１個のニューロン、中間層には１５個のニューロンが存在する。入力と出力は、それぞれ１つだけのパラメータを有する。モデルパラメータの計算はレーベンバーグ・マルカート法を使用して実行する。
【０１１２】
一旦モデルが計算されると、例えば初期１．５秒間のセンサ信号の間にオンラインで使用可能である。計算はセンサ自動変更のものと非常に類似しており、以下の式で示す。
【数２８】
0006738177.tif 000029
ここでｘoffsetは観測した共分散Ｒから取り除くべきオフセット、Ｇはゲインである。これらオフセットとゲインは、前述のオフライン手順により計算済である。
【０１１３】
次に中間層出力が次式のように計算される。
【数２９】
0006738177.tif 000030
ここでＢhiddenは中間層のバイアスベクトルであり、中間層のニューロン数、Ｑに等しい長さを有する。
【０１１４】
Ｗhiddenは中間層の係数ベクトルであり、そのサイズは［Ｑ×１］である。ｙsigmoidは非線形のシグモイド関数であり、その式は前述している。中間層出力であるｙhiddenは、長さＭのベクトルである。
【０１１５】
最後のステップとして、次式で示す中間層出力を計算する。
【数３０】
0006738177.tif 000031
ここでＢoutは、長さが１の出力層バイアスベクトルである。Ｗoutは、出力層の係数ベクトルであり、そのサイズは［１×Ｑ］である。出力のｄfはＩＭＭ−ＥＫＦブロック２０で直接使用される。
【０１１６】
残りのノイズ共分散推定値は、現タイムフレームでのセンサ信号の分散である。
【０１１７】
２．４  ＩＭＭ−ＥＫＦ推定
ｘiバーからｘNバーまでのＮ個の状態ベクトル、および周波数と振幅と位相との３つのパラメータが存在し、モデル２１、２２の全てにおいて相異なっている。全てのセンサに対するＩＭＭ−ＥＫＦの一般式は、次式で与えられる。
【数３１】
0006738177.tif 000032
ここでｘkはサンプルｋfにおける線形状態空間方程式、ｙkは非線形観測方程式である。ｖkおよびｗkは、それぞれプロセスノイズおよび観測ノイズであり、指数のｉはセンサ番号を示す。Ｆk-1は時刻ｋ−１におけるサンプルの線形状態空間モデルであり、ｈk-1は時刻ｋ−１におけるサンプルの非線形観測モデルである。
【０１１８】
Ｆは線形状態空間遷移行列であり、本ケースでは次のように表される。
【数３２】
0006738177.tif 000033
状態空間行列は、図１に示すモデル切替ブロック２３によりモデル切り替え時に変更することも可能であるが、前述のハーストによるセンサ自動変更においては必要としない。
【０１１９】
例えば観測方程式は次のように書き換えることが可能である。
【数３３】
0006738177.tif 000034
観測方程式は、特許文献１−３に説明されているように、正弦波の和として表すことが可能である。ＩＭＭ−ＥＫＦでは３つのステップがあり、それらの方程式を以下に示す。
【０１２０】
予測
【数３４】
0006738177.tif 000035
ここでｘk|k-1バーは、前のパラメータ値に基づいた現状態パラメータ予測値である。ｘk-1|k-1バーは前の状態パラメータ、Ｐk-1|k-1バーは前の共分散に基づいた現共分散である。Ｑはプロセスノイズ共分散である。
【０１２１】
更新
【数３５】
0006738177.tif 000036
ＳkiおよびＫkiは次式で表される。
【数３６】
0006738177.tif 000037
ｈは非線形関数なので、線形化する必要がある。したがって、Ｈikバーは、センサｉに対応する非線形関数ｈを局所的に線形化したものであり、ｘk|k-1iバーにおいて評価したヤコビアンとして定義される。前述のケースにおいては、次式のように表される。
【数３７】
0006738177.tif 000038
【０１２２】
モデル切替
全てのセンサに対する状態ベクトル推定値が求まったので、センサ自動切替ブロック３０の決定に従い、最終の状態パラメータ推定値は次のように表される。
【数３８】
0006738177.tif 000039
ここでｐはセンサ自動切替ブロック３０で予測されたセンサ番号である。
【０１２３】
図１４は、１４個（Ｎ＝１４の場合）のモデル２１、２２を有しかつモデル切替が可能な場合、ＩＭＭ−ＥＫＦブロック２０のＩＭＭ−ＥＫＦ処理方法を示す。図１に示した発明の実施形態によれば、モデル切替ブロック２３がモデル２１とモデル２２の切替を実行する。黒点が最終状態推定値であり、図１に示す生体パラメータ推定＆追跡ブロック２４により、所定のタイムフレームＴごとに脈拍周波数推定値を提供するために使用される。図１４に示す「フィルタ」ブロックはＩＭＭ−ＥＫＦブロック２０による更新処理に対応しており、生体パラメータ推定および追跡ブロック２４において更新処理を実行する。
【０１２４】
センサ自動切替とＩＭＭ−ＥＫＦとを組合せることにより、身体が大きく動く場合、さらに一般的にあらゆる身体動作においても、頑健な脈拍推定が達成される。
【０１２５】
図１５は脈拍推定結果を示す。ｘ軸は時間、ｙ軸は毎分の心臓鼓動数を示す。破線のボックスは、車が動いている（高速ターンをしながら運転中）場合に対応している。その他は車が静止している場合の結果である。１０％の許容差も示す（外側の一点鎖線の曲線）。中央の破線が実際の脈拍値を、実線が推定結果を示す。
【０１２６】
図１５からわかるように、脈拍推定は非常に精確であり、脈拍推定結果に見られる振動は、脈拍に伴う呼吸変調を表している。
【０１２７】
本発明によれば、センサ自動変更により、ノイズと身体動作依存とを考慮しつつ、タイムフレームＴごとに脈拍推定に使用する一つのセンサ（好ましくは考えられるうちで最良のセンサ）を予測して決定する。
【０１２８】
さらに本発明によれば、センサ変更情報を使用しての脈拍推定は、必要だと判断した場合、モデルを切り替えるＩＭＭ−ＥＫＦの内部で実行する。
【０１２９】
以上の実施形態によれば、２個以上のセンサをそれぞれ圧力変化測定用として有し、さらに１個以上の加速度計を接続した支持部材上の対象の生体パラメータを推定する。２個以上のセンサからの信号はそれぞれかかるセンサにより測定される圧力変化に対応しており、かかる信号に基づいて各２個以上のセンサにセンサ固有モデルを設定する。選択プロセスにおいては、タイムフレームＴごとに、２個以上のセンサから１個のセンサを選択し、１個以上の加速度計からの信号に基づいて対象の生体パラメータを推定するために使用する。推定プロセスにおいては、タイムフレームＴごとに、選択したセンサに設定したセンサ固有のモデルを用いて対象の生体パラメータを推定する。
【０１３０】
図１に示す装置１の機能は、適切な場合には、ソフトウェア、ファームウェア、ハードウェアの任意の組み合わせまたはいずれかによって実現してもよい。一般的には、本発明の実施形態は、メモリーに格納されプロセッサにより実施可能なコンピュータソフトウェアにより実現されるか、またはハードウェアにより実現するか、またはソフトウェアとファームウェアの両方もしくはいずれかとハードウェアとの組み合わせにより実現される。
【０１３１】
メモリはローカルの技術環境に適した任意の種類のものでよく、半導体ベースのメモリ装置、磁気メモリ装置およびシステム、光メモリデバイスおよびシステム、固定メモリ、および取り外し可能メモリなど、適切なデータストレージ技術を使用して実現することができる。プロセッサはローカルの技術環境に適した任意の種類のものでよく、非限定的な例として汎用コンピュータ、特殊目的コンピュータ、マイクロプロセッサ、デジタル信号プロセッサ（ＤＳＰ）、およびマルチコアプロセッサアーキテクチャに基づくプロセッサのうち、いずれか一つまたは複数を含んでもよい。
【０１３２】
さらにこの点に関して、前述の様々な論理ステップは、プログラムのステップ、または論理回路と論理ブロックと論理機能とを相互接続したもの、またはプログラムのステップおよび論理回路と論理ブロックと論理機能との組み合わせにより実現してもよい。
【０１３３】
上記の説明は本発明を例示するものであり、本発明を限定するものとして解釈してはならないということを理解しておく必要がある。技術に優れたものであれば、添付の特許請求の範囲によって定義される本発明の適用範囲から逸脱することなく、様々な修正および応用を考え得る可能性がある。
【図１】
0006738177.tif 000040
【図２】
0006738177.tif 000041
【図３】
0006738177.tif 000042
【図４】
0006738177.tif 000043
【図５】
0006738177.tif 000044
【図６】
0006738177.tif 000045
【図７】
0006738177.tif 000046
【図８】
0006738177.tif 000047
【図９】
0006738177.tif 000048
【図１０】
0006738177.tif 000049
【図１１】
0006738177.tif 000050
【図１２】
0006738177.tif 000051
【図１３】
0006738177.tif 000052
【図１４】
0006738177.tif 000053
【図１５】
0006738177.tif 000054