JP6056980B2

JP6056980B2 - ピーク検出方法

Info

Publication number: JP6056980B2
Application number: JP2015535280A
Authority: JP
Inventors: 陽野田; 弘明小澤
Original assignee: Shimadzu Corp
Current assignee: Shimadzu Corp
Priority date: 2013-09-09
Filing date: 2013-09-09
Publication date: 2017-01-11
Anticipated expiration: 2033-09-09
Also published as: JPWO2015033478A1; CN105518455A; WO2015033478A1; EP3029459A1; CN105518455B; US10198630B2; US20160224830A1; EP3029459A4

Description

本発明は、クロマトグラフ分析により得られたクロマトグラムや質量分析により得られたプロファイルスペクトルなどの信号波形に現れるピークを検出するピーク検出方法に関する。

ガスクロマトグラフィや液体クロマトグラフィ等のクロマトグラフ分析により得られるクロマトグラムには、試料に含まれる成分に対応したピークが現れる。通常、そのピークの出現位置（時間）は成分の種類に依存し、ピークの大きさ（つまりその高さや面積）は該ピークに対応する成分の量や濃度に依存する。そのため、クロマトグラムを利用して試料中の成分を同定するにはピークの位置を正確に特定することが重要であり、試料中の成分の量や濃度を把握するにはピークの高さや面積を正確に求めることが重要である。いずれにしても、クロマトグラムに現れるピークの位置やその高さ又は面積値などを求めるには、該クロマトグラムの波形に基づいて有意な成分由来のピークを的確に検出することが必要となる。

クロマトグラムに対する従来一般的なピーク検出の手法には、例えば非特許文献１に記載の方法のように、クロマトグラム波形の傾きを基準としたものが多い。しかしながら、こうしたピーク検出方法では、図５（ａ）に示すようにベースラインの変動があったり、図５（ｂ）に示すように多くのノイズが重畳していたりする場合に、的確にピークを検出することが難しいという問題があった。もちろん、ピーク検出に先立ってベースライン補正やノイズ除去のためのスムージング処理などの前処理を行うことで、ベースライン変動やノイズの影響を或る程度除くことはできるものの、その効果は必ずしも十分であるとはいえない。

また、真のピークを見出すためのクロマトグラム波形の傾きはそのピーク幅に依存するため、的確にピークを検出するには、判定基準である傾きの閾値等のピーク検出パラメータを試料毎に設定する必要がある。しかしながら、代謝解析やバイオマーカ探索といった多量の試料をほぼ連続的に解析しなければならないような場合には、試料毎にピーク検出パラメータを設定することは困難であり、様々な試料に対してピーク検出の精度を安定的に上げることが難しいという問題もあった。

こうした問題を解決し得る、従来とは全く別のピーク検出方法として、ウェーブレット係数空間での隆起線を用いる方法（以下「ウェーブレット隆起線検出方法」という）が提案されている（非特許文献２参照）。このウェーブレット隆起線検出方法を概略的に説明する。

ウェーブレット隆起線検出方法では、質量分析により得られた生のスペクトルデータ（プロファイルスペクトルデータ）に対して直接的に、つまりはベースライン補正やノイズ除去のための前処理を行うことなく、連続ウェーブレット変換（Continuous Wavelet Transform）を適用し、スケールファクタ（Scale factor）を変化させたときのウェーブレット係数を求める。スケールファクタはマザーウェーブレットをウェーブレット関数に変換する際の二つのパラメータのうちの一つであり、マザーウェーブレットを拡大又は縮小するパラメータである。一般に、ウェーブレット係数は、スケールファクタ等のパラメータの下で与えられたウェーブレット関数の成分が元の信号波形にどの程度含まれているかを相対的に示す。このウェーブレット隆起線検出方法では、質量電荷比毎にスケールファクタを変化させたときのウェーブレット係数を計算し、横軸を質量電荷比、縦軸をスケールファクタ、ウェーブレット係数の強度をそれら横軸及び縦軸に共に直交する第３の軸とした３次元の係数空間に描出する。すると、元のプロファイルスペクトル波形上の真のピークに対応する位置に局所的な最大値を示す特徴的な隆起線が観測される。そこで、この隆起線を利用してプロファイルスペクトル波形上のピークを検出する。

プロファイルスペクトル波形のベースライン変動は、狭い時間範囲においては奇関数であるとみなせる。そのため、マザーウェーブレットとして偶関数を用いることによって、ベースライン変動による奇関数成分を相殺することができ、ベースライン補正を予め行うことなく的確にピークを検出することができる。また、スケールファクタを変化させた様々な幅のウェーブレット関数を用いて得られるウェーブレット係数の強度を比較評価することにより、様々なピーク幅を持つピークを的確に検出できるという特徴を有する。

なお、上記のような連続ウェーブレット変換によるウェーブレット係数の算出や３次元表示、さらにはウェーブレット係数の極大値をプロットした隆起線の表示などは、非特許文献３に記載されたような既存のソフトウエアを用いて行うことができる。

「ピーク波形処理を確認しましょう」、株式会社島津製作所、［平成２５年６月７日検索］、インターネット＜URL : http://www.an.shimadzu.co.jp/hplc/support/lib/lctalk/23/23lab.htm＞パン・ドゥ（Pan Du）、ほか２名、「インプルーブド・ピーク・デテクション・イン・マス・スペクトラム・バイ・インコーポレイティング・コンティニュアス・ウェーブレット・トランスフォーム-ベースド・パターン・マッチング（Improved peak detection in mass spectrum by incorporating continuous wavelet transform-based pattern matching）」、バイオインフォメトリクス（Bioinformatics）、2006年、22巻、17号、pp.2059-2065 「Wavelet Toolbox 概説書」、米国マスワークス社（The MathWorks Inc.）、［平成２５年６月１８日検索］、インターネット＜URL : http://www.mathworks.com/tagteam/58032_TT031_Wavelet_Tlbx_Manual.pdf＞

しかしながら、本願発明者の検討によれば、従来提案されているウェーブレット隆起線検出方法では、次のようないくつかの問題がある。
（１）従来のウェーブレット隆起線検出方法は正のピークのみが現れる信号波形を扱うことを前提としているが、クロマトグラフ装置の検出器の種類によっては、クロマトグラムに負のピークが現れることもよくある。こうした負ピークに対応するために、単純に、負の隆起線を扱うという対処法も考えられる。しかしながら、通常の正のピークの裾部付近においても負の隆起線が発生する場合があるため、単純に負の隆起線を扱うようにしただけでは真の負のピークと偽の負のピークとを識別できない。

（２）二つのピークが隣接していて特に一方のピークの裾部と他方のピークの裾部との重なりが大きい場合に、隣接しているピークがウェーブレット係数に相互に影響し、ピーク幅が広く信号強度も大きい側のピークに対応する隆起線に他方のピークに対応する隆起線が引張られてしまう傾向にある。そのため、隆起線に基づいて推定されるピークトップの位置が実際のピークの位置からずれてしまう場合がある。図６にこうした例を示す。この場合、隆起線からはピークトップ位置は時刻ｔ2であると推定されるものの、実際のピークトップ位置は時刻ｔ1であり、分析者がクロマトグラムを見たときに明らかにピークらしくない位置にピークトップが検出されてしまうことになる。

（３）ピークの検出性能を上げるためにはピークとノイズとを分離する必要があり、従来のウェーブレット隆起線検出方法では、ピーク幅にその判別閾値を設定してノイズを除去したり、ＳＮ比を用いてノイズを除去したりすることが提案されている。しかしながら、クロマトグラフ装置においては該装置から出力される検出信号に対し既にローパスフィルタリング処理が施されている場合がよくある。こうした信号波形に対し上述したようなノイズ推定を適用しても、図５（ｂ）に示すようにノイズの中にピークが存在する場合と、同図（ｃ）に示すように真のピークが大量に現れる場合とを識別することは難しい。

本発明はこうした課題を解決するためになされたものであり、その目的とするところは、上述した従来のウェーブレット隆起線検出方法の諸欠点を解消し、与えられた信号波形上のピークを的確に検出することができるピーク検出方法を提供することにある。

上記課題を解決するために成された本発明の第１の態様は、第１のディメンジョンに沿った信号強度の変化を示す信号波形上のピークを検出するピーク検出方法であって、
a)前記信号波形に対し連続ウェーブレット変換を行い、スケールファクタを所定範囲で変化させたときのウェーブレット係数を前記第１のディメンジョン毎に求め、該ウェーブレット係数の強度を第３のディメンジョンとして３次元空間に描画したウェーブレット係数画像上に現れる隆起線に基づいて、前記信号波形に現れる正及び負のピークの候補を抽出するピーク候補抽出ステップと、
b)前記ピーク候補抽出ステップにおいて抽出された負ピーク候補について、該負ピーク候補から推定されるベースラインを基準とした該負ピークの傾き度合いの判定、又は、該負ピークを挟んだ両側の正ピーク候補の有無の判定の少なくともいずれかによって、偽の負ピークを識別して前記負ピーク候補から除外する偽負ピーク除去ステップと、
を備えることを特徴としている。

また上記課題を解決するために成された本発明の第２の態様は、第１のディメンジョンに沿った信号強度の変化を示す信号波形上のピークを検出するピーク検出方法であって、
a)前記信号波形に対し連続ウェーブレット変換を行い、スケールファクタを所定範囲で変化させたときのウェーブレット係数を前記第１のディメンジョン毎に求め、該ウェーブレット係数の強度を第３のディメンジョンとして３次元空間に描画したウェーブレット係数画像上に現れる隆起線に基づいて、前記信号波形に現れる正及び負のピークの候補を抽出するピーク候補抽出ステップと、
b)ピーク候補抽出ステップにおいて抽出されたピーク候補について、前記信号波形を該ピーク候補のピークトップにおいて偶関数化して得られた波形と予め定めたモデルピーク波形との相関係数を計算し、該相関係数に基づいて偽のピークを識別して前記ピーク候補から除外する偽ピーク除去ステップと、

を備えることを特徴としている。

本発明に係るピーク検出方法において、処理対象の信号波形は、典型的には各種クロマトグラフ装置により得られたクロマトグラムである。この場合、第１のディメンジョンとは時間（保持時間）である。また、処理対象の信号波形は、質量分析計で得られるプロファイルスペクトルであってもよい。この場合、第１のディメンジョンとは質量電荷比である。また、処理対象の信号波形は、分光光度計などで得られる吸光スペクトルや反射スペクトルなどであってもよい。この場合、第１のディメンジョンとは波長又は波数である。もちろん、これら例示した以外の信号波形であってもよい。

本発明に係るピーク検出方法において、ピーク候補抽出ステップでは、基本的に非特許文献２に記載されているのと同様の手法により、つまりウェーブレット隆起線検出法により、与えられた信号波形に対するピーク候補を抽出する。ただし、従来のウェーブレット隆起線検出法では正のピークのみを抽出していたのに対し、本発明に係るピーク検出方法では、正のピークと同様の手法により負のピークの候補も抽出する。

なお、非特許文献２に記載のウェーブレット隆起線検出法では、連続ウェーブレット変換に用いるマザーウェーブレットとしてメキシカンハットを用いているが、本発明に係るピーク検出方法では、マザーウェーブレットとして単峰凸関数を正規化したものを用いるのが好ましい。具体的には、ガウス関数を標準偏差σの２倍程度の範囲で打ち切った形状を有する波形など、ピークトップ付近に比べて裾部分を粗くみる波形のマザーウェーブレットを用いるとよい。こうしたマザーウェーブレットでは負値部分が存在しないため、連続ウェーブレット変換に伴う負ピークのアーティファクトの発生を抑えることができるという利点がある。

ピーク候補抽出ステップにより正負両方のピーク候補が抽出されたあと、ピークのアーティファクトによる偽のピークを識別してピーク候補から除外するが、第１の態様と第２の態様とではその手法が相違する。即ち、第１の態様においては、負ピークのみを対象として、偽負ピーク除去ステップにより、負ピーク候補から推定されるベースラインを基準とした該負ピークの傾き度合いの判定、又は、該負ピークを挟んだ両側の正ピーク候補の有無の判定の少なくともいずれかを行うことによって、偽の負ピークを識別する。

例えば、負ピークの傾き度合いが或る閾値を超えている場合には、正のピークの裾部におけるベースラインの傾きが大きいことに起因して生じる偽の負ピークであると判断する。また、上述したように、例えばクロマトグラフ装置では検出器の種類によってはクロマトグラム上に負のピークが現れるものの、正のピークの出現頻度に比べれば負のピークの出現頻度はかなり低い。したがって、例えば負ピークを挟んですぐ両側に正ピーク候補が存在する場合、つまり負ピーク候補が二つの正ピーク候補に挟まれている場合には、二つの正ピーク候補により擬似的に負ピークにみえることで生じる偽の負ピークであると判断する。

一方、第２の態様においては、正若しくは負ピークの一方又は両方を対象として、偽ピーク除去ステップにより、元の信号波形をピーク候補のピークトップにおいて偶関数化して得られた波形と予め定めたモデルピーク波形との相関係数（以下「偶関数化相関係数」という）を計算し、該偶関数化相関係数に基づいて偽のピークを識別する。例えばモデルピーク波形としては適当なガウス関数などを用いればよい。例えば上述したようにベースラインの傾きによって偽のピーク候補が生じている場合、信号波形上でこのピーク候補のピークトップを挟んだ波形形状を偶関数化すると、極端に裾が重くなったり、裾が正負に波打つ形状になる。そのため、上記偶関数化相関係数は小さくなる筈であり、該偶関数化相関係数が所定の閾値以下であれば、これは偽のピークであると判断することができる。一方、傾いたベースライン上に真のピークがある場合は、偶関数化処理後にピーク中の左右対称な成分が抽出される。一般に、抽出後の波形はガウス関数に近い形をしているため、ガウス状関数との相関は高い。本発明では、本偶関数化相関係数は上記課題によって生じるアーティファクトの判定のほか、ピーク形状らしさの特徴量として用いる。

このように本発明に係るピーク検出方法では、ウェーブレット係数の３次元表示に現れる隆起線に基づいて正ピークのみならず負ピークの候補も検出したあとに、偽のピーク候補を的確に除去することができる。それによって、クロマトグラム等の信号波形に出現する正及び負のピークを的確に検出することができる。

また本発明に係るピーク検出方法において、好ましくは、上記信号波形における所定の特徴量に基づいて該信号波形を複数の区間に区分し、該区分毎に求めた信号のＳＮ比を利用して偽ピークを識別して除去するＳＮ比基準偽ピーク除去ステップをさらに有するとよい。また本発明に係るピーク検出方法では、ピーク候補のピーク幅分布を求め、該分布から逸脱するピーク幅を示すピークを偽ピークとして識別するピーク幅基準偽ピーク除去ステップをさらに有するものとするとよい。

これらピーク検出方法によれば、例えば上記第１の態様においても、偽の負ピークだけでなく、偽の正ピークも的確に除去することができる。また、第１、第２態様のいずれにおいても、前述の方法では除去しきれなかった偽ピークを除去できる可能性が高まる。

また本発明に係るピーク検出方法では、上記ピーク候補抽出ステップにより求まるピーク候補の推定位置を、上記信号波形における該推定位置の両側の信号値に基づいて修正するピーク位置修正ステップをさらに有することが好ましい。或いは、上記ピーク候補抽出ステップにより求まるピーク候補の推定位置を、上記信号波形における該推定位置の両側の信号値の二次微分値に基づいて修正するピーク位置修正ステップをさらに有することが好ましい。また、ピークの位置のみならず、ピーク幅についても、上記ピーク位置修正ステップにより修正された位置におけるウェーブレット係数の大きさに基づいて修正するピーク幅修正ステップをさらに有するとよい。

即ち、或るピーク候補について、隆起線から得られるスケールファクタに基づく窓幅（つまりはピーク幅）及びピーク位置（例えば時間）のそれぞれ近傍であるという制約、より具体的には、初期的に推定されるピーク幅及びピーク位置のそれぞれ近傍の所定範囲であるという制約の下で、予め定めたモデルピーク波形との一致具合が最も良くなるようにピーク幅とピーク位置とを調整する。例えば、信号波形上で、最初に推定されたピーク位置よりも右方又は左方において該ピーク位置における強度よりも信号強度が大きい場合に、その方向にピーク位置を所定分だけずらし、その位置での最大のウェーブレット係数を与えるスケールファクタからピーク幅を求める。そして、ウェーブレット係数の強度が或る程度大きな値で収束したときのピーク幅及びピーク位置を、そのピーク候補のピーク幅及びピーク位置であるとすればよい。

これにより、例えば二つのピークが近接していて裾部が重なっているような場合でも、一方のピークにおけるウェーブレット係数に対する他方のピークの波形形状の影響が排除されるので、それぞれのピークの幅及び位置を的確に求めることが可能となる。

また本発明に係るピーク検出方法では、
上記ピーク候補抽出ステップにより求まるピーク候補の中の、任意の一つの着目するピーク候補を中心とする所定範囲に存在する複数のピーク候補の特徴量分布に基づいて、上記着目するピーク候補が、ローパスフィルタ処理されたノイズ波形の中のピークであるか、又は、複数の真のピークからなるピーク群波形の中のピークであるかを判断するピーク種別判定ステップをさらに有し、
該ピーク種別判定ステップによる判定結果に基づいて偽ピークを識別する際の閾値を変更するとよい。

或いは、ピーク種別判定ステップは、一部のピーク候補の特徴量分布ではなく、上記ピーク候補抽出ステップにより求まる全ピーク候補の特徴量分布に基づいて、各ピーク候補が、ローパスフィルタ処理されたノイズ波形の中のピークであるか、又は、複数の真のピークからなるピーク群波形の中のピークであるかを判断するようにしてもよい。

ここでいう「特徴量」とは、例えば、或るピークの裾部に別のピークが重なっている確率（度合い）、ピークの密度、ピークとモデルピーク波形との一致度合いを示す偶関数化相関係数のヒストグラム、などを用いることができる。また、上記ＳＮ比基準偽ピーク除去ステップを有する場合には、このＳＮ比基準偽ピーク除去ステップにおいてＳＮ比に基づいて除去されたピークの度合いを「特徴量」としてもよい。

上述したように、クロマトグラフ装置において検出信号にローパスフィルタ処理を実行した信号波形（クロマトグラム）が出力される構成となっている場合があるが、こうしたローパスフィルタ処理されたノイズはピークの特徴量が所定範囲に収まる傾向にある。そこで、全ピーク候補の特徴量分布、又は所定範囲に存在する複数のピーク候補の特徴量分布に基づいて、或るピーク候補がローパスフィルタ処理されたノイズ波形の中のピークであるか、又は、複数の真のピークからなるピーク群波形の中のピークであるかを、かなり的確に判定することができる。

これにより、例えば白色ノイズにローパスフィルタ処理が施されていてピーク様に見える場合でも、これを真のピークと区別することができる。その結果、ピークの検出精度が向上する。

本発明に係るピーク検出方法によれば、従来のウェーブレット隆起線検出方法では対応不可能である負ピークも的確に検出することができる。また、複数のピークが近接しているような場合や、或いは、信号波形にローパスフィルタ処理が施されていて本来の白色ノイズがピーク様にみえるような場合であっても、ピークを正確に検出することができる。

本発明に係るピーク検出方法を実施するピーク検出装置の一実施例のブロック構成図。本実施例のピーク検出装置で実施されるピーク検出方法の処理の流れを示すフローチャート。図２中のピーク位置及びピーク幅の補正のためのサブルーチンを示すフローチャート。ウェーブレット係数算出に用いるマザーウェーブレットの周波数特性を示す図。様々なピークの状況を示すクロマトグラムの例であり、（ａ）はベースライン変動上にピークが存在する例を示す図、（ｂ）はノイズ中に一つのピークが存在する例を示す図、（ｃ）はピーク群の中にピークが存在する例を示す図。ウェーブレット係数の極値の位置が隣接ピークの影響によってずれる場合の例を示す図。白色ノイズにローパスフィルタ処理を行ったときに得られる信号波形を示す図。クロマトグラムにおいて負ピークのアーティファクトが発生する例を示す図。クロマトグラムにおいて、負ピークが近傍のピークのベースラインに比べて負方向に大きく突出する例を示す図。

以下、本発明に係るピーク検出方法の一例について、添付図面を参照して説明する。
図１は本発明に係るピーク検出方法を実施するピーク検出装置の一実施例のブロック構成図、図２は該装置で実施されるピーク検出方法の処理の流れを示すフローチャートである。

図１に示すように、本実施例のピーク検出装置は、処理対象である信号波形が入力されるウェーブレット変換処理部１、ピーク候補抽出部２、ピーク位置・幅補正部３、偽負ピーク除去部４、ピーク種別判定部５、及び偽ピーク除去部６を備え、信号波形上で検出されたピークの情報が出力される。

この例において、ピーク検出の対象である信号波形は、一般的な液体クロマトグラフ、ガスクロマトグラフ、液体クロマトグラフ質量分析計、液体クロマトグラフ質量分析計などで得られたクロマトグラム波形であるが、本発明に係るピーク検出方法は、クロマトグラム以外の信号波形にも適用可能である。例えば、信号波形は、質量分析計で得られるプロファイルスペクトル、分光光度計などで得られる吸光スペクトルや反射スペクトルなどであってもよい。

本実施例のピーク検出装置は、一般的には、クロマトグラフ装置で収集されたデータをリアルタイムで又はバッチ的に処理するデータ処理装置に包含され、多くの場合、パーソナルコンピュータをハードウエア資源とし、該コンピュータに予めインストールされた専用のデータ処理用ソフトウエアをコンピュータ上で実行することで具現化される。

図２に示すフローチャートに従って、本実施例のピーク検出装置におけるピーク検出動作を説明する。
本実施例のピーク検出装置は、図示しない記憶装置から処理対象であるクロマトグラムデータを読み込み（ステップＳ１）、該データに対して連続ウェーブレット変換を行い、スケールファクタを所定範囲内で且つ所定ステップ幅で以て変化させたときのウェーブレット係数を各時間において算出する（ステップＳ２）。

いま、時間ｔに依存して信号強度が変化するクロマトグラムがｆ(ｔ)であるとし、連続ウェーブレット変換のためのマザーウェーブレットをψ(ｔ)とすると、連続ウェーブレット変換は次の(1)式で定義できる。
Ｃ(ａ，ｂ)＝∫ｆ(ｔ)ψ((ｔ−ｂ)／ａ)ｄｔ …(1)
ここで、Ｃはウェーブレット係数であり、スケールファクタａとシフト係数ｂとの関数として得られる。∫は対象とする時間全体、つまりクロマトグラムの始点から終点までの時間全体の積分である。なお、(1)式は本発明に特有の事項ではなく、ごく一般的な連続ウェーブレット変換の定義である。

上述した非特許文献２に記載の方法では、マザーウェーブレットψ(ｔ)として図４中に示す周波数特性を持つメキシカンハットと呼ばれる関数が用いられている。これに対し、本実施例のピーク検出方法では、図４中に示すガウス（Gauss）関数がマザーウェーブレットψ(ｔ)として用いられる。ただし、ここで用いられるのは完全なガウス関数ではなく、標準偏差σの２倍である２σの範囲で切り出した波形である。なお、ガウス関数でなくても、単峰である凸関数であればよい。

図４から分かるように、メキシカンハットの周波数特性をみると、強度値が両側よりも低い、いわゆる負値部分（図４中のＡ部）が存在する。そのため、マザーウェーブレットψ(ｔ)としてメキシカンハットを用いた場合、このメキシカンハットの負値部分が非常に大きな正ピークに乗じられると、ウェーブレット係数上ではそこに負のピークがあるかのように見えるアーティファクト（負ピークアーティファクト）が出現してしまう。もともと負のピークを考慮しない従来のウェーブレット隆起線検出方法では、負ピークアーティファクトの出現は何ら問題とならないが、本実施例のピーク検出方法では正ピークだけでなく負ピークも検出するため、負ピークアーティファクトによる偽ピークが真の負ピークと識別できなくなると大きな問題である。これに対し、本実施例のピーク検出方法では、負値部分が存在しないガウス関数をマザーウェーブレットψ(ｔ)に用いることで、ウェーブレット変換に起因して負ピークアーティファクトが生じることを回避することができる。

上記のようなウェーブレット変換を行うためには、例えば非特許文献３に記載した「Wavelet toolbox」などを利用することができる。この場合、例えば、マザーウェーブレット及びスケールファクタのレンジを指定し、対象とするクロマトグラムデータに対する連続ウェーブレット変換を実行すると、横軸を時間（シフト係数）とし、縦軸をスケールファクタとした２次元グラフ上に、ウェーブレット係数の強度をカラー表示で表した係数空間における３次元的な画像が得られる。

次に、ピーク候補抽出部２は、ウェーブレット係数の３次元表示において隆起線を検出することで、正ピーク及び負ピークの候補を抽出し、それらピーク候補の位置（時間）とピーク幅の初期値を求める（ステップＳ３、Ｓ４）。ウェーブレット係数の３次元表示における隆起線検出の手法は、非特許文献２に記載されている手法を単純に正・負に拡張したものとする。即ち、大きなスケールファクタ（つまり低周波である）に対応するウェーブレット係数から順に極大値を辿る隆起線により正ピークの候補を求め、一方、同じく極小値を辿る隆起線により負ピークの候補を抽出する。これによって、通常、多数のピーク候補が抽出される。

理想的には、ウェーブレット係数に基づく隆起線はピークの中心（ピークトップ）の位置（この場合には時間）を指す筈である。しかしながら、複数のピークが近接している場合には、その近接するピークによるウェーブレット係数成分が加算された結果、ウェーブレット係数の極値の位置が実際のピーク中心位置からずれてしまい、図６に示したようにピーク位置のずれが発生する。非特許文献２に記載のウェーブレット隆起線検出方法では、隆起線として複数検出されたピーク位置のうち、スケールファクタが最も小さい（つまりはピーク幅が最も小さい）極値がある位置をピーク中心と定めることで、隣接ピークによる極値に対する引き込みの影響を最小限にしている。しかしながら、こうした手法では、高周波ノイズが重畳している場合にピーク位置を正しく推定することができない。

そこで、本実施例のピーク検出方法では、高周波ノイズがあった場合でもより安定的にピーク位置を求めるために、隆起線上で最も大きなウェーブレット係数を与えるスケールファクタを採用する。即ち、理想的にはピーク幅と同じ幅をもつウェーブレット関数となるスケールファクタを採用する。これにより、単一ピークであれば、より周波数の高いウェーブレット関数を用いる場合に比べてピーク位置の精度・安定性は高くなる。しかしながら、前述のように、近接ピークによる隆起線の引き込みの影響があるため、ピーク候補が求まり、ピーク候補毎にピーク位置とピーク幅の初期値が定まったならば、ピーク位置・幅補正部３は、以下の手順で各ピーク候補のピーク位置及びピーク幅を補正する（ステップＳ５）。

図３は図２中のステップＳ５の処理のサブルーチンを示すフローチャートである。
或る一つのピーク候補について、初期位置が与えられると、クロマトグラム波形上でその初期位置の両側の信号強度をみて、正ピークであれば信号強度が高くなる方向（負ピークであれば逆に信号強度が低くなる方向）へと所定距離だけ移動するようにピーク位置を補正する（ステップＳ５１）。例えば図６の例では、ｔ2が初期位置であるとすると、その初期位置から後方側よりも前方側で信号強度は高いから、前方側に移動させるようにピーク位置を補正する。

ただし、ベースラインの傾きが極端に大きい場合には、ピークの端部の信号強度が一番高くなる（つまりはピークトップの信号強度よりも高くなる）ことがあり、その場合には上記のような信号強度に依存する方法を用いることができない。そこで、初期推定位置とピーク幅とに基づき、ピーク位置補正後の位置が不自然に遠いと判断される場合、具体的には、例えば、ピーク位置補正後の位置が初期推定位置とピーク幅とに基づいて求まる所定位置よりも遠くにある場合には、ピーク位置補正を実行しないように制限を掛けるとよい。なお、このとき、好ましくは、信号波形の二次微分を利用し、その変極点を超えるようなピーク位置の補正を行わないように、さらに制限をかけるようにするとよい。

次に、補正後のピーク位置において、最も大きなウェーブレット係数を与えるスケールファクタをピーク幅とするようにピーク幅を補正する（ステップＳ５２）。ただし、実際にはピークの上に液体クロマトグラフにおける送液ポンプの脈動由来のポンプノイズが重畳している場合や、大きなベースライン変動にピークが重畳している場合もあるので、単純に各スケールファクタのウェーブレット係数を比べるのではなく、各ウェーブレット係数にそれぞれ適宜の重みを加えて比較し、最大のスケールファクタを見つけるようにしてもよい。ピーク位置とピーク幅とを補正したならば、例えば信号強度の変化量が或る範囲に収束した、或いは、ピーク幅の変化量が或る閾値以下になった、等を判断することにより、ピーク位置及びピーク幅が或る程度十分に補正されたか否かを判定する（ステップＳ５３）。そして、補正が不足している場合には、ステップＳ５１へと戻り、ステップＳ５１、Ｓ５２の処理を繰り返す。一方、ピーク位置及びピーク幅が或る程度十分に補正されたと判定されると補正を終了する。

通常、ステップＳにおいて抽出されたピーク候補には本来ピークではない偽ピークが多く含まれる。そこでまず、偽負ピーク除去部４は負ピークアーティファクトによる偽負ピークを除去する（ステップＳ６）。
負ピークアーティファクトは通常、図８に示すように、正ピークの裾部分や正ピークと正ピークとの間の谷の部分で発生する。図８左の例は、正ピークがあることによって大きく傾いたベースラインが存在すると誤って認識された結果、正ピークの裾部分がベースラインよりも下方に突出した負ピークであると認識された例である。図８右の例は、ごく近接した信号強度がほぼ同程度の二つの正ピークによって、ピークトップ付近に略水平のベースラインが存在すると誤って認識された結果、正ピークの谷の部分がベースラインよりも下方に突出した負ピークであると認識された例である。

図８左に示した負ピークアーティファクトによる偽ピークは、ベースラインに対しかなり傾いた形状のピークとして検出される。そこで、偽負ピーク除去部４は、ピークが著しく傾いていることを認識することにより、或いは、信号波形上の該ピークをピークトップ付近で偶関数化し、これと予め用意したモデルピークとの相関係数を計算してその係数値を判断することにより、上記のような偽負ピークを検出する。なお、こうした判断では偽の負ピークのみならず偽の正ピークも検出可能であるから、正負にかかわらず偽のピークを除去するようしてもよい。

一方、図８右に示すように、二つの正ピークが近接して存在している場合には、原理的には真の負ピークである可能性もあるものの、クロマトグラムにおいては負ピークが出現する可能性に比べて正ピークが出現する可能性が非常に高いということが分かっているので、一方的に負ピークアーティファクトであると判断すればよい。或いは、偽負ピーク除去部４は、図９に示すように、負ピーク候補が周りのピークのベースライン及びピーク高さに比べて一定強度以上下向きに突出する場合のみ負ピークとして検出し、それに該当しない負ピーク候補は偽であるとみなして除去してもよい。

非特許文献２に記載のウェーブレット隆起線検出方法では、ウェーブレット係数の隆起線を検出した後、高周波数成分の信号エネルギを基準としてＳＮ比を算出し、隆起線で示されるピークが真のピークであるか或いはノイズに起因する偽ピークであるのかを判断している。しかしながら、データを取得するクロマトグラフ装置によっては、ローパスフィルタ処理が施された、つまり高周波成分が低減された信号が出力されることがある。このような場合、実際には周波数の高い白色ノイズが信号波形に重畳されていたとしても、図７に示すように、ピーク波形が多量に集まったようにみえることになる。また、クロマトグラフ装置から出力される信号波形には、ピークとの区別が困難であるポンプノイズや光学的手法を用いた検出器の光源のゆらぎによるノイズ等が重畳していることもあり、こうしたノイズもピーク群にみえる。このような信号に対して従来のウェーブレット隆起線検出方法を適用すると、本来はノイズであるにもかかわらず、全てをピークとして検出してしまうことになる。

これに対し、本実施例のピーク検出装置ではピーク種別判定部５が、入力された信号波形上のピークが、ノイズがローパスフィルタ処理されたものであるか否か、或いは、ポンプノイズ等のピーク群に類似した形状のノイズであるか否か、を判定する。
白色ノイズにローパスフィルタ処理が施されて得られるピーク様の波形や、ポンプノイズ、光源ノイズなどのピーク様の波形（つまりは偽のピーク波形）は、後述するようなピークの特徴量が一定の範囲に収まるという傾向がある。そこで、実際の処理対象として想定される信号波形に基づいてピークの特徴量分布を予め求めておき、それに基づいて設定した適当な閾値で以て処理対象の信号波形から求めた特徴量を判断して、与えられた信号波形中の或る時間範囲がノイズ区間であるか、或いはピーク群区間であるかを弁別している。

ここで利用可能であるピークの特徴量としては、以下のようなものが挙げられる。これら特徴量はいずれも、上述したようなノイズに起因する偽のピーク群である場合には、所定の範囲に収まり、ノイズ由来ではない真のピークであれば、その所定の範囲から逸脱する可能性が高いといえる。
（１）ピークの上に別のピークが重なる確率
本実施例では経験的に、この確率が１５％以下である領域をピーク群領域、５０％以上である領域をノイズ領域と設定している。
（２）ピークの密度（ピーク本数／データ点数、或いはピーク幅積算値／データ点数）
本実施例では経験的に、８０％以上の部分がピークで埋められている場合にその領域がノイズ領域であり得ると設定している。
（３）ピークとモデルピークとの形状の一致度合いを示す偶関数化相関係数のヒストグラム
本実施例では経験的に、偶関数化相関係数が約０．７以上であるヒストグラムの面積が全面積の２０％を下回る場合にその領域がノイズ領域であり得ると設定している。
（４）高周波ノイズ成分で以て正規化したピークの高さのヒストグラム
本実施例では経験的に、ノイズ成分分布から２σ以上離れていて小さいピーク高さと、同様にノイズ成分分布から２σ以上離れていて高いピーク高さとの比が１：３よりも大きい場合に、その領域がノイズ領域であり得ると設定している。
（５）変極点の密度
本実施例では経験的に、変曲点が高周波ノイズの１波長中に１５〜７０％程度あるときにそれがノイズである可能性が高いとみなしている。
もちろん、上述したノイズ判定等のための数値は一例であり、適宜に変更可能である。

ピーク種別判定部５は、具体的には、入力されたクロマトグラム波形を時系列順に適宜の時間範囲に区切り、その時間範囲毎に、そこに含まれるピーク群について上記のような特徴量の少なくとも一つを調べ、その特徴量が予め設定されている閾値範囲内であるか否かを判定することにより、その区間が上記のような特徴的なノイズが存在するノイズ区間であるか、多数のピークが存在するピーク群区間であるか、或いはいずれでもない区間（例えばノイズが小さく明瞭な孤立ピークが存在するような区間）であるのかを識別する。これにより、クロマトグラム波形から抽出された各ピーク候補を、ノイズ中のピーク、ピーク群中のピーク、又はそれ以外の一般のピークのいずれに分類する（ステップＳ７）。

そうして、ピーク候補が分類されたあとに、偽ピーク除去部６は各ピーク候補について偽ピークであるか、真のピークであるのかを判定し、偽ピークを除去する（ステップＳ８）。具体的には、ノイズ区間、ピーク群区間、及び、いずれでもない区間の三つの区間について、その区間に含まれるピーク候補やノイズに基づいて次のような特徴量を求め、その特徴量についての閾値に基づいて偽ピークを識別する。
（Ａ）ピーク高さのヒストグラム
（Ｂ）ピーク幅のヒストグラム
（Ｃ）高周波ノイズレベル
（Ｄ）偶関数化相関係数（ピークとモデルピークとの波形形状の一致度を示す指標）、及び（Ａ）〜（Ｄ）の特徴量を理由に偽ピークであると判断されたピークのヒストグラム

これら特徴量を判断するための閾値は、（Ａ）、（Ｂ）、（Ｄ）などのヒストグラムについては、ピークが分布すると推定される範囲をヒストグラム形状から求め、その範囲を外れるものを偽ピークであるとみなすように、閾値を定めればよい。また、（Ｃ）のノイズレベルについては、一般なＳＮ比に基づく偽ピーク検出のようにノイズレベルの所定数倍のレベルを閾値とすればよい。例えば偽ピークのノイズの高さは理想的にはカイ（χ）二乗分布となるので、高さが低いほうから分布の幅を推定して適宜の閾値を定めればよい。こうして、三つの区間毎に偽ピーク識別の閾値が定められ、その閾値に従って各ピーク候補が偽ピークであるか否かを判断する。

ただし、上記（Ａ）、（Ｂ）のヒストグラムは真のピークのヒストグラムであることが好ましく、偽ピークが除去されて、残ったピーク候補が真である確率が高まるとヒストグラム自体が変化する。また、（Ｄ）のヒストグラムも同様に変化する。即ち、一部ピーク候補が偽ピークであると判断された結果を利用して、さらに偽ピーク、真ピーク毎にピーク幅及びピーク高さのヒストグラム、つまり、より確からしいヒストグラムを作成することができる。
そこで、好ましくは、こうした偽ピーク、真ピーク毎のヒストグラム情報を追加情報として、再び、同様の偽ピーク判定処理を実行することにより、より高精度な偽ピーク判定を行うことができる。即ち、一度の偽ピーク判定処理で以て全ての偽ピークを除去するのではなく、まず粗く偽ピーク判定処理を行って、偽である可能性が特に高いものをまず除去し、その結果を利用して、偽ピーク判定処理自体の精度を高めることで、偽であるか真であるか判断しにくいピーク候補がいずれであるかを判断すればよい。それによって、真のピークを誤って除去してしまうことを防止しつつ、偽ピークを的確に除去することができる。

そうして、偽ピークを除去したならば、残ったピーク候補について、真のピークであるとして、ピーク位置やピーク幅などの情報とともに出力する（ステップＳ９）。
以上のように、本実施例のピーク検出装置では、ウェーブレット隆起線検出を利用して、従来よりも的確にピークを検出することができる。

なお、上記実施例は本発明の一例にすぎず、本発明の趣旨の範囲で適宜、変更、修正、追加を行っても本願特許請求の範囲に包含されることは当然である。

１…ウェーブレット変換処理部
２…ピーク候補抽出部
３…ピーク位置・幅補正部
４…偽負ピーク除去部
５…ピーク種別判定部
６…偽ピーク除去部

Claims

第１のディメンジョンに沿った信号強度の変化を示す信号波形上のピークを検出するピーク検出方法であって、
a)前記信号波形に対し連続ウェーブレット変換を行い、スケールファクタを所定範囲で変化させたときのウェーブレット係数を前記第１のディメンジョン毎に求め、該ウェーブレット係数の強度を第３のディメンジョンとして３次元空間に描画したウェーブレット係数画像上に現れる隆起線に基づいて、前記信号波形に現れる正及び負のピークの候補を抽出するピーク候補抽出ステップと、
b)前記ピーク候補抽出ステップにおいて抽出された負ピーク候補について、該負ピーク候補から推定されるベースラインを基準とした該負ピークの傾き度合いの判定、又は、該負ピークを挟んだ両側の正ピーク候補の有無の判定の少なくともいずれかによって、偽の負ピークを識別して前記負ピーク候補から除外する偽負ピーク除去ステップと、
を備えることを特徴とするピーク検出方法。
第１のディメンジョンに沿った信号強度の変化を示す信号波形上のピークを検出するピーク検出方法であって、
a)前記信号波形に対し連続ウェーブレット変換を行い、スケールファクタを所定範囲で変化させたときのウェーブレット係数を前記第１のディメンジョン毎に求め、該ウェーブレット係数の強度を第３のディメンジョンとして３次元空間に描画したウェーブレット係数画像上に現れる隆起線に基づいて、前記信号波形に現れる正及び負のピークの候補を抽出するピーク候補抽出ステップと、
b)ピーク候補抽出ステップにおいて抽出されたピーク候補について、前記信号波形を該ピーク候補のピークトップにおいて偶関数化して得られた波形と予め定めたモデルピーク波形との相関係数を計算し、該相関係数に基づいて偽のピークを識別して前記ピーク候補から除外する偽ピーク除去ステップと、
を備えることを特徴とするピーク検出方法。
請求項１又は２に記載のピーク検出方法であって、
前記ピーク候補抽出ステップでは、連続ウェーブレット変換に用いるマザーウェーブレットとして単峰凸関数を正規化したものを用いることを特徴とするピーク検出方法。
請求項１〜３のいずれかに記載のピーク検出方法であって、
前記信号波形における所定の特徴量に基づいて該信号波形を複数の区間に区分し、該区分毎に求めた信号のＳＮ比を利用して偽ピークを識別して除去するＳＮ比基準偽ピーク除去ステップをさらに有することを特徴とするピーク検出方法。
請求項１〜３のいずれかに記載のピーク検出方法であって、
ピーク候補のピーク幅分布を求め、該分布から逸脱するピーク幅を示すピークを偽ピークとして識別するピーク幅基準偽ピーク除去ステップをさらに有することを特徴とするピーク検出方法。
請求項１〜５のいずれか記載のピーク検出方法であって、
前記ピーク候補抽出ステップにより求まるピーク候補の推定位置を、前記信号波形における該推定位置の両側の信号値に基づいて修正するピーク位置修正ステップをさらに有することを特徴とするピーク検出方法。
請求項１〜５のいずれか記載のピーク検出方法であって、
前記ピーク候補抽出ステップにより求まるピーク候補の推定位置を、前記信号波形における該推定位置の両側の信号値の二次微分値に基づいて修正するピーク位置修正ステップをさらに有することを特徴とするピーク検出方法。
請求項６又は７に記載のピーク検出方法であって、
前記ピーク候補抽出ステップにより求まるピーク候補の推定ピーク幅を、前記ピーク位置修正ステップにより修正された位置におけるウェーブレット係数の大きさに基づいて修正するピーク幅修正ステップをさらに有することを特徴とするピーク検出方法。
請求項１〜８のいずれかに記載のピーク検出方法であって、
前記ピーク候補抽出ステップにより求まるピーク候補の中の、任意の一つの着目するピーク候補を中心とする所定範囲に存在する複数のピーク候補の特徴量分布に基づいて、前記着目するピーク候補が、ローパスフィルタ処理されたノイズ波形の中のピークであるか、又は、複数の真のピークからなるピーク群波形の中のピークであるのかを判断するピーク種別判定ステップをさらに有し、
該ピーク種別判定ステップによる判定結果に基づいて偽ピークを識別する際の閾値を変更することを特徴とするピーク検出方法。
請求項１〜８のいずれかに記載のピーク検出方法であって、
前記ピーク候補抽出ステップにより求まる全ピーク候補の特徴量分布に基づいて、各ピーク候補が、ローパスフィルタ処理されたノイズ波形の中のピークであるか、又は、複数の真のピークからなるピーク群波形の中のピークであるのかを判断するピーク種別判定ステップをさらに有し、
該ピーク種別判定ステップによる判定結果に基づいて偽ピークを識別する際の閾値を変更することを特徴とするピーク検出方法。
請求項９又は１０に記載のピーク検出方法であって、
前記ピーク種別判定ステップでは、一つのピークの上に他のピークが重畳する度合いを特徴量として複数のピーク候補の又は全ピーク候補の特徴量分布を求めることを特徴とするピーク検出方法。
請求項９又は１０に記載のピーク検出方法であって、
前記信号波形における所定の特徴量に基づいて該信号波形を複数の区間に区分し、該区分毎に求めた信号のＳＮ比を利用して偽ピークを識別して除去するＳＮ比基準偽ピーク除去ステップを有し、
前記ピーク種別判定ステップでは、前記ＳＮ比基準偽ピーク除去ステップにおいてＳＮ比に基づいて除去されたピークの度合いを特徴量として複数のピーク候補の又は全ピーク候補の特徴量分布を求めることを特徴とするピーク検出方法。