JPWO2017110753A1 - 個数分析方法、個数分析装置および個数分析用記憶媒体 - Google Patents

Abstract

【課題】粒子状または分子状の分析物の種別に応じた個数ないし個数分布を高精度に分析することのできる個数分析方法、個数分析装置および個数分析用記憶媒体を提供することである。
【解決手段】検体の粒子の貫通孔１２の通過に応じて、ナノポアデバイス８によって検出した粒子通過検出信号のデータ群をベースにして、コンピュータ制御プログラムを実行することにより、個数導出手段である粒子種分布推定プログラムの実行によって、該粒子通過検出信号として得られる粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出することができる。
【選択図】図２

Description

本発明は、例えばウイルス、細菌等の微小物体、あるいは微細な粉塵等の粒子状または分子状の分析物の個数ないし個数分布を分析する個数分析方法、個数分析装置および個数分析用記憶媒体に関する。

従来より、細菌などを対象にした微生物検査は生化学的手法によって行われている。生化学的手法による検査においては、細菌の数や種類の同定のために培養・染色が行われる。

ＷＯ２０１３−１３７２０９号公報

上記従来の生化学的手法による検査においては、検査時間が数日程度（例えば、大腸菌の培養時間は１〜２日）要し、また検査実施者にも専門的な技術も要求されていた。このため、人やモノの行来が増大する現代社会においては、食の安全、パンデミック防止、ＰＭ２．５のような微小粒子状物質による大気汚染による健康被害の防止のためには、その場で誰もが実施可能で、迅速かつ簡便に低コストで行える検査手法が求められている。

特許文献１には、マイクロ・ナノメートルスケールの微細貫通孔（マイクロポア、以下ナノポアと総称する。）を有するマイクロ・ナノポアデバイスを用いた微小物体（細菌、ウイルス等）の電気的検出技術が開示されている。

マイクロ・ナノポアデバイスの動作原理は以下の通りである。ナノポアの上下と、ナノポアは電解質溶液で充填された状態で、ナノポアを挟むようにして配置された電極対に電圧印加を行うと、ポア径、イオン濃度および印加電圧に比例し、かつポア深さに反比例した一定の電流が計測される。細菌などの検体（分析物）がポア（貫通孔）を通過する際、一部のイオン電流が検体により阻害されるため、パルス状の電流変化が現れる。この電流変化を観測することにより、電解質溶液中に存在する検体を検出することができる。

溶液中の検体の種類が判明している場合には、該電流変化の検出数が検体の総数となる。しかしながら、実際の検査においては、種別不明の分析物を検査対象とする場合があるため、ナノポアデバイスを実用的に使用するうえで、単に電流変化の抽出のみでは分析物種別の詳しい検査に適用できないという問題があった。

本発明の目的は、粒子状または分子状の分析物の種別に応じた個数ないし個数分布を高精度に分析することのできる個数分析方法、個数分析装置および個数分析用記憶媒体を提供することである。

本発明は、上記課題に鑑み、ナノポアがポア径に対してポア厚が十分に小さい低アスペクトの場合には、検出した阻害電流の波形には通過する分析物の形態が反映される点に着目し、阻害電流波形の特徴を数学的に抽出して得られるデータに対して統計的コンピュータ解析を加えることにより分析物種別に応じた個数ないし個数分布を導出し得るという知
見に基づいてなされた発明である。

本発明に係る第１の形態は、
貫通孔を形成した隔壁と、前記貫通孔を介して前記隔壁の表裏側に配置された電極とを配置し、
粒子状または分子状の分析物を含む流動性物質を前記隔壁の一面側に供給し、
前記分析物が前記貫通孔を通過することにより生ずる電極間の通電変化を検出した検出信号のデータに基づいてコンピュータ制御プログラムの実行によって分析物種別の個数を分析する個数分析方法であって、
前記コンピュータ制御プログラムは、
前記検出信号として得られた分析物通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、分析物種別の個数を導出する個数導出手段を少なくとも有することを特徴とする個数分析方法である。

本発明に係る第２の形態は、
前記特徴量は、
前記パルス状信号の波形の局所的特徴を示す第１類型と、
前記パルス状信号の波形の全体的特徴を示す第２類型のいずれかである個数分析方法である。

本発明に係る第３の形態は、
前記第１類型の特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さを表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかである個数分析方法である。

本発明に係る第４の形態は、
前記第２類型の特徴量は、
パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかである個数分析方法である。

本発明に係る第５の形態は、
前記コンピュータ制御プログラムは、
前記検出信号のデータないしそれに含まれる揺らぎ成分から、分析物通過のないときの基準ラインを抽出する基準ライン抽出手段と、
前記基準ラインを基準にして所定範囲を超える信号データを前記パルス状信号のデータとして抽出するパルス抽出手段と、
抽出した前記パルス状信号のデータから前記特徴量を抽出する特徴量抽出手段と、を含む個数分析方法である。

本発明に係る第６の形態は、
貫通孔を形成した隔壁と、前記貫通孔を介して前記隔壁の表裏側に配置された電極とを配置し、
粒子状または分子状の分析物を含む流動性物質を前記隔壁の一面側に供給し、
前記分析物が前記貫通孔を通過することにより生ずる電極間の通電変化を検出した検出信号のデータに基づいてコンピュータ制御プログラムの実行によって分析物種別の個数を分析する個数分析装置であって、
前記検出信号のデータを記憶する記憶手段と、
前記検出信号のデータに基づく前記コンピュータ制御プログラムの実行制御を行う制御手段と、を有し、
前記コンピュータ制御プログラムは、
前記検出信号として得られた分析物通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、分析物種別の個数を導出する個数導出手段を少なくとも有することを特徴とする個数分析装置である。

本発明に係る第７の形態は、
前記特徴量は、
前記パルス状信号の波形の局所的特徴を示す第１類型と、
前記パルス状信号の波形の全体的特徴を示す第２類型のいずれかである個数分析装置である。

本発明に係る第８の形態は、
前記第１類型の特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さを表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかである個数分析装置である。

本発明に係る第９の形態は、
前記第２類型の特徴量は、
パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクト
ルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかである個数分析装置である。

本発明に係る第１０の形態は、
前記コンピュータ制御プログラムは、
前記検出信号のデータないしそれに含まれる揺らぎ成分から、分析物通過のないときの基準ラインを抽出する基準ライン抽出手段と、
前記基準ラインを基準にして所定範囲を超える信号データを前記パルス状信号のデータとして抽出するパルス抽出手段と、
抽出した前記パルス状信号のデータから前記特徴量を抽出する特徴量抽出手段と、を含む個数分析装置である。

本発明に係る第１１の形態は、前記個数導出手段により導出した個数データを分析物種別に所定の出力態様で出力する出力手段を有する個数分析装置である。

本発明に係る第１２の形態は、第１の形態に係る前記コンピュータ制御プログラムを記憶したことを特徴とする個数分析用記憶媒体である。

第１の形態に係る個数分析方法によれば、上記のナノポアデバイスにより検出した検出信号のデータ群をベースにして、前記コンピュータ制御プログラムを実行することにより、前記個数導出手段によって、該検出信号として得られた分析物通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、例えば、細菌、微小粒子状物質、分子状物質等の分析物の個数を導出することができるので、該分析物の種別に応じた個数ないし個数分布を高精度に分析することができ、分析検査における簡易化および低コスト化を実現することができる。

第２の形態によれば、前記特徴量は、パルス状信号由来のパラメータとして、該パルス状信号の波形の局所的特徴を示す第１類型と、該パルス状信号の波形の全体的特徴を示す第２類型のいずれかであるので、これらのうち１または２以上の特徴量を使用して個数分析を行うことによって、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、分析検査における簡易化および低コスト化に寄与することができる。

本発明に係る個数分析方法において、第１類型または第２類型の特徴量の少なくとも１つ以上の特徴量を使用して個数分析する場合に限らず、第１類型の特徴量と第２類型の特徴量をそれぞれ少なくとも１つ以上、組み合わせて使用して個数分析を行うことができる。

第３の形態によれば、
前記第１類型の特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さを表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかであるので、これらのうち１または２以上の特徴量を使用して個数分析を行うことによって、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、分析検査における簡易化および低コスト化に寄与することができる。

第４の形態によれば、
前記第２類型の特徴量は、
パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかであるので、これらのうち１または２以上の特徴量を使用して個数分析を行うことによって、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、分析検査における簡易化および低コスト化に寄与することができる。

第５の形態によれば、前記基準ライン抽出手段により、前記検出信号のデータないしそれに含まれる揺らぎ成分から、分析物通過のないときの基準ラインを抽出し、前記パルス抽出手段により、前記基準ラインを基準にして所定範囲を超える信号データを前記パルス状信号のデータとして抽出し、前記特徴量抽出手段により、抽出した前記パルス状信号のデータから特徴量を抽出するので、パルス状信号由来の特徴量に基づいて個数分析を実行して、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、分析検査における簡易化および低コスト化に寄与することができる。

第６の形態によれば、第１の形態に係る個数分析方法に基づく個数分析をコンピュータ解析により実行可能であるので、第１形態で説明したコンピュータ制御プログラムによる全ての効果を有し、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、簡便かつ安価に個数分析を行える個数分析装置を提供することができる。

第７の形態によれば、前記特徴量は、パルス状信号由来のパラメータとして、該パルス状信号の波形の局所的特徴を示す第１類型と、該パルス状信号の波形の全体的特徴を示す第２類型のいずれかであるので、これらのうち１または２以上の特徴量を使用して個数分析を行うことによって、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、簡便かつ安価に個数分析を行える個数分析装置を提供することができる。

本発明に係る個数分析装置において、第１類型または第２類型の特徴量の少なくとも１つ以上の特徴量を使用して個数分析する場合に限らず、第１類型の特徴量と第２類型の特
徴量をそれぞれ少なくとも１つ以上、組み合わせて使用して個数分析を行うことができる。

第８の形態によれば、
前記第１類型の特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さを表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかであるので、これらのうち１または２以上の特徴量を使用して個数分析を行うことによって、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、簡便かつ安価に粒子種分析を行える個数分析装置を提供することができる。

第９の形態によれば、
前記第２類型の特徴量は、
パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかであるので、これらのうち１または２以上の特徴量を使用して個数分析を行うことによって、分析物種別に応じた個数ないし個数分布を高精度に分析することができ、簡便かつ安価に粒子種分析を行える個数分析装置を提供することができる。

第１０の形態によれば、前記基準ライン抽出手段により、前記検出信号のデータないしそれに含まれる揺らぎ成分から、分析物通過のないときの基準ラインを抽出し、前記パルス抽出手段により、前記基準ラインを基準にして所定範囲を超える信号データを前記パルス状信号のデータとして抽出し、前記特徴量抽出手段により、抽出した前記パルス状信号のデータから特徴量を抽出するので、パルス状信号由来の特徴量に基づいて個数分析を実行して、分析物種別に応じた個数ないし個数分布を高精度にコンピュータ分析することができ、簡便かつ安価に個数分析を行える個数分析装置を提供することができる。

第１１の形態によれば、前記個数導出手段により導出した個数データを分析物種別に所定の出力態様で出力することができるので、例えば、検査現場において即座に検査結果を表示出力するなどして、迅速な検査結果を報知することができ、迅速簡便かつ安価に個数分析を行える個数分析装置を提供することができる。

第１２の形態によれば、第１の形態に係る前記コンピュータ制御プログラムを記憶した個数分析用記憶媒体を提供することができる。したがって、本形態に係る記憶媒体は、第１の形態で説明した前記コンピュータ制御プログラムによる効果を有するので、個数分析用記憶媒体に記憶したコンピュータ制御プログラムをコンピュータにインストールして該コンピュータに個数分析動作させることによって簡便かつ安価に個数分析を行うことができる。

本発明における記憶媒体としては、フレキシブルディスク、磁気ディスク、光ディスク、ＣＤ、ＭＯ、ＤＶＤ、ハードディスク、モバイル端末等、コンピュータにより読み取り可能な記憶媒体のいずれかを選択することができる。

本発明によれば、コンピュータ端末を利用して、例えば、細菌、微小粒子状物質、分子状物質等の分析物の個数ないし個数分布の分析を安価で、簡易かつ高精度に行うことができる。

本発明の一実施形態である個数分析装置の概略ブロック図である。 マイクロ・ナノポアデバイスの概略構成を示す概略側断面図である。 前記個数分析装置のＰＣ１により実行可能な個数分析処理の説明に必要な処理プログラム構成を示す図である。 実施例の大腸菌と枯草菌につき実測した粒子通過によるパルス波形例を示す図である。 本発明に係る各種特徴量を説明するためのパルス波形図である。 カルマンフィルターを説明するための図である。 カルマンフィルターの各因子を実際の計測電流データで説明するための図である。 カルマンフィルターの予測（８Ａ）と更新（８Ｂ）の繰返しの詳細を示す図である。 ＢＬ推定処理プログラムに基づくＢＬ推定処理を示すフローチャートである。 カルマンフィルターの因子調整に使用したビーズモデルの波形図である。 大腸菌２２と枯草菌２３が電解質溶液２４中に混在する様子を模式的に示した貫通孔１２周辺の拡大図である。 調整因子のｍ、ｋ、αの組合せに応じてビーズモデルの波形から拾われたパルスの数を示す表である。 特徴量抽出プログラムの実行処理内容の概要を示すフローチャートである。 粒子種推定処理を示すフローチャートである。 １つの波形データに関する各特徴量（１５Ａ）と、大腸菌と枯草菌の粒子種における確率密度関数のイメージ図（１５Ｂ）とを示す図である。 大腸菌と枯草菌の粒子種の個々より得られた確率密度分布の重ね合わせのイメージ図である。 ｋ個の粒子種別の粒子総数と、粒子種別の出現確率と、データ全体の出現頻度の期待値との関係を示すイメージ図である。 ラグランジュ未定乗数法により最適化を行う制約付き対数尤度最大化式の導出過程を説明するための図である。 データファイル作成処理を示すフローチャートである。 確率密度関数の推定処理を示すフローチャートである。 粒子数の推定処理を示すフローチャートである。 Ｈａｓｓｅｌｂｌａｄ反復法による粒子数推定処理を示すフローチャートである。 ＥＭアルゴリズムによる処理手順を示すフローチャートである。 本実施形態に係る個数分析装置により分析した結果の一例を示す図である。 特徴量としてパルス波長、波高を使用した検証例と、特徴量としてパルス波長、ピーク位置比を使用した検証例の各推定結果データを示す表である。 特徴量としてピーク付近波形の広がり、パルス波長を使用した検証例と、特徴量としてピーク付近波形の広がり、波高を使用した検証例の各推定結果データを示す表である。 特徴量として尖度と、パルス波高を使用した場合における個数推結果を示す図である。 ＢＬ推定処理プログラムに基づくＢＬ推定処理を示すフローチャートである。 大腸菌と枯草菌の混合比をそれぞれ、１：１０、２：１０、３：１０、３５：１００とした場合における各個数推定結果を示すヒストグラムである。 大腸菌と枯草菌の混合比をそれぞれ、４：１０、４５：１００、１：２とした場合における各個数推定結果を示すヒストグラムである。 特徴量としてパルス波長、パルス波高を使用した場合における各粒子の散布状態を合成した図である。 特徴量としてピーク付近波形の広がり、パルス波長を使用した場合、特徴量としてピーク付近波形の広がり、ピーク位置比を使用した場合、ピーク付近波形の広がり、パルス波高を使用した場合における各粒子の散布状態を合成した図である。 マイクロ・ナノポアデバイス８を用いて、３種の粒子３３ａ、３３ｂ、３３ｃが貫通孔１２を通過して得られる検出信号の波形例と、特徴量に基づいて得られる確率密度関数の導出例を示す図である。 俯角および面積の特徴量を説明するためのパルス波形図である。 波高ベクトルの取得の仕方を説明するための図である。 ｄ次元の波高ベクトルとデータサンプリングとの関係を説明するための図である。 時間（波長）および波幅に関する慣性モーメントの特徴量を説明するためのパルス波形図である。 ｄ_w次元の波幅ベクトルとデータサンプリングと関係を説明するための図である。 波幅に関する慣性モーメントを波幅ベクトルにより取得する取得過程を説明するための図である。 複数の方向に分割した場合の特徴量作成用波形ベクトルの一例を説明するための図である。 特徴量抽出の処理内容を示すフローチャートである。 １ＭＨｚ、５００ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である ２５０ｋＨｚ、１２５ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 ６３ｋＨｚ、３２ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 １６ｋＨｚ、８ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 ４ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価表である。 全サンプリングデータに対する各特徴量組合せに関する推定評価表である。 １ＭＨｚ〜１２５ｋＨｚでの高密度サンプリングしたときの各特徴量組合せに関する推定評価表である。 ６３ｋＨｚ〜４ｋＨｚでの低密度サンプリングしたときの各特徴量組合せに関する推定評価表である。 全サンプリングデータを使用したとき（５０Ａ）および高密度にサンプリングしたとき（５０Ｂ）に高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフである。 低密度にサンプリングしたときに高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ａ）と、全サンプリングデータを使用したときの４種類の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ｂ）である。 ４種類の各特徴量組合せに対する、特徴量作成に要する計算時間とＨａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間との合計計算時間を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ａ）と、各特徴量組合せに対する特徴量作成に要する計算時間を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ｂ）である。 ４種類の各特徴量組合せに対する、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間を示すサンプリング周波数−所要計算時間（秒）のグラフである。

本発明の一実施形態に係る個数分析装置を図面を参照して以下に説明する。本実施形態においては、分析物の一例として細菌等の微生物粒子を分析する粒子個数分析形態で説明する。

図１は本実施形態の個数分析装置の概略構成を示す。この個数分析装置は、パーソナルコンピュータ（以下、ＰＣという。）１により構成され、ＰＣ１にはＣＰＵ２、ＲＯＭ３、ＲＡＭ４およびデータファイル記憶部５を有する。ＲＯＭ３には本発明に係る個数分析プログラムが格納されている。個数分析プログラムは、同プログラムを記憶した記憶媒体（ＣＤ、ＤＶＤ等）からインストールされて格納されている。ＰＣ１にはキーボード等の入力手段６および液晶ディスプレイ等の表示手段７が入出力可能に接続されている。データファイル記憶部５には個数分析用データが格納可能になっている。

図２はマイクロ・ナノポアデバイス８を用いた粒子検出装置の概略構成を示す。

粒子検出装置は、マイクロ・ナノポアデバイス８およびイオン電流検出部により構成されている。マイクロ・ナノポアデバイス８は、チャンバー９と、チャンバー９を上下の収容空間に区画する隔壁１１と、隔壁１１の表裏側に配置された一対の電極１３、１４とを有する。隔壁１１は基板１０上に形成されている。隔壁１１の中央付近には、微小径の貫通孔１２が穿設されている。貫通孔１２の下方には、基板１０の一部を下向きに凹状に取り除いた凹部１８を設けている。

マイクロ・ナノポアデバイス８は半導体デバイス等の製造技術（例えば、電子線描画法やフォトリソグラフィ）を用いて作成される。すなわち、基板１０はＳｉ材で構成され、表面上にＳｉ₃Ｎ₄膜による隔壁１１が薄膜形成されている。凹部１８は基板１０の一部をエッチングにより除去して形成されている。

隔壁１１は、大きさ１０ｍｍ角で、厚さ０．６ｍｍのＳｉ基板に５０ｎｍのＳｉＮ膜を積層して形成されている。Ｓｉ₃Ｎ₄膜にレジストを塗布し、電子線描画法により、直径３
μｍの円形の開口パターンを形成して、貫通孔１２が穿設されている。貫通孔１２の裏側においては、ＫＯＨによるウェットエッチングを行い５０μｍ角の開口を形成して凹部１８を設けている。凹部１８の形成は、ウェットエッチングに限らず、例えば、ＣＦ₄系ガスによるドライエッチングなどによる等方的なエッチング等により行うことができる。

隔壁１１用の膜には、ＳｉＮ膜の他に、ＳｉＯ₂膜、Ａｌ₂Ｏ₃膜、ガラス、サファイア、セラミック、樹脂、ゴム、エラストマー等の絶縁性膜を使用することができる。基板１０の基板材料もＳｉに限定されるものではなく、ガラス、サファイア、セラミック、樹脂、ゴム、エラストマー、ＳｉＯ₂、ＳｉＮ、Ａｌ₂Ｏ₃等を使用することができる。

貫通孔１２は、上記の基板上の薄膜に形成する場合に限らず、例えば、貫通孔１２を形成した薄膜状シートを基板上に接着することによって、貫通孔を有する隔壁を形成するようにしてもよい。

イオン電流検出部は、電極１３、１４の電極対と、電源１５、増幅器１６および電圧計２０により構成されている。電極１３、１４は貫通孔１２を介して対向配置されている。増幅器１６はオペアンプ１７と帰還抵抗１９とにより構成されている。オペアンプ１７の（−）入力端子と電極１３は接続されている。オペアンプ１７の（＋）入力端子は接地されている。オペアンプ１７の出力側と電源１５の間に電圧計２０が接続配置されている。電源１５により電極１３、１４間には、０．０５〜１Ｖの印加電圧を使用できるが、本実施例では０．０５Ｖが印加されるようになっている。増幅器１６は電極間に流れる電流を増幅して電圧計２０側に出力する。電極１３、１４の電極材料としては、例えば、Ａｇ／ＡｇＣｌ電極、Ｐｔ電極、Ａｕ電極等を使用でき、好ましくはＡｇ／ＡｇＣｌ電極である。

チャンバー９は、マイクロ・ナノポアデバイス８周囲を密閉状に囲む流動性物質収容容器であり、電気的および化学的に不活性な材質、例えば、ガラス、サファイア、セラミック、樹脂、ゴム、エラストマー、ＳｉＯ₂、ＳｉＮ、Ａｌ₂Ｏ₃等により形成されてよい。

チャンバー９内には注入口（図示せず）から検体２１を含む電解質溶液２４が充填される。検体２１は、例えば、細菌、微小粒子状物質、分子状物質等の分析物である。検体２１を、流動性物質である電解質溶液２４に混入して、マイクロ・ナノポアデバイス８による計測が行われる。イオン電流検出部による計測終了時には排出口（図示せず）より充填溶液は排出可能になっている。電解質溶液には、例えば、リン酸緩衝生理食塩水（ＰＢＳ）、Ｔｒｉｓ−ＥＤＴＡ（ＴＥ）バッファーやそれらの希釈液の他、これらと同様なすべての電解質溶液剤を使用することができる。計測は、検体含有電解質溶液をチャンバー９内に導入、充填するごとに行う場合に限らず、溶液溜から簡易ポンプ装置により検体含有電解質溶液（流動性物質）を汲み出して注入口よりチャンバー９内に充填し、計測後に排出口から排出させ、さらに、別の溶液溜あるいは新たな溶液を溶液溜に貯留し、新たに汲み出して次の計測を行う連続計測システムを構成するようにしてもよい。

電解質溶液２４をチャンバー９内に充填した状態で、貫通孔１２の上下の電極１３、１４間に電源１５の電圧印加を行うと、貫通孔１２に比例した一定のイオン電流が電極間に流れる。電解質溶液２４中の細菌等の検体が貫通孔１２を通過する際には、一部のイオン電流が検体により阻害されるため、電圧計２０によりパルス状のイオン電流減少を計測することができる。したがって、マイクロ・ナノポアデバイス８を用いた粒子検出装置によれば、計測電流の波形変化を検出することにより、検体（例えば、粒子）毎の貫通孔１２通過による、流動性物質中に含まれる粒子個々の存在を高精度に検出することができる。計測態様には、流動性物質を強制的に流動させながら計測する場合に限らず、流動性物質を非強制的に流動させながら計測する場合を含むことができる。

電圧計２０によるイオン電流の計測出力は外部出力可能になっている。この外部出力は、変換回路装置（図示せず）によりデジタル信号データ（計測電流データ）に変換されて記憶装置（図示せず）に一旦保存された後、データファイル記憶部５に格納される。

本実施形態に係る個数分析装置は、分析対象として１種または２種以上の粒子（分析物の一例）を含む流動性物質（電解質溶液２４）を隔壁１１上側の一面側に供給し、粒子が貫通孔１２を通過することにより生ずる電極１３、１４間の通電変化を検出した検出信号のデータ（計測電流データ）に基づいてコンピュータ制御プログラム（個数分析プログラム）の実行によって粒子種別の個数ないし個数分布を分析する個数分析方法であって、該検出信号に含まれ粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出する個数分析方法に基づいて、粒子種別個数の自動分析を行う。ＰＣ１は、ＣＰＵ２の制御によりＲＯＭ３に格納した個数分析プログラムを実行することにより、データファイル記憶部５に格納、記憶した計測電流データに対する個数分析処理を行うことができる。

図３は，ＰＣ１により実行可能な個数分析処理の説明に必要な処理プログラム構成を示す。分析対象のデータ実施例としては、分析物としての２種の粒子（大腸菌と枯草菌）を含む電解質溶液２４を用いて抽出した計測電流データ（各粒子のパルス抽出データ）を元データに使用している。

個数分析プログラムには、検出信号として得られた、貫通孔１２の粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度関数を求める確率密度関数モジュールプログラムと、確率密度推定の結果から粒子種別の個数を導出する粒子種分布推定プログラムとを含む。さらに、個数分析プログラムには、データ群から抽出したベースラインを基準にして、パルス状信号の波形形態の特徴を示す特徴量を抽出する特徴量抽出プログラムと、抽出した特徴量に基づいて得られる粒子毎のパルス特徴量データによるデータファイルを作成するデータファイル作成プログラムとを含む。確率密度推定処理および粒子種別の個数導出処理は、データファイル作成プログラムにより作成されたデータファイルから取込んだ特徴量のスカラーデータに対して実行可能になっている。特徴量抽出プログラムには、元の計測電流データから該ベースラインを抽出するベースライン推定処理プログラムを含む。

粒子種分布推定の前提として真の確率密度関数の形式が未知であるから、確率密度関数モジュールプログラムの実行により、カーネル法と呼ばれるノンパラメトリック（関数形式を指定しない）確率密度推定が行われる。推定対象の元データは、パルス状信号から得られた、例えば、波高ｈ・時間幅Δｔ・出現数等を含むパルス出現分布データである。元の計測データ分布の各データを計測誤差不確定性を導入したガウス分布で表し、各ガウス分布の重ね合わせにより確率密度関数が得られる。確率密度関数モジュールプログラムの実行により確率密度推定処理を行い、元データを該元データに基づいた未知の複雑な確率密度関数（例えば、特徴量のパルス波高・パルス幅・出現確率）で表すことができる。

図３３は、マイクロ・ナノポアデバイス８を用いて、３種の粒子３３ａ、３３ｂ、３３ｃが貫通孔１２を通過して得られる検出信号の波形例と、特徴量に基づいて得られる確率密度関数の導出例を示す。同図（３３Ａ）は、マイクロ・ナノポアデバイス８を用いた粒子検出装置を模式的に示す。同図（３３Ｂ）〜（３３Ｄ）は、各検出信号の波形データを示す。同図（３３Ｅ）〜（３３Ｇ）は、各波形データから得られた確率密度関数の３次元分布図を示す。（３３Ｅ）〜（３３Ｇ）における、ｘ軸、ｙ軸、ｚ軸は、それぞれ、特徴量のパルス波高、パルス幅および確率密度推定により得られた確率密度を示す。

本実施形態においては、上記のように、ノンパラメトリックな密度関数の推定法の一つであるカーネル法に基づいて確率密度推定処理が行われる。カーネル法は、１つのデータ点にある関数（カーネル関数）を当てはめ、これを全てのデータ点について行い、配置された関数を重ね合わせる推定法であり、滑らかな推定値を得るに適する。

確率密度関数モジュールプログラムの実行により、計測電流波形のパルス波高、パルス幅等のデータから多変数多次元確率密度とみなして２次元以上に拡張して加重の最適推定を行い粒子種別個数分布の推定処理が行われる。加重の最適推定にはＨａｓｓｅｌｂａｌｄ反復法に基づいて実行されるＥＭアルゴリズムソフトウエアが使用される。ＥＭアルゴリズムは、ＰＣ１にあらかじめインストールされている。粒子種別個数分布の推定処理により得られた粒子種別個数分布結果は、表示手段７に粒子種別に対する出現頻度（粒子個数）のヒストグラムで表示出力可能になっている。

本発明に係る特徴量は、パルス状信号由来のパラメータとして、該パルス状信号の波形の局所的特徴を示す第１類型に属するものと、該パルス状信号の波形の全体的特徴を示す第２類型に属するもののいずれかである。これらのうちの１または２以上の特徴量を使用して個数分析を行うことによって、粒子種等の分析物種別に応じた個数ないし個数分布を高精度に分析することができる。

図１１は、大腸菌２２と枯草菌２３の２種の粒子が電解質溶液２４中に混在する様子を模式的に示した貫通孔１２周辺の拡大図である。

＜特徴量について＞
図４は、実施例の大腸菌と枯草菌につき実測した粒子通過によるパルス波形例を示す。図４の（４−１）〜（４−９）は、大腸菌の実測パルス波形例（９種類）を示し、（４−１０）〜（４−１８）は、枯草菌の実測パルス波形例（９種類）を示す。両者を外観で比較すると、両者間に波高や波長には差異はあまりないが、ピーク位置や波形尖度等の粒子通過パルス波形形態の属性に顕著な相違がみられる。例えば、大腸菌の場合、ピークが時間経過に伴い前倒し傾向にあり、全体的に波形が尖っている（波形尖度が大きい）。枯草菌の場合、ピークが時間経過に伴い後倒し傾向にあり、波形尖度が小さい。

本発明者らは、上記の粒子通過パルス波形形態の属性の違いに基づいて、確率分布作成のベースに用いる特徴量をパルス波形データから粒子種（大腸菌と枯草菌）別に抽出することに着目した。

図５は、本発明に係る各種特徴量を説明するためのパルス波形図である。図５において、横軸は時間、縦軸はパルス波高を示す。

第１類型の特徴量は、
所定の時間幅内における波形の波高値、
パルス波長ｔ_a、
パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
該波形の鋭さ（ピーク波形の広がり）を表す尖度、
パルス開始からパルスピークに至る傾きを表す俯角、
波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかである。

図５の５ａ〜５ｄは、それぞれ、パルス波長、波高値、ピーク位置比、尖度を示す。図
５のＢＬは、パルス波形データから抽出（後述のＢＬ抽出処理参照）した基準ライン（以下、ベースラインという。）を示す。これら４種類のパルス特徴量は、図５に基づいて示すと以下の（１）〜（４）で定義される。

（１）波長（パルス幅）Δｔ： Δｔ＝ｔ_e−ｔ_s（ｔ_sはパルス波形の開始時間、ｔ_eはパルス波形の終了時間、Δｔ＝ｔ_a）
（２）波高｜ｈ｜： ｈ＝ｘ_p−ｘ_o（ＢＬのｘ_oを基準にしてパルスピークＰＰのｘ_pまでのパルス波形の高さ）
（３）ピーク位置比ｒ： ｒ＝（ｔ_p−ｔ_s）／（ｔ_e−ｔ_s）（パルス波長（＝Δｔ）と、パルス開始からパルスピークｐｐに至るまでの時間ｔ_b（＝ｔ_p−ｔ_s）との比）
（４）ピーク尖度κ： 波高｜ｈ｜＝１、ｔ_s＝０、ｔ_e＝１となるように正規化し、パルスピークＰＰから波高３０％の水平線と交差する時刻の時刻集合[Ｔ]＝[[ｔ_i]｜ｉ＝１，・・・，ｍ]を収集して、下記数１に示すように、時刻集合[Ｔ]のデータの分散をパルス波形広がりとしてκが求められる。

図３４は、俯角、面積および面積比の特徴量を説明するためのパルス波形図である。図３４において横軸は時間、縦軸はパルス波高を示す。これら３種類のパルス特徴量は、図３４に基づいて示すと以下の（５）、（６）、（７）で定義される。

（５）俯角θは、（３４Ａ）に示すように、パルス開始からパルスピークに至る傾きであり、下記数２により定義される。

（６）面積ｍは、下記数３に示すように、単位べクトル[ｕ]と波高ベクトル[ｐ]との内積による面積[ｍ]で定義される。なお、以下の説明において、変数Ａのベクトル表記は[Ａ]で示される。例えば、（３４Ｂ）の１０分割例に示すように、面積ｍは、一つの波形を所定の時間毎に１０分割したときの時間区分面積ｈ_i（幅ｈ_x、高さｈ_yとしたとき、ｈ_i＝ｈ_x×ｈ_y、i＝１〜１０）の総和を表す面積である。

ここで、特徴量計算の準備として、以下に定義するｄ次元波高ベクトル[ｐ]（＝（ｈ₁，ｈ₂，・・・，ｈ_d））をあらかじめ計算して求めておく必要がある。

図３５は、波高ベクトルの取得の仕方を説明するための図である。

（３５Ａ）に示すように、一つの波形データにつき、波長をｄ等分してｄ個のデータグループの分化が行われる。ついで、（３５Ｂ）に示すように、各グループ（各分割区間）ごとに波高の値を平均化して、例えば、１０等分するときには平均値Ａ１〜Ａ１０が求められる。この平均化には、波高値を規格化しない場合と、波高値を規格化する場合とを含むことができる。数３で表記した面積[ｍ]は、規格化しない場合を示す。このようにして求めた平均値を成分とするｄ次元ベクトルが「波高ベクトル」と定義される。

図３６は、ｄ次元の波高ベクトルとデータサンプリングと関係を説明するための図である。

（３６Ａ）に示すように、パルスデータの取得に関わるサンプリングレートが大きい場合、パルス部分におけるステップ数（データ数）Ｔがベクトルの次元数ｄを上回るので、上述の取得手順により各区分の平均値を成分した波高ベクトルを得ることができる。一方、サンプリングレートを下げていくと、パルス部分におけるステップ数Ｔがベクトルの次元数ｄ（＞Ｔ）を下回る事態が生ずる。Ｔ＜ｄの場合、上述の取得手順により各区分の平均値を取得できないので、３次スプライン補間によってｄ次元の波高ベクトルを取得することができる。

特徴量抽出プログラムには、波高ベクトルデータを取得するための波高ベクトル取得プログラムが含まれている。波高ベクトル取得プログラムの実行により、パルスステップ数Ｔがベクトルの次元数ｄを上回るか（Ｔ＞ｄ）同等の場合（Ｔ＝ｄ）、時間方向にｄ等分した各区分の平均値を求め、該平均値を成分とするｄ次元波高ベクトルを取得し、パルスステップ数Ｔがベクトルの次元数ｄを下回る場合（Ｔ＜ｄ）、３次スプライン補間を実行してｄ次元波高ベクトルを取得するようになっている。すなわち、３次スプライン補間法を用いた補間処理を行うことにより、パルスステップ数が少ない場合にもベクトルの次元数を一定にすることができる。

（７）面積比ｒ_mは、（３４Ｂ）に図示した時間区分面積ｈ_iをパルス開始からパルスピークに至るまでの区間での和の、全波形面積に対する面積比で定義される。下記数４は、面積比ｒ_mを示す。

第１類型の特徴量は、パルス波高、パルス波長、パルス面積等のパルス状信号の波形に一義的に由来し、局所的特徴を示す特徴量である。第２類型の特徴量は、第１類型の局所的特徴に対し全体的特徴を示す特徴量である。

第２類型の特徴量は、
パルス開始時点を中心にして時間区分面積を質量に、かつ該中心から時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかである。

図３７は、時間（波長）および波幅に関する慣性モーメントの特徴量を説明するためのパルス波形図である。図３７において横軸は時間、縦軸はパルス波高を示す。これらのパルス特徴量は、図３７に基づいて示すと、以下の（８）、（９）、（１０）、（１１）、（１２）、（１３）で定義される。

（８）時間慣性モーメントは、（３４Ｂ）と同様に一つの波形を所定の時間毎にｉ次元で等分割したときの時間区分面積ｈ_iを質量に、かつ該中心から時間区分面積ｈ_iに至る時間を回転半径に擬制したときに定まる特徴量である。すなわち、時間慣性モーメントの特徴量は、下記数５に示すように、ベクトル[ｖ]と波高ベクトル[ｐ]との内積による[Ｉ]で定義される。ここで、ベクトルの次元をｎとしたとき、[ｖ]＝（１²，２²，３²，・・・ｎ²）および[ｐ]＝（ｈ₁，ｈ₂，・・・，ｈ_d）である。例えば、時間慣性モーメントは、（３７Ａ）の１０分割例に示すように、（３４Ｂ）と同様に一つの波形を所定の時間毎に１０分割したとき、時間区分面積ｈ_i（幅ｈ_x、高さｈ_yとしたとき、ｈ_i＝ｈ_x×ｈ_y、i＝１〜１０）を質量に、かつ該中心から時間区分面積ｈ_iに至る時間を回転半径に擬制したときに定まる特徴量であり、（６）の面積ｍと同様に、波高ベクトルにより求めることができる。

（９）規格化された時間慣性モーメントは、（８）で示した時間区面積を作成した波形に対して波高が基準値の「１」になるように波高方向に規格化した波形を用いて、（８）と同様にして作成した波高ベクトルｈ_iにより数５で定義される特徴量である。

（１０）波幅平均値慣性モーメントは、（３７Ｂ）の１０分割例に示すように、（８）と同様に一つの波形を波高方向にｉ次元で等分割し、パルスピーク前後それぞれにおいて各分割単位(分割領域ｗ_i）毎に時刻値の平均値を算出し、分割領域ｗ_iの同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布ｈ_i（ベクトルの次元数をｎとしたとき、ｉ＝１〜ｎ）と擬制して波形裾野の時間軸Ａｔを回転中心にした場合の慣性モーメントとして定義される特徴量である。定義式は、数５と同じであり、（１０）の特徴量は、ベクトル[ｖ]と質量分布ｈ_iとの内積により求めることができる。

（１１）規格化された波幅平均値慣性モーメントは、（１０）で示した分割領域ｗ_iを作成した波形に対して波長が基準値の「１」になるように波長方向に規格化した波形を用いて、（１０）と同様にして作成した質量分布ｈ_iにより数５で定義される特徴量である。

（１２）波幅分散慣性モーメントは、波幅平均値慣性モーメントと同様に、
一つの波形を波高方向にｉ次元で等分割し、パルスピーク前後それぞれにおいて各分割単位(分割領域ｗ_i）毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布ｈ_i（ベクトルの次元数をｎとしたとき、ｉ＝１〜ｎ）と擬制して波形裾野の時間軸Ａｔを回転中心にした場合の慣性モーメントとして定義される特徴量であり、波幅平均値慣性モーメントと同様に、数５で定義される。

（１３）規格化された波幅分散慣性モーメントは、（１２）で示した分割領域ｗ_iを作成した波形に対して波長が基準値の「１」になるように波長方向に規格化した波形を用いて、（１２）と同様にして作成した質量分布ｈ_iにより数５で定義される特徴量である。

波幅平均値慣性モーメントおよび波幅分散慣性モーメントは、上記のように、数５で定義される特徴量であり、該定義中のベクトル[ｐ]は、波幅平均値慣性モーメントの場合、時刻値の平均値の差のベクトルであり、波幅分散慣性モーメントの場合、時刻値の分散ベクトルである。以下の説明で、（１０）〜（１３）の波幅に関する慣性モーメントにおけるベクトル[ｐ]を[ｐ_w]と表す。

（１０）〜（１３）の波幅に関する慣性モーメントのデータ作成演算には、図３６に示した波高ベクトルの縦横軸を交換した波幅ベクトル[ｐ_w]（＝[ｐ₁，ｐ₂，・・・，ｐ_dw]）を用いて行われる。波幅ベクトルは、（１０）〜（１３）の特徴量の定義で示した平均値の差ベクトルまたは分散ベクトルである。波幅ベクトルを密度分布として捉えることによって、（１０）と（１１）の波幅平均値慣性モーメントおよび（１２）と（１３）の波幅分散慣性モーメントを求めることができる。波幅ベクトルは、パルス波形データを波高
方向にｄ_w等分し、各区分毎に求めた波高値の平均値の差または分散を成分にもつｄ_w次元ベクトルである。（３７Ｂ）の場合、波幅ベクトルの次元は１０次元である。（３７Ｂ）に示した時間軸Ａｔは、ベースラインＢＬと異なり、波幅ベクトルから得られるパルス裾野周りの回転軸ラインである。

図３８は、ｄ_w次元の波幅ベクトルとデータサンプリングとの関係を説明するための図である。

特徴量抽出プログラムには、以下のｄ_w次元の波幅ベクトルの作成演算処理により波幅ベクトルを取得するための波幅ベクトル取得プログラムが含まれている。

パルス波形データは、波高方向に様々な間隔で分布しているため、波高方向に分割した区画においてデータ点が存在しない不存在領域Ｂｄを１または２以上含む場合が生ずる。（３８Ａ）においては、不存在領域Ｂｄの１例を矢印で示している。不存在領域Ｂｄは、データ間隔が粗くなってデータ点が存在せず、上記の数６で定義される波幅に関する慣性モーメントの成分を得ることができない。そこで、前述のパルス波形広がりの場合と同様に、パルスピークまでの波高をｄ_w等分したときの各波高の時刻集合[Ｔｋ]＝[[ｔ_i]｜ｉ＝１，・・・，ｍ]を収集することにより波幅ベクトルの成分が作成される。このとき、データ点が存在しない不存在領域Ｂｄにおいては、線形補間によって成分データの取得が行われる。該線形補間は、パルスピークの（１０ｋ＋５）％（ｋ＝０，１，２，３、・・・）の値をまたぐ連続した２つのデータに対して行われる。（３８Ｂ）は、データ点ｔ_iとｔ_i+1の間に生じた不存在領域Ｂｄの高さｋについての線形補間点ｔ_kの一例を示している。なお、波幅ベクトルの作成に際して、（３８Ｃ）に示すように、パルス波形データの裾野領域ＵＲに波高データの食い違いが生じている場合、パルスピークに近い方に揃えるように、パルスピークに遠い側の波高データＤｕは切り捨てられる。波幅ベクトル取得プログラムの実行処理には、不存在領域Ｂｄに対する線形補間処理と、波高データの食い違いに対する波高データＤｕの切り捨て処理が含まれている。

図３９は、波幅に関する慣性モーメントを波幅ベクトルにより取得する取得過程を説明するための図である。

（３９Ａ）は、波形を波高方向に１０等分した例であり、一つの波形３９ａに対して、上記の線形補間処理および切り捨て処理を行うことにより得られた、波幅ベクトルの分割領域３９ｂおよび回転軸ライン３９ｃを示す。

（３９Ｂ）に示すように、各分割単位ごとに、パルスピーク前後それぞれにおいて時刻値の平均値を算出し、分割領域の同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルの波幅ベクトルを取得することができる。この平均値の差ベクトルを質量分布と擬制して回転軸ライン３９ｃ（時間軸）を回転中心にした（１０）の波幅平均値慣性モーメントを作成することができる。また、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを取得することができる。この分散ベクトルを質量分布と擬制して回転軸ライン３９ｃ（時間軸）を回転中心にした（１２）の波幅分散慣性モーメントを作成することができる。

第２類型の特徴量には、上記の６種類に限らず、例えば、時刻値の平均値を算出して定まる平均値ベクトルおよびそれを規格化した特徴量を含めることができる。平均値ベクトルは、（３７Ｂ）の１０分割例に示すように、一つの波形を波高方向にｉ次元で等分割し、パルスピーク前後それぞれにおいて各分割単位(分割領域ｗ_i）毎に時刻値の平均値を算出し、分割領域ｗ_iの同一波高位置の平均値をベクトルの成分とする特徴量である。平均値ベクトルを規格化した特徴量は、該平均値ベクトルに対し波長が基準値になるように規
格化した場合の特徴量である。また、特徴量の作成に用いた波高ベクトルおよび波幅ベクトルのベクトル次元数は、分割数に拘泥する必要はなく任意に設定することができる。波高ベクトルおよび波幅ベクトルは、波長ないし波高の一方向に細分化した場合であるが、特徴量の作成には、複数方向に細分化したベクトルを使用することができる。

図４０は、複数の方向に分割した場合の特徴量作成用波形ベクトルの一例を説明するための図である。

（４０Ａ）は、一つの波形データをメッシュ状に分割したデータマップ４０ａを示す。データマップ４０ａは、波形データを横軸の時間軸方向にｄ_n分割し、縦軸の波高方向にｄ_w分割してマトリクス状にデータ点の数の分布状態を示す。（４０Ｂ）は、マトリクス状の区画（格子）の一部を拡大した分布状態を示す。（４０Ｂ）の分布状態において、１１×１３個の格子に０〜６個のデータ点数が分布している。このマトリクス分割によって、各格子内のデータ点数／総データ点数をｄ_n×ｄ_w次元ベクトルの成分とした波形ベクトルを、マトリクス配列のデータ群を走査状に並べ替えたベクトルに変換することにより波高ベクトルや波幅ベクトルに代えて特徴量の作成に使用することができる。

＜ベースラインの推定について＞
一般に細菌等は微細に異なる形態を有する微小物体である。例えば、平均的な大腸菌の場合、２〜４μｍの体長で外径が０．４〜０．７μｍである。平均的な枯草菌の場合、２〜３μｍの体長で外径が０．７〜０．８μｍである。さらに、大腸菌などには２０〜３０ｎｍの鞭毛が付随している。

検体粒子として細菌等を使用する場合には、パルス波形データから僅かの違いを見逃すと、個数判定精度の低下をもたらしてしまう。このため、特徴量を正確に算出して確率分布の推定基礎とするために、粒子通過パルス波高を正確に把握する必要があり、これには計測信号のベースラインの推定を行う必要がある。しかし、計測信号の元データのベースラインには、ノイズデータや微弱な計測電流による揺らぎが含まれているために、この揺らぎ成分等を除いたベースラインを確定してからパルス波高等を検出する必要がある。ベースラインの推定（以下、ＢＬの推定という。）は実用上、コンピュータによりオンラインで（即時的に）行うのが好ましい。

ＢＬの推定をコンピュータ上で行う手法として、離散的な誤差のある観測から時々刻々変化する量を推定するに好適なカルマンフィルターを使用すれば、これにより外乱（システム雑音や観測雑音）を取り除いてベースラインＢＬを推定することができる。

カルマンフィルターとは、離散的制御過程が図６の（６Ａ）に示す線形差分方程式により定義され、更新可能な状態ベクトル[ｘ]の時刻[ｔ]における値を推定する手法である。カルマンフィルターにおいては、状態ベクトル[ｘ]およびシステム制御入力[ｕ_t]の値は直接観測できないものとされる。

状態ベクトル[ｘ]は、図６の（６Ｂ）に示す観測モデルにより、間接的に推定されるものとする。システム制御入力[ｕ_t]については、その統計的変動幅[σ_u,t]のみをパラメータとして仮定する。

本実施形態における計測電流データ[Ｘ]はベクトルではなくスカラーであり、さらに各種行列もスカラーであり、[Ｆ]＝[Ｇ]＝[Ｈ]＝[１]とみなすことができる。したがって、時刻ｔの実際の電流値のベースラインレベル、時刻ｔで計測された電流、時刻ｔの観測ノイズをそれぞれ、[ｘ_t]、[ｙ_t]、[ν_t]とすると、[ｘ_t]および[ｙ_t]は図６の（６Ｃ）に示すように表される。[ｘ_t]、[ｕ_t]、[ν_t]は観測不可能な因子であり、[ｙ_t]は観測可能
な因子である。イオン電流検出部による計測周波数をｆ（Ｈｚ）とすると時刻データは１／ｆ（秒）刻みとなる。システム制御入力[ｕ_t]の影響は実際上非常に小さいものと仮定してベースラインの推定を行うことができる。

図７は上記の各因子を実際の計測電流データで示した図である。イオン電流検出部による実際の計測の際には、粒子が貫通孔１２に詰まったりして、ベースラインの歪みが生ずるが、計測時には歪みの発生時点で中断して、歪み原因を除去してから計測が行われるので、元のデータ集合には歪みのないベースラインを含むデータのみが収集されている。

カルマンフィルターによる推定は予測と更新の繰返しにより行われる。ベースラインの推定にもカルマンフィルターによる予測と更新を繰返して実行される。

図８は、カルマンフィルターにおける予測（８Ａ）と更新（８Ｂ）の繰返しの詳細を示す図である。図８において、ベクトル表記に付加した「ハット」記号は推定値を示している。添え字の「ｔ｜ｔ−１」は（ｔ−１）時点の値に基づく、ｔ時点の値の推定値であることを示している。

図９はＢＬ推定処理プログラムに基づくＢＬ推定処理を示す。ＢＬ推定処理にいては、ＢＬの推定と、それに基づくパルス波高値の抽出が行われる。

ＢＬ推定処理の実行に際しては、カルマンフィルターにおける予測と更新の処理に必要な調整因子の開始時刻ｍ、定数ｋ、αの値は推定対象のデータ属性に応じて適切な値にあらかじめ調整（チューニング）して決めておく必要がある。αの値はベースラインの推定値の分散を調整するための値である。ｋの値は図８に示したカルマンフィルターにおける更新Ａの実行回数に関係する値である（図９のステップＳ５７、Ｓ６２参照）。開始時刻ｍは計測サンプリングの１個分を１ステップとして計算されたステップ数分の時間データである。

図１０は該調整に使用したビーズモデルの波形図を示す。図２においては、粒子として細菌等と同程度の大きさの微小ビーズ玉を混入させた場合（ビーズモデル）の溶液状態を示している。図１０の（１０Ａ）はイオン電流検出部によってサンプリング周波数９０００００Ｈｚで取得した波形データである。（１０Ａ）に示すビーズモデルの波形はなだらかに減衰していく波形を示している。（１０Ａ）の右端部分に激しい落ち込みが生じており、それを拡大して（１０Ｂ）に示している。

ビーズモデルの波形から（１０Ｂ）に示すベースラインの段差部分（１０Ｃ）が検出された場合に、その直前期間が初期値計算期間となる。例えば、ｍ＝１０００００とした場合、当該初期値計算期間を除いた期間において有意性のあるパルスが１１〜１２個目視で確認することができる。

図１２は、調整因子のｍ、ｋ、αの組合せに応じてビーズモデルの波形から拾われたパルスの数を示す表である。

図１２の（１２Ａ）は、ｍ＝１００００の場合のｋ値（１０、３０、５０、７０、９０）、α値（２、３、４、６）の組合せによるパルス数を示す。同図（１２Ｂ）は、ｍ＝５００００の場合のｋ値（１０、３０、５０、７０、９０）、α値（２、３、４、６）の組合せによるパルス数を示す。（１２Ｃ）は、ｍ＝１０００００の場合のｋ値（１０、３０、５０、７０、９０）、α値（２、３、４、６）の組合せによるパルス数を示す。

図１２の３種のシミュレーション結果を比較すると、（１２Ａ）と（１２Ｂ）の場合、
計測されるべきパルス数は１２になり、（１２Ｃ）では１１となっている。したがって、実施例では、パルス数の最大値の最も小さい（１２Ｃ）を採用して、ｍ＝１０００００、ｋ＝５０、α＝６のチューニング設定を行っている。これらのチューニング設定データはあらかじめＲＡＭ２３の設定エリアに記憶、設定されている。

図９のＢＬ推定処理は上記チューニング設定下、図８に示したカルマンフィルターによるＢＬ推定で行われる。まず、ステップＳ５１にて、時刻ｍにおけるカルマンフィルターの初期値がＲＡＭ２３のワークエリアに設定される。このとき、データファイル記憶部５に格納したパルス波形データはＲＡＭ２３のワークエリアに読み込まれる。ついで、時刻（ｍ＋１）におけるカルマンフィルターの予測と更新（図８のＡおよびＢ）が実行される（ステップＳ５２）。予測と更新においては、図８に示したカルマンフィルターの各演算が実行され、ＲＡＭ２３に記憶される。以降、所定の単位時間毎に予測と更新（ＡおよびＢ）が繰り返し実行され、時刻ｔにおけるカルマンフィルターの予測と更新Ａが行われたとき、下記数６の条件が満たされたか否かが判断される（ステップＳ５３、Ｓ５４）。単位時間は元データのサンプリング周波数により定まる値であり、あらかじめＲＡＭ２３にセットされている。

数６の条件が満たされない場合、時刻ｔにおけるカルマンフィルターの更新Ｂが実行され、単位時間経過したデータ毎にステップＳ５３〜Ｓ５５の処理が繰り返される。上記数６の条件が満たされた場合、その回数値が１回ごとにＲＡＭ２３のカウントエリアに累積記憶される（ステップＳ５４、Ｓ５６）。ついで、該カウント値に基づき、数６の条件が時刻ｓを起点としてｋ回連続して満たされたか否かが判断される（ステップＳ５７）。ｋ回連続していない場合はステップＳ５５に進み、更新Ｂが行われる。

ｋ回連続した場合はステップＳ５８に進み、ＢＬ確定のためのホールド必要期間が開始したと判定される。このとき、ホールド必要期間のホールド開始時刻をｓとしてＲＡＭ２３に記憶されるとともに、時刻（ｓ＋１）〜時刻（ｓ＋ｋ−１）の間のカルマンフィルターの演算結果は記憶されずに捨てられる。

ホールド必要期間の開始により、時刻ｔにおけるパルスの落ち込み最大値がＲＡＭ２３に更新可能に記憶される（ステップＳ５９）。ついで、ステップＳ５４と同様に、ホールド必要期間における、下記数７の条件が満たされているか否かの判断が行われる（ステップＳ６０）。

上記数７の条件が満たされない場合、パルスの落ち込み最大値の更新が行われる（ステップＳ５９、Ｓ６０）。数７の条件が満たされた場合、その回数値が１回ごとにＲＡＭ２３のカウントエリアに累積記憶される（ステップＳ６０、Ｓ６１）。ついで、該カウント値に基づき、数７の条件が時刻ｓ２を起点としてｋ回連続して満たされたか否かが判断される（ステップＳ６２）。ｋ回連続していない場合はステップＳ５９に戻る。

ｋ回連続した場合はステップＳ６３に進み、このとき更新記憶されたパルスの落ち込み最大値がパルス波高値の推定値としてＲＡＭ２３に記憶される。パルス波高値の推定値はパルス開始時刻およびパルス終了時刻のデータとともに記憶される。パルス波高値の推定を終えると、ホールド必要期間は終了と判定される。この終了によりホールド必要期間のホールド終了時刻がｓ２としてＲＡＭ２３に記憶される（ステップＳ６４）。次に、ステップＳ６５に進み、時刻ｓの値はカルマンフィルターの演算処理の再開時の初期値として、時刻ｓ２〜時刻（ｓ＋ｋ−１）の期間について遡及してカルマンフィルターの演算が実行される。ステップＳ６５の後は、全パルス波形データのＢＬ推定処理を行ったか否かが判断されて（ステップＳ６６）、全パルス波形データの推定完了で終了し、残データがあるときはステップＳ５３に移る。

＜特徴量抽出について＞
図１３は、特徴量抽出プログラムの実行処理内容の概要を示す。

特徴量抽出処理は、図９の上記ＢＬ推定処理の実行によってパルス波高値（波高｜ｈ｜）の抽出データがあることを条件に実行可能になる（ステップＳ４１）。パルス波高値の抽出データがある場合、前述の波高ベクトル取得プログラムおよび波幅ベクトル取得プログラムが実行されて、各種ベクトルのデータ作成演算が実行される（ステップＳ４２）。波高ベクトルおよび波幅ベクトルの全てのデータ取得を終えると、該ベクトルデータが保存される（ステップＳ４３、Ｓ４４）。ついで、各種の特徴量の抽出処理が実行される（ステップＳ４５）。波高ベクトルおよび波幅ベクトルのデータ取得に際しては、３次スプライン補間法を用いた補間処理、線形補間処理および切り捨て処理が随時行われる。

図４１は、特徴量の抽出処理（ステップＳ４５）の実行処理内容を示す。ステップＳ７１〜Ｓ８３は、それぞれ上記（１）〜（１３）で定義された第１類型および第２類型の特徴量の算出と、算出された特徴量の記憶、保存の処理を示す。

第１類型の特徴量は、ステップＳ７１〜Ｓ７６において算出される。波長（パルス幅）Δｔは、パルス波高値の抽出データ群に対し時系列的に順次算出されて記憶される（ステップＳ７１）。算出された特徴量は、ＲＡＭ４の特徴量記憶用メモリエリアに記憶される。パルス幅は、Δｔ（＝ｔ_e−ｔ_s;ｔ_sはパルス波形の開始時刻、ｔ_eはパルス波形の終了時刻）を演算して求められる。ピーク位置比ｒは、パルス波高値の抽出データ群に対し時系列的に順次算出されて記憶される（ステップＳ７２）。ピーク位置比ｒは、ｒ＝（ｔ_p−ｔ_s）／（ｔ_e−ｔ_s）（パルス幅Δｔと、パルス開始からパルスピークｐｐに至るまで
の時間（＝ｔ_p−ｔ_s）との比）を演算して求められる。

ピーク尖度κは、パルス波高値の抽出データ群に対し時系列的に順次算出されて記憶される（ステップＳ７３）。パルス波高値｜ｈ｜＝１、ｔ_s＝０、ｔ_e＝１となるように正規化し、パルスピークＰＰから波高３０％の水平線と交差する時刻の時刻集合Ｔ＝[[ｔ_i]｜ｉ＝１，・・・，ｍ]を収集して、時刻集合Ｔのデータの分散を演算してパルス波形広がりとしてκが求められる。

俯角θは、パルス開始からパルスピークまでの時刻と波高のデータと、前掲の数２の演算とに基づいて求められる（ステップＳ７４）。面積ｍは、波高ベクトルのデータにより求められ、時間区分面積ｈ_iを区分数（あらかじめ設定された区分数：１０）に応じて求め、それらの総和を求めることにより算出、記憶される（ステップＳ７５）。面積比ｒ_mは、全波形面積と、時間区分面積ｈ_iをパルス開始からパルスピークに至るまでの区間での部分和とをそれぞれ求め、部分和の全波形面積に対する面積比を算出して記憶される（ステップＳ７６）。

第２類型の特徴量は、ステップＳ７７〜Ｓ８２において算出される。時間慣性モーメントは、波高ベクトルのデータにより求められ、区分数に応じて求めた時間区分面積ｈ_iと、前掲の数５の演算とに基づいて算出され、記憶される（ステップＳ７７）。（９）の規格化された時間慣性モーメントは、ステップＳ７７において得られた時間慣性モーメントに対し波高が基準値の「１」になるように波高方向に規格化処理（波高ベクトルと規格化ベクトルの内積）した規格化データとして記憶される（ステップＳ７８）。波幅平均値慣性モーメントは、ステップＳ４２〜Ｓ４４で求めた波幅ベクトル（平均値の差ベクトル）のデータから、パルスピーク前後それぞれにおいて分割単位(あらかじめ設定された分割数：１０）毎に算出した時刻値の平均値の差と、前掲の数６の演算とに基づいて算出され、記憶される（ステップＳ７９）。（１１）の規格化された波幅平均値慣性モーメントは、ステップＳ７９において得られた波幅平均値慣性モーメントに対し波長が基準値「１」になるように波長方向に規格化処理（平均値の差ベクトルと規格化ベクトルの内積）した規格化データとして記憶される（ステップＳ８０）。波幅分散慣性モーメントは、波幅ベクトル（分散ベクトル）のデータから、分割単位毎に算出した時刻値の分散と、前掲の数６の演算とに基づいて算出され、記憶される（ステップＳ８１）。（１３）の規格化された波幅分散慣性モーメントは、ステップＳ８１において得られた波幅分散慣性モーメントに対し波長が基準値「１」になるように波長方向に規格化処理（分散ベクトルと規格化ベクトルの内積）した規格化データとして記憶される（ステップＳ８２）。

全データからの特徴量の抽出を終了すると、各データのファイル保存が行われ、別のデータ群があるか否かが判断される（ステップＳ８３、Ｓ８４）。引き続き別ファイルのデータ群があれば、上記処理（ステップＳ７１〜Ｓ８２）が繰返し実行可能になっている。処理すべきデータがなくなれば特徴量の抽出処理は終了する（ステップＳ８５）。上記の抽出処理においては、第１類型および第２類型の特徴量を全て求めるようにしているが、入力手段６の指定入力により所望の特徴量を指定可能にし、該指定による特徴量だけを抽出可能にすることができる。

図１４は、粒子種分布推定プログラムに基づいて実行される粒子種推定処理を示す。
＜確率密度関数の推定について＞
同種の粒子であっても計測されるパルス波形が一定とは限らないので、粒子種分布推定のための準備として、テストデータから予め粒子種別のパルス波形の確率密度関数の推定が行われる。確率密度関数の推定により導出される確率密度関数によって各パルスの出現確率を表すことができる。

図１５の（１５Ｂ）は大腸菌と枯草菌の粒子種において、パルス波形の特徴量としてパルス幅とパルス波高が用いて得られたパルス波形に対する確率密度関数のイメージ図であり、図中の濃淡によりパルスの出現確率を表している。図１５の（１５Ａ）は１つの波形データに関する第１類型の特徴量の一部を示す。

パルス幅Δｔとパルス波高ｈの真の密度関数は未知であるので、ノンパラメトリックな確率密度関数の推定を行う必要がある。本実施形態では、カーネル関数としてガウス関数を採用したカーネル密度推定を用いている。

カーネル密度推定とは、計測データにカーネル関数で与えられた確率密度分布を想定し、それらの分布を重ね合わせた分布を確率密度関数とみなす手法である。カーネル関数としてガウス関数を用いた場合は、各データに対して正規分布を想定し、それらを重ね合わせた分布を確率密度関数とみなすことができる。

図１６は、大腸菌と枯草菌の粒子種の個々より得られた確率密度分布の重ね合わせのイメージ図である。同図（１６Ｃ）は、パルス幅Δｔとパルス波高ｈの特徴量データ（１６Ａ）から、各粒子につき求めた確率密度分布（１６Ｂ）を重ね合わせた状態を示している。

入力データ[ｘ]に対する確率密度関数[ｐ（ｘ）]は教師データ数[Ｎ]、教師データ[μ_i]、分散共分散行列[Σ]を用いて、下記数８で表される。

さらに、確率密度関数[ｐ（ｘ）]は下記数９に示すように、各次元のガウス関数の積で表すことができる。

数９からわかるように、各パルス属性が正規分布に従う独立な確率変数であると仮定していることに相当し、これは３次元以上にも同様に拡張可能である。したがって、本実施形態においては２種以上の粒子種個数の分析が可能である。

確率密度関数モジュールプログラムは、２種の特徴量に対する確率密度関数を演算して求める機能を有する。すなわち、２つの特徴量[（β，γ）]による推定対象データを使う場合、カーネル関数としてガウス関数を採用したカーネル密度推定における確率密度関数[ｐ（β，γ）]は下記数１０で表される。

数１０に基づいて確率密度関数モジュールプログラムにより実行される確率密度関数推定処理は後述の図２０により詳述するように、２つの特徴量における確率密度関数の推定処理を行う。

図１７はｋ個の粒子種別の粒子総数と、粒子種別の出現確率と、データ全体の出現頻度の期待値との関係を示すイメージ図である。同図（１７Ａ）はデータ全体の出現頻度を示す。同図（１７−１）〜（１７−ｋ）は粒子種別の出現頻度を示す。パルス[ｘ]が計測される出現頻度の期待値は、粒子種別の確率密度関数に従ってパルス[ｘ]が計測される出現頻度の期待値の和となる。図１７に示すように、粒子種別の粒子総数[ｎ_i]と、粒子種別出現確率[ｐ_i（ｘ）]とから粒子種別の期待値の和として下記数１１で表すことができる。

本実施形態においては、あらかじめ求められた粒子種別の確率密度関数の推定を行った確率密度関数データ（数９参照）が分析参照データとしてＲＡＭ２３に記憶されている。粒子種別個数分析は、数１１に基づいて、分析対象の全体データの出現頻度を各分析データから適合する粒子種別の個数を割り出すことにより行われる。個数分析は異なる粒子種のヒストグラム（粒子種に対する出現頻度（粒子数））を推定することにより行われる。

図１４の粒子種推定処理においては、データの編集によって特徴量によるデータファイルを作成するデータファイル作成処理（ステップＳ１）と、粒子数の推定処理（ステップＳ２）と、推定粒子種分布の算出処理（ヒストグラム作成処理）（ステップＳ３）が行われる。粒子数の推定処理においては、最尤法、ラグランジュ未定乗数法およびＨａｓｓｅｌｂｌａｄ反復法による推定手法が用いることができる。
＜最尤法（ｍａｘｉｕｍ ｌｉｋｅｌｉｈｏｏｄ ｅｓｔｉｍａｔｉｏｎ：統計学において与えられたデータからそれが従う確率分布の母数を点推定する方法である。）について
＞
今、実際のパルス推定結果として、データセット[Ｄ]＝[ｘ₁，ｘ₂，ｘ₃，・・・ｘ_N]が得られれているとする。推定されたｊ番目のパルス波高データが出現する尤度（尤もらしさ）は下記数１２で表わされる。

するとデータセットＤが出現する尤度は下記数１３で表される。

数１３の尤度を最大化する様な粒子種分布の値セット[ｎ]＝[ｎ₁，・・・，ｎ_k]^Tが最も尤もらしい粒子種分布である。
＜ラグランジュ未定乗数法(束縛条件のもとで最適化を行う解析学的方法であり、各束縛条件に対して未定乗数を用意し、これらを係数にする線形結合を新しい関数（未定乗数も新たな変数とする）として捉えることによって束縛問題を通常の極値問題として解く方法である。）について＞
データセットＤが出現する尤度を最大化することは、データセット[Ｄ]が出現する対数尤度を最大化することに等しい。下記数１４はラグランジュ未定乗数法の適否を調べるための対数尤度を導出する過程を示す。

数１４において、途中の係数１／Ｎ^Nは最終式では省略している。

ここで、粒径個数分布の値セットｎ＝[ｎ₁，・・・，ｎ_k]^Tには、「合計がＮである」という制約（下記数１５参照）がある。

したがって、最も尤もらしい粒子種分布を得るという命題は、制約付き対数尤度最大化の問題になるので、ラグランジュ未定乗数法により最適化を行うことが可能である。ラグランジュ未定乗数法により最適化を行う制約付き対数尤度最大化式を下記数１６で表すことができる。

数１６に示す制約付き対数尤度最大化式からは、図１８に示す数学的導出過程を経て下記数１７に示す[ｋ]個の連立方程式を導き出すことができる。

数１７に示す連立方程式を数値的に解くには、Ｈａｓｓｅｌｂｌａｄが提唱した反復法を用いて行うことができる。Ｈａｓｓｅｌｂｌａｄ反復法によれば、下記数１８の反復計算を行えばよい。この反復法の詳細は提唱論文（Ｈａｓｓｅｌｂｌａｄ Ｖ ．，１９６６，Ｅｓｔｉｍａｔｉｏｎ ｏｆ ｐａｒａｍｅｔｅｒｓ ｆｏｒ ａ ｍｉｘｔｕｒｅ
ｏｆ ｎｏｒｍａｌ ｄｉｓｔｒｉｂｕｔｉｏｎｓ． Ｔｅｃｈｎｏｍｅｒｉｃｓ，８，ｐｐ．４３１−４４４）に記述されている。

数１８の反復計算には、市販のＥＭアルゴリズムのソフトウエアを利用して行われる。ＥＭアルゴリズムは、命名の由来から明らかなように、確率分布のパラメータを、尤度関数を最大化することで計算する方法、つまり尤度関数である確率分布の期待値（Ｅｘｐｅｃｔａｔｉｏｎ）を最大化（Ｍａｘｉｍｉｚａｔｉｏｎ）することができるアルゴリズムである。ＥＭアルゴリズムによれば、求めたいパラメータの初期値を設定して、その値から尤度(期待値)を計算して、多くの場合、尤度関数の偏微分が０になる条件を使って、繰り返し計算で最大尤度のパラメータを計算することができる。ＥＭアルゴリズムを使用して行うＨａｓｓｅｌｂｌａｄ反復法の演算処理は、求めるパラメータの初期値を設定して、その値から尤度(期待値)を計算し、さらに尤度関数の偏微分が０になる条件を使って繰り返し計算を行って最大尤度のパラメータを計算する工程を有する。

＜粒子種推定処理について＞
図１４に示した粒子種推定処理において実行可能な、データファイル作成処理（ステップＳ１）、確率密度関数の推定処理（ステップＳ２）、粒子数の推定処理（ステップＳ３）および推定粒子種分布の算出処理（ステップＳ４）を以下に詳述する。

図１９はデータファイル作成プログラムにより実行されるデータファイル作成処理（ステップＳ１）を示す。

ＰＣ１の入力手段６を使用して、データファイルを作成するｋ個（実施例では２個）ずつの特徴量の指定操作を行うことができる。指定された特徴量の組合せ入力がＲＡＭ２３に設定される（ステップＳ３０）。特徴量の設定毎の特徴量データファイルのデータがＲＡＭ２３のワークエリアに読み込まれる（ステップＳ３１）。特徴量データファイルは、図９のＢＬ推定処理および図１３の特徴量推定処理で推定され抽出され、ファイル保存されている特徴量（パルス波高値等）データである。

個数推定に使う特徴量をｋ個指定することにより、Ｎ行ｋ列の行列データが作成される（ステップＳ３２）。作成された行列データは、粒子種分布推定用データファイルに出力され、指定特徴量別に保存される（ステップＳ３３）。指定特徴量に対するすべてのデータファイルの生成を終えると終了する（ステップＳ３４）。

図２０は、確率密度関数モジュールプログラムにより実行される確率密度関数の推定処理（ステップＳ２）を示す。確率密度関数推定処理は、数５に基づいて２つの特徴量における確率密度関数の推定処理が行われる。

データファイル作成処理（ステップＳ１）において作成された、確率密度関数推定対象のデータファイルのデータを読み込んで、Ｎ行２列の行列[Ｄ]が作成される（ステップＳ２０、Ｓ２１）。行列[Ｄ]の列ごとの下記数１９に示す分散

が算出される（ステップＳ２２）。ついで、下記数２０に示す分散パラメータが標準偏差係数ｃを用いて下記数２１に示すように設定される（ステップＳ２３）。

分散パラメータおよび行列[Ｄ]の各行を下記数２２に示す教師データとして代入されて確率密度関数が求められ、ＲＡＭ２３の所定エリアに記憶される（ステップＳ２４、Ｓ２５）。上記のステップＳ２０〜Ｓ２５の処理は全ての処理対象データからの確率密度関数の導出を行うまで行われる（ステップＳ２６）。

図２１は粒子数の推定処理（ステップＳ３）を示す。

まず、上記のステップＳ２０、Ｓ２１と同様に、データファイル作成処理において作成された、粒子数推定対象のデータファイルのデータを読み込んで、Ｎ行２列の行列[Ｄ]が作成される（ステップＳ１０、Ｓ１１）。行列[Ｄ]データに対して、Ｈａｓｓｅｌｂｌａｄ反復法による推定処理が実行される（ステップＳ１２）。

図２２は、ＥＭアルゴリズムによって実行される、Ｈａｓｓｅｌｂｌａｄ反復法による粒子数推定処理を示す。図２３はＥＭアルゴリズムによる処理手順を示す。

まず、初期値の設定（処理１９Ａ）を行った後、確率密度関数に基づく個数計算が順次実行される（処理１９Ｂ）（ステップＳ１２ａ、Ｓ１２ｂ）。個数計算の反復は（１９Ｃ）に示す収束条件を満たすまで実行される（ステップＳ１２ｃ）。ＥＭアルゴリズムの実行結果（粒子種ごとの推定個数データ）はＲＡＭ２３の所定エリアに格納される（ステップＳ１２ｄ）。

粒子数推定処理により得られた粒子種ごとの推定個数データは、ステップ４において、粒子種別の個数分布データに編集され、表示指定に応じて表示手段７にヒストグラム表示出力可能になる。図１４では省略しているが、本実施形態においては、分散図出力の指定を受けた場合には、特徴量データによる粒子種別の分散図を表示出力可能になっている。

図２４は本実施形態に係る粒子種個数分析装置により分析した結果の一例を示す。同図（２４Ａ）および（２５Ｂ）は分析対象の粒子種である大腸菌、枯草菌の顕微鏡拡大写真
である。（２４Ｃ）および（２５Ｄ）は特徴量としてパルス波高およびパルス尖度を注して粒子数推定処理の実行により得られた粒子種ごとの推定個数データのヒストグラム、分散図を示す。

＜特徴量による粒子種個数の分析精度の検証１について＞
本発明者らは、上記実施例の大腸菌と枯草菌の計測電流データを用いて下記の評価条件下で粒子種個数の分析性能の検証１を行った。

検証１の評価条件は以下の通りである。

（１）大腸菌と枯草菌の１０００ｋＨｚ実験測定データで評価を行う。

（２）特徴量には、波長Δｔ、波高ｈ、ピーク位置比ｒ、ピーク尖度ｋの４つの第１類型の特徴量を算出して用いる。

（３）各特徴量の組合せについて個数推定処理を実施する。

（４）大腸菌と枯草菌の実測データをランダムに学習用とテスト用に分けて推定評価する。この推定評価を１０回反復して実施し、それらの平均精度と標準偏差を算出する。この場合、実際に近い精度を評価する交差検定法（ｃｒｏｓｓ ｖａｌｉｄａｔｉｏｎ）により行う。

（５）検証粒子（大腸菌と枯草菌）の実測データの一部を個別に個数分析し、残りを所定の混合比δによりランダムに混合して検証用として、個数分析結果を比較する。ランダムデータ混合用データ混合プログラムをＲＯＭ３に格納しＰＣ１を利用してデータのランダム混合を実行して、そのランダム混合したデータに対する個数推定を行う。すなわち、図１９のステップＳ３２における行列データには、データ混合プログラムにより作成されたＮ行ｋ列のランダム置換行列データを使用する。混合比δには、大腸菌の混合率として１０、２０、３０、３５、４０、４５、５０％の７種類を使用する。ＢＬ推定用のパラメータ(調整因子）ｍ、ｋ、αの値はそれぞれ、１０００００、４００、６を使用し、確率密度関数の推定用の標準偏差係数ｃには０．１を設定する。粒子種個数推定時の収束条件αは０．１に設定する。なお、評価に使用した上記調整因子の値には、図１２で示したシミュレーション例と同様にしてより厳密な調整を行って得られた値が使用されている。

図２５の（２５Ａ）および（２５Ｂ）は、特徴量としてパルス波長、波高を使用した検証例と、特徴量としてパルス波長、ピーク位置比を使用した検証例の各推定結果データを示す。

本検証により得られた全パルスの数は大腸菌で１４６個、枯草菌で４０５個であった。

図２６の（２６Ａ）および（２６Ｂ）は特徴量としてピーク付近波形の広がり、パルス波長を使用した検証例と、特徴量としてピーク付近波形の広がり、波高を使用した検証例の各推定結果データを示す。

図２６の（２６Ａ）および（２６Ｂ）は特徴量としてピーク付近波形の広がり、パルス波長を使用した検証例と、特徴量としてピーク付近波形の広がり、波高を使用した検証例の各推定結果データを示す。

粒子種別個数の評価は、図２７の（２７Ｂ）に示す数式で表わされる「重み付き平均相対誤差」により行うことができる。「重み付き平均相対誤差」は各粒径の相対誤差にその
粒径の真の個数割合を掛けたものを、全粒径について足した数値である。

図２７の（２７Ａ）は特徴量として尖度と、パルス波高を使用した場合における個数推結果を示す。

図２８の（２８Ａ）および（２８Ｂ）は特徴量としてパルス波長、パルス波高を使用した場合における各混合比δ別の個数推定結果と、特徴量としてパルス波長、ピーク位置比を使用した場合における各混合比δ別の個数推定結果を示す。

図２９の（２９Ａ）〜（２９Ｄ）は大腸菌と枯草菌の混合比をそれぞれ、１：１０、２：１０、３：１０、３５：１００とした場合における各個数推定結果を示すヒストグラムである。

図３０の（３０Ａ）〜（３０Ｃ）は、大腸菌と枯草菌の混合比をそれぞれ、４：１０、４５：１００、１：２とした場合における各個数推定結果を示すヒストグラムである。

図３１の（３１Ａ）および（３１Ｂ）は、特徴量としてパルス波長、パルス波高を使用した場合における各粒子の散布状態を合成した図である。

図３２の（３２Ａ）、（３２Ｂ）および（３２Ｃ）は、特徴量としてピーク付近波形の広がり、パルス波長を使用した場合、特徴量としてピーク付近波形の広がり、ピーク位置比を使用した場合、ピーク付近波形の広がり、パルス波高を使用した場合における各粒子の散布状態を合成した図である。

上記の性能評価実験から、以下の評価結果が得られた。

（１）図３１および図３２のデータ散布図では、４つの特徴量に関して、大腸菌と枯草菌の特徴は大きく重なるが、明らかな相違があることが認められる。

（２）図２７の（２７Ａ）等に示す種別個数分布の推定結果からは、この評価検証の特徴量の中でパルス波高とピーク尖度の特徴量を組み合わせた場合が最も精度がよく、重み付き平均相対誤差の評価で４〜１２％の分析精度を得ることができる。上記実施形態においては４種類全部の特徴量を抽出しているが、上記検証結果を踏まえて一部の特徴量（例えば、パルス波高とピーク尖度）だけを抽出して個数分析するようにしてもよい。
＜特徴量による粒子種個数の分析精度の検証２について＞
本発明者らは、上記実施例の大腸菌と枯草菌の計測電流データを用いて、検証１とは別の粒子種個数の分析性能の検証２を行った。検証２においては、検証１とは異なり、第１類型および第２類型の特徴量（（１）〜（１３）の１３種類）を算出して用い、これらの組合せに係る特徴量とサンプリングデータ数との関連性および各組合せの分析性能を検証した。

図４２の（４２Ａ）および（４２Ｂ）は、それぞれ、全データのうち、１ＭＨｚ、５００ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４３の（４３Ａ）および（４３Ｂ）は、それぞれ、全データのうち、２５０ｋＨｚ、１２５ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４４の（４４Ａ）および（４４Ｂ）は、それぞれ、全データのうち、６３ｋＨｚ、３２ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４５の（４５Ａ）および（４５Ｂ）は、それぞれ、全データのうち、１６ｋＨｚ、８ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４６は、４ｋＨｚでサンプリングしたときの各特徴量組合せに関する推定評価結果を示す。これらの表中の各組合せご
との推定評価結果は、検証１の（４）と同様に交差検定法により得られた、上側に記載の平均精度と、下側に括弧書きで示した標準偏差を表す。表中の慣性Ｉ、慣性Ｉ（規格化）、慣性Ｉ ｗ、慣性Ｉ ｗｖ、慣性Ｉ ｗ（規格化）、慣性Ｉ ｗｖ（規格化）は、それぞれ、（８）の時間慣性モーメント、（９）の規格化された時間慣性モーメント、（１０）の波幅平均値慣性モーメント、（１２）の波幅分散慣性モーメント、（１１）の規格化された波幅平均値慣性モーメント、（１３）の規格化された波幅分散慣性モーメントの特徴量を示す。

図４７は、全サンプリングデータにおける各特徴量組合せに関する推定評価結果を示す。図４８は、全データのうち１ＭＨｚ〜１２５ｋＨｚでの高密度サンプリングしたときの各特徴量組合せに関する推定評価結果を示す。図４９は、全データのうち６３ｋＨｚ〜４ｋＨｚでの低密度サンプリングしたときの各特徴量組合せに関する推定評価結果を示す。

図５０は、全サンプリングデータを使用したとき（５０Ａ）および高密度にサンプリングしたとき（５０Ｂ）に高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフである。図５０における上位５種の特徴量の組合せは、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、ピーク位置比ｒ−慣性Ｉ、俯角θ−慣性Ｉ、慣性Ｉ−慣性Ｉ ｗ（規格化）である。

図５１は、低密度にサンプリングしたときに高い個数推定精度が得られる上位５種の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ａ）と、全サンプリングデータを使用したときの４種類の特徴量の組合せに関するサンプリング周波数−重み付き平均相対誤差（平均値）のグラフ（５１Ｂ）である。図５０および図５１の縦軸の値は、５０回の交差検定を行って得られた重み付き平均相対誤差の平均値である。５１Ａにおける上位５種の特徴量の組合せは、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、ピーク位置比ｒ−面積ｍ、俯角θ−面積ｍ、面積ｍ−慣性Ｉ ｗｖ（規格化）である。５１Ｂにおける４種類の特徴量の組合せは、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、尖度ｋ−波高｜ｈ｜、尖度ｋ−ピーク位置比ｒである。

検証２から得られた結果は、以下の通りである。
（Ｒ１）図４７および図５０に示すように、全サンプリングデータを使用したときには、上位５種の組合せ、すなわち、波長Δｔ−慣性Ｉ、波長Δｔ−面積ｍ、ピーク位置比ｒ−慣性Ｉ、俯角θ−慣性Ｉ、慣性Ｉ−慣性Ｉ ｗ（規格化）の特徴量の場合、高い個数推定精度を得ることができる。これらの特徴量の組合せによる個数推定精度（重み付き平均相対誤差）は、例えば、波長Δｔ−慣性Ｉで２５０〜１０００ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−面積ｍで１２５〜２５０ｋＨｚのサンプリング領域において約９〜１０％、、波長Δｔ−慣性Ｉで１６〜６３ｋＨｚのサンプリング領域において約１３〜１５％である。

（Ｒ２）図４８に示すように、全サンプリングデータより少ないが高密度なサンプリングデータを使用したときに高い個数推定精度が得られる特徴量は、上位５種の組合せで示せば、波長Δｔ−慣性Ｉ、波長Δｔ−面積ｍ、ピーク位置比ｒ−慣性Ｉ、慣性Ｉ−慣性Ｉ ｗ、俯角θ−慣性Ｉの５種である。これらの特徴量の組合せによる個数推定精度（重み付き平均相対誤差）は、例えば、波長Δｔ−慣性Ｉで２５０〜１０００ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−面積ｍで１２５〜２５０ｋＨｚのサンプリング領域において約９〜１０％、、波長Δｔ−慣性Ｉで１６〜６３ｋＨｚのサンプリング領域において約１３〜１５％である。

（Ｒ３）図４９に示すように、高密度サンプリングデータと比べてさらに少ない低密度サンプリングデータを使用した場合に高い個数推定精度が得られる特徴量は、上位５種の組
合せで示せば、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、俯角θ−面積ｍ、面積ｍ−慣性Ｉ
ｗｖ（規格化）、ピーク位置比ｒ−面積ｍの５種である。これらの特徴量の組合せによる個数推定精度（重み付き平均相対誤差）は、波長Δｔ−慣性Ｉで２５０〜１０００ｋＨｚのサンプリング領域において約９〜１０％、波長Δｔ−面積ｍで１２５〜２５０ｋＨｚのサンプリング領域において約９〜１０％、、波長Δｔ−慣性Ｉで１６〜６３ｋＨｚのサンプリング領域において約１３〜１６％である。

（Ｒ４）（Ｒ１）〜（Ｒ３）からわかるように、第１類型と第２類型の特徴量の組合せを使用しても高精度の個数推定を行うことができる。さらに、本発明に係る個数分析方法によれば、サンプリング数が十分に多くなくても、所定のサンプリング数が得られれば十分あるときと同程度の精度で個数分析を行うことができる。例えば、検証１で調べた尖度ｋとピーク位置比ｒとの組合せでは、１２％の最大誤差を生じていたが、例えば、波長Δｔ−慣性Ｉの特徴量による場合には、全データを使用しなくとも１ＭＨｚ〜１２５ｋＨｚでの高密度サンプリングデータを使用して、つまり部分的なデータであっても個数推定処理を約９％の高精度で行うことができる。したがって、本実施形態に係る個数分析装置は、定常的な個数分析にとどまらず、例えば、緊急性を要する検疫検査や医療現場において、菌類等の粒子有無や個数の判別に即応的実施に好適な検査ツールとして使用することができる。

＜個数分析処理時間の検証３について＞
個数推定には、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する所要計算時間がかかるので、この所要計算時間とサンプリング周波数との関係について特徴量の比較検討を検証３で検証した。検証３の比較検討例には、図５１の（５１Ｂ）に示した、波長Δｔ−面積ｍ、波長Δｔ−慣性Ｉ、尖度ｋ−波高｜ｈ｜、尖度ｋ−ピーク位置比ｒの４種類の特徴量の組合せを使用した。これらの組合せは、他の組合せと比較して交差検定精度の良い組合せである。個数分析の計算に要する時間には、特徴量作成に要する時間と、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間とが含まれるので、特徴量作成に要する計算時間ＣＴ１、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する計算時間ＣＴ２およびそれらの合計計算時間ＣＴ３（＝ＣＴ１＋ＣＴ２）について比較検討を行った。この場合も、それぞれの所要計算時間は、５０回の交差検定を行って得られた各計算時間の平均値である。

図５２は、４種類の各特徴量組合せに対する合計計算時間ＣＴ３を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ａ）と、各特徴量組合せに対する特徴量作成に要する計算時間ＣＴ１を示すサンプリング周波数（ｋＨｚ）−所要計算時間（秒）のグラフ（５２Ｂ）である。図５３は、各特徴量組合せに対する計算時間ＣＴ２を示すサンプリング周波数−所要計算時間（秒）のグラフである。

（５２Ａ）に示すように、波長Δｔ−面積ｍと波長Δｔ−慣性Ｉの特徴量組合せＧ１は、ほぼ同じ合計計算時間になっており、尖度ｋ−波高｜ｈ｜と尖度ｋ−ピーク位置比ｒの特徴量組合せＧ２は、ほぼ同じ合計計算時間になっている。（５２Ｂ）に示すように、特徴量組合せＧ１のそれぞれの特徴量作成に要する計算時間は同じであり、特徴量組合せＧ２のそれぞれの特徴量作成に要する計算時間は同じになっている。図５３に示すように、Ｈａｓｓｅｌｂｌａｄ法による反復計算に要する時間は、特徴量組合せＧ１、Ｇ２のいずれにおいても、１ＭＨｚ〜１６ｋＨｚでのサンプリング領域において約３，５秒以下の短時間で処理可能になっている。

検証３の特徴量組合せＧ１、Ｇ２の比較結果から明らかに、第１類型と第２類型の同一類型との組合せであっても異なる混合組合せであっても特徴量を使用して所要計算時間の短縮化を図ることができる。したがって、本実施形態に係る個数分析装置によれば。定常
的な個数分析にとどまらず、例えば、緊急性を要する検疫検査や医療現場において、菌類等の粒子有無や個数の判別処理を迅速に行うことができる。

以上の性能評価からわかるように、本実施形態によれば、ナノポアデバイス８により検出した検出信号のデータ群をベースにして、コンピュータ制御プログラム（個数分析プログラム）における、個数導出手段である粒子種分布推定プログラムの実行によって、該検出信号として得られた粒子通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出することができる。したがって、本実施形態に係る個数分析装置を用いることによって、例えば、細菌や微小粒子状物質等の分析物種別に応じた個数ないし個数分布を高精度に分析することができ、個数分析検査における簡易化および低コスト化を実現することができる。ナノポアデバイス８による検出信号を個数分析装置に直接的に取り込んでデータ保存可能にすることにより、検査・分析を統合した粒子種統合分析システムを構築するようにしてもよい。

本実施形態に係る個数分析装置において、特徴量に基づくデータ群から確率密度推定を行い、粒子種別の個数を導出した結果を、出力手段である表示手段７に表示出力あるいはプリンタにプリント出力することができる。したがって、本実施形態によれば、高精度の導出結果（粒子個数、粒子個数分布、推定精度等）を、例えば、ヒストグラムや散布図の出力形態で認知可能に即応的に報知することができるので、例えば、迅速な対応を要する医療現場や検疫場における有用な検査ツールとして本実施形態に係る個数分析装置を使用することができる。

本発明は、個数分析プログラムを搭載した特定のＰＣ等のコンピュータ端末に限らず、該個数分析プログラムの一部ないし全部を記憶した個数分析用記憶媒体に適用することができる。すなわち、所定のコンピュータ端末に該個数分析用記憶媒体に記憶した個数分析プログラムをインストールして所望のコンピュータに個数分析動作させることができるので、簡便かつ安価に個数分析を行うことができる。本発明の適用可能な記憶媒体には、フレキシブルディスク、磁気ディスク、光ディスク、ＣＤ、ＭＯ、ＤＶＤ、ハードディスク、モバイル端末等、コンピュータにより読み取り可能な記憶媒体のいずれかを選択して使用することができる。

尚、本発明は上記実施形態に限定されるものではなく、本発明の技術的思想を逸脱しない範囲における種々変形例、設計変更などをその技術的範囲内に包含するものであることは云うまでもない。

本発明によれば、例えば、細菌や微小粒子状物質等の粒子種別に応じた個数ないし個数分布を高精度に分析することができ、分析物の個数分析検査における簡易化、迅速化および低コスト化に寄与することができる。本発明に係る個数分析技術は、例えば、一刻を争う医療検査に有用であるとともに、感染症状が顕在化する前の微量な細菌・ウイルスの検査にも適応でき、近年注目されている先制医療に対する有力な分析技術にもなり、さらに、迅速な検査結果を必要とする食品の出荷前検査や検疫検査など広範囲に適用することができる。

１ パーソナルコンピュータ
２ ＣＰＵ
３ ＲＯＭ
４ ＲＡＭ
５ データファイル記憶部
６ 入力手段
７ 表示手段
８ マイクロ・ナノポアデバイス
９ チャンバー
１０ 基板
１１ 隔壁
１２ 貫通孔
１３ 電極
１４ 電極
１５ 電源
１６ 増幅器
１７ オペアンプ
１８ 凹部
１９ 帰還抵抗
２０ 電圧計
２１ 検体
２２ 大腸菌
２３ 枯草菌
２４ 電解質溶液

Claims (12)

1. 貫通孔を形成した隔壁と、前記貫通孔を介して前記隔壁の表裏側に配置された電極とを配置し、
   粒子状または分子状の分析物を含む流動性物質を前記隔壁の一面側に供給し、
   前記分析物が前記貫通孔を通過することにより生ずる電極間の通電変化を検出した検出信号のデータに基づいてコンピュータ制御プログラムの実行によって分析物種別の個数を分析する個数分析方法であって、
   前記コンピュータ制御プログラムは、
   前記検出信号として得られた分析物通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、分析物種別の個数を導出する個数導出手段を少なくとも有することを特徴とする個数分析方法。
2. 前記特徴量は、
   前記パルス状信号の波形の局所的特徴を示す第１類型と、
   前記パルス状信号の波形の全体的特徴を示す第２類型のいずれかである請求項１に記載の個数分析方法。
3. 前記第１類型の特徴量は、
   所定の時間幅内における波形の波高値、
   パルス波長ｔ_a、
   パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
   該波形の鋭さを表す尖度、
   パルス開始からパルスピークに至る傾きを表す俯角、
   波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
   パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかである請求項２に記載の個数分析方法。
4. 前記第２類型の特徴量は、
   パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
   前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
   波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
   前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
   波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
   前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかである請求項２に記載の個数分析方法。
5. 前記コンピュータ制御プログラムは、
   前記検出信号のデータないしそれに含まれる揺らぎ成分から、分析物通過のないときの基準ラインを抽出する基準ライン抽出手段と、
   前記基準ラインを基準にして所定範囲を超える信号データを前記パルス状信号のデータ
   として抽出するパルス抽出手段と、
   抽出した前記パルス状信号のデータから前記特徴量を抽出する特徴量抽出手段と、を含む請求項１〜４のいずれかに記載の個数分析方法。
6. 貫通孔を形成した隔壁と、前記貫通孔を介して前記隔壁の表裏側に配置された電極とを配置し、
   粒子状または分子状の分析物を含む流動性物質を前記隔壁の一面側に供給し、
   前記分析物が前記貫通孔を通過することにより生ずる電極間の通電変化を検出した検出信号のデータに基づいてコンピュータ制御プログラムの実行によって分析物種別の個数を分析する個数分析装置であって、
   前記検出信号のデータを記憶する記憶手段と、
   前記検出信号のデータに基づく前記コンピュータ制御プログラムの実行制御を行う制御手段と、を有し、
   前記コンピュータ制御プログラムは、
   前記検出信号として得られた分析物通過に対応するパルス状信号の波形形態の特徴を示す特徴量に基づくデータ群から確率密度推定を行い、分析物種別の個数を導出する個数導出手段を少なくとも有することを特徴とする個数分析装置。
7. 前記特徴量は、
   前記パルス状信号の波形の局所的特徴を示す第１類型と、
   前記パルス状信号の波形の全体的特徴を示す第２類型のいずれかである請求項６に記載の個数分析装置。
8. 前記第１類型の特徴量は、
   所定の時間幅内における波形の波高値、
   パルス波長ｔ_a、
   パルス開始からパルスピークに至るまでの時間ｔ_bとｔ_aとの比ｔ_b／ｔ_aで表わされるピーク位置比、
   該波形の鋭さを表す尖度、
   パルス開始からパルスピークに至る傾きを表す俯角、
   波形を所定の時間毎に区分した時間区分面積の総和を表す面積、および
   パルス開始からパルスピークに至るまでの時間区分面積の和の、全波形面積に対する面積比のいずれかである請求項７に記載の個数分析装置。
9. 前記第２類型の特徴量は、
   パルス開始時点を中心にして前記時間区分面積を質量に、かつ該中心から前記時間区分面積に至る時間を回転半径に擬制したときに定まる時間慣性モーメント、
   前記時間慣性モーメントに対し波高が基準値になるように規格化した場合の規格化された時間慣性モーメント、
   波形を波高方向に等分割し、パルスピーク前後それぞれにおいて各分割単位毎に時刻値の平均値を算出し、同一波高位置の平均値の差をベクトルの成分とする平均値の差ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅平均値慣性モーメント、
   前記波幅平均値慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅平均値慣性モーメント、
   波形を波高方向に等分割し、分割単位毎の時刻値から分散を求め、該分散をベクトルの成分とする分散ベクトルを質量分布と擬制して波形裾野の時間軸を回転中心にしたときに定まる波幅分散慣性モーメント、および
   前記波幅分散慣性モーメントに対し波長が基準値になるように規格化した場合の規格化された波幅分散慣性モーメントのいずれかである請求項７に記載の個数分析装置。
10. 前記コンピュータ制御プログラムは、
    前記検出信号のデータないしそれに含まれる揺らぎ成分から、分析物通過のないときの基準ラインを抽出する基準ライン抽出手段と、
    前記基準ラインを基準にして所定範囲を超える信号データを前記パルス状信号のデータとして抽出するパルス抽出手段と、
    抽出した前記パルス状信号のデータから前記特徴量を抽出する特徴量抽出手段と、を含む請求項６〜９のいずれかに記載の個数分析装置。
11. 前記個数導出手段により導出した個数データを分析物種別に所定の出力態様で出力する出力手段を有する請求項６〜１０のいずれかに記載の個数分析装置。
12. 請求項１に記載の前記コンピュータ制御プログラムを記憶したことを特徴とする個数分析用記憶媒体。

Images