基于噪声估计的谱峰识别方法及系统

1.本发明涉及数据处理领域，特别涉及一种基于噪声估计的谱峰识别方法及系统。


背景技术：

2.质谱分析作为高端定量检测分析，在检测的灵敏度、特异性、分析速度、多指标同时检测等方面有非常强的优势。发达国家大力发展医用质谱技术，质谱技术在医用实验室已有30年的应用历史，国外某些大型医用实验室拥有多达几十台质谱仪(如美国arup实验室)，基于质谱的临床检测的项目已达400余项，涉及产前检查、新生儿筛查、滥用药物监测、代谢物检查(氨基酸、脂肪酸)、类固醇激素检测(内分泌)、维生素族检测以及微生物鉴定等领域。国内医用质谱正步入加速发展的快车道，目前已应用于新生儿遗传筛查、维生素d检测、药品检测等检测领域，市场空间预测百亿以上。医用质谱仪中，以液相色谱三重四极杆串联质谱联用仪(lc-ms/ms)应用最为广泛，开展的医学检测项目最多。
3.谱峰检测是lc-ms/ms数据处理最为核心的部分，直接影响仪器的最终分析结果。对谱峰检测方法的优化能以最小的成本带来仪器性能指标的提升。提升谱峰检测方法对于微弱谱峰和重叠峰的检测能力，能够有效提高仪器的灵敏度和分辨率指标。所以，现在有必要提供一种可靠的谱峰检测方法。


技术实现要素：

4.本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种基于噪声估计的谱峰识别方法及系统。本发明的主要原理是色谱峰信号要远远高于噪声信号，通过比较色谱峰信号与噪声的强度值以及信噪比估计值，来确定峰的顶点以及峰的起始点。
5.为实现上述目的，本发明采用的技术方案是：一种基于噪声估计的谱峰识别方法，包括以下步骤：
6.1)输入原始色谱数据r0，r0经预处理后的数据r1再采用s-g平滑算法进行平滑处理得到数据r2；
7.2)针对数据r2，采用寻峰算法识别峰顶点，得到峰顶点集合p；
8.3)计算预处理后的数据r1中的每个点的信噪比估计值；
9.4)寻找峰顶点集合p中每个峰的峰起点：
10.对于峰顶点集合p中的峰顶点pi，在数据r1中寻找与峰顶点pi最近的点，作为数据r1中选择的峰顶点，记为pi
current
；
11.遍历数据r1中点pi
current
左边所有的点，当某一个点pis同时满足以下条件a1、b1、c1时，则将点pis作为峰顶点pi对应的峰起点：
12.a1、点pis的信号强度比峰顶点pi右边1/2个峰宽w范围内的所有点的信号强度均低；
13.b1、点pis到点p
current
的保留时间的差的绝对值小于预先设定的峰宽w；
14.c1、点pis的信噪比估计值大于输入的信噪比阈值t1；
15.遍历峰顶点集合p中所有的峰顶点，按照上述方法寻找得到p中每个峰顶点对应的峰起点；
16.5)寻找每个峰的峰终点：
17.遍历数据r1中点p
current
右边所有的点，当某一个点pid同时满足以下条件a2、b2、c2时，则将点pid作为峰顶点pi对应的峰终点：
18.a2、点pid的信号强度比峰顶点pi左边所有点的信号强度均低；
19.b2、点p
current
到点pid的保留时间的差的绝对值小于预先设定的峰宽w；
20.c2、点pid的信噪比估计值大于输入的信噪比阈值t1；
21.遍历峰顶点集合p中所有的峰顶点，按照上述方法寻找得到p中每个峰顶点对应的峰终点；
22.6)根据每个峰顶点p
current
以及对应的峰起点、峰终点得到数据r1中所有的色谱峰。
23.优选的是，所述步骤1)中，对原始色谱数据r0进行预处理的方法为：
24.对于原始色谱数据r0中的每个数据点，采用每个数据点的前一个和后一个数据点的平均值来代替，处理后的数据即为数据r1，计算公式为：
25.f[r0(x)]＝[r0(x-1)
×
λ1 r0(x)
×
λ2 r0(x 1)
×
λ3]/(λ1 λ2 λ3)；
[0026]
其中，r0(x)表示原始色谱数据r0中的数据点，f[r0(x)]表示预处理后的数据点，λ1、λ2、λ3均为权重系数。
[0027]
优选的是，其中，λ1＝0.5，λ2＝1，λ3＝0.5；
[0028]
优选的是，所述步骤2)中采用的寻峰算法具体包括以下步骤：
[0029]
2-1)输入数据r1，设定信噪比阈值t2；
[0030]
2-2)计算预处理后的数据r1中的每个点的信噪比估计值；
[0031]
2-3)对于数据r1，当其中的一个点同时满足以下条件
①‑⑤
时，则将当前点作为峰顶点：
[0032]
①
当前点的信号强度值同时大于左边第一个点的信号强度值和右边第一个点的信号强度值；
[0033]
②
当前点、当前点左边的第一个点、当前点右边的第一个点的信噪比估计值均大于信噪比阈值t2；
[0034]
③
当前点左边的第二个点的信号强度值小于当前点左边的第一个点的信号强度值；
[0035]
④
当前点右边的第二个点的信号强度值小于当前点右边的第一个点的信号强度值；
[0036]
⑤
当前点左边的第二个点和右边的第二个点的信噪比估计值均大于信噪比阈值t2；
[0037]
2-4)遍历数据r1中的所有点，按照所述步骤2-3)的方法，找出数据r1中的所有峰顶点，构建得到峰顶点集合p。
[0038]
优选的是，所述步骤2)和步骤3)中计算信噪比估计值的方法包括以下步骤：
[0039]
s1、将数据r1分为若干个数据单元，对于每一个数据单元，按照以下方法计算数据单元中的所有数据点的信噪比估计值；
[0040]
s2、计算数据单元中所有数据点的期望e(x)和标准差stdev(x)，其中x表示构成数
据点的向量，记数据选择阈值为ins
max
，ins
max
＝e(x) ηstdev(x)，η为常数；
[0041]
将数据r1的直方图的箱数设置为n
bin
，则直方图被划分为n
bin
段，每一段的范围长度ins
size
为：
[0042][0043]
从而将数据r1中超过ins
max
的数据去除，将数据r1中小于ins
max
的数据点统计到直方图中，直方图的分段区间为：
[0044]
[0,ins
size
)，[ins
size
,2ins
size
)，
……
，[(n
bin-1)ins
size
,n
bin
ins
size
)；
[0045]
将数据r1中小于ins
max
的所有数据点统计到该分段区间中，计算落入到每个分段区间内的数据点的数量；
[0046]
s3、对n
bin
段分段区间按照落入其中的数据点的数量进行排列，筛选出数据点的数量的中位数对应的分段区间[(n
m-1)ins
size
,nmins
size
)，则噪声的估计初始值n0为：
[0047]
n0＝(n
m-0.5)ins
size
；
[0048]
噪声需大于等于1，所以修正之后的噪声估计值n为：
[0049]
n＝max{1,(n
m-0.5)ins
size
}；
[0050]
s4、计算该数据单元中每个点的信噪比估计值ynj：
[0051]
ynj＝yj/n；
[0052]
其中，yj表示该数据单元中的点j的信号值；
[0053]
s4、计算数据r1中的所有数据单元中的数据点的信噪比估计值。
[0054]
优选的是，其中，η＝3。
[0055]
优选的是，其中，n
bin
＝30。
[0056]
本发明还提供一种基于噪声估计的谱峰识别系统，其采用如上所述的方法进行谱峰识别。
[0057]
本发明还提供一种存储介质，其上存储有计算机程序，该程序被执行时用于实现如上所述的方法。
[0058]
本发明还提供一种计算机设备，包括存储器、处理器以及存储在所述存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上所述的方法。
[0059]
本发明的有益效果是：
[0060]
本发明将基于直方图统计的信噪比估计用到峰识别算法中，并改进了直方图统计，采用中值直方图方式，减少了计算复杂度；
[0061]
本发明采用基于噪声估计算法相比较于传统的小波变换算法，具有算法复杂度低、峰识别速度更快且识别精度相当的优点，但更适用于商业应用。
附图说明
[0062]
图1为本发明的实施例中的基于噪声估计的谱峰识别结果。
具体实施方式
[0063]
下面结合实施例对本发明做进一步的详细说明，以令本领域技术人员参照说明书
文字能够据以实施。
[0064]
应当理解，本文所使用的诸如“具有”、“包含”以及“包括”术语并不排除一个或多个其它元件或其组合的存在或添加。
[0065]
实施例1
[0066]
本实施例的一种基于噪声估计的谱峰识别方法，包括以下步骤：
[0067]
1)输入原始色谱数据r0，r0经预处理后的数据r1再采用s-g平滑算法进行平滑处理得到数据r2；
[0068]
其中，对原始色谱数据r0进行预处理的方法为：
[0069]
对于原始色谱数据r0中的每个数据点，采用每个数据点的前一个和后一个数据点的平均值来代替，处理后的数据即为数据r1，计算公式为：
[0070]
f[r0(x)]＝[r0(x-1)
×
λ1 r0(x)
×
λ2 r0(x 1)
×
λ3]/(λ1 λ2 λ3)；
[0071]
其中，r0(x)表示原始色谱数据r0中的数据点，f[r0(x)]表示预处理后的数据点，λ1、λ2、λ3均为权重系数。本实施例中，λ1＝0.5，λ2＝1，λ3＝0.5；
[0072]
savitzky-golay(s-g)平滑算法是由abraham savitzky和marcel j.e.golay在1964年提出的，是一种基于多项式最小二乘拟合的时域滤波方法。其基本思想是设置一个滤波窗口，采用多项式拟合窗内的原始谱数据，使拟合后的数据与原始数据的拟合残差最小。本发明中采用s-g平滑算法进行处理，能获得满意的处理结果。
[0073]
2)针对数据r2，采用寻峰算法识别峰顶点，得到峰顶点集合p；
[0074]
寻峰算法具体包括以下步骤：
[0075]
2-1)输入数据r1，设定信噪比阈值t2；
[0076]
2-2)计算预处理后的数据r1中的每个点的信噪比估计值；
[0077]
2-3)对于数据r1，当其中的一个点同时满足以下条件
①‑⑤
时，则将当前点作为峰顶点：
[0078]
①
当前点的信号强度值同时大于左边第一个点的信号强度值和右边第一个点的信号强度值；
[0079]
②
当前点、当前点左边的第一个点、当前点右边的第一个点的信噪比估计值均大于信噪比阈值t2；
[0080]
③
当前点左边的第二个点的信号强度值小于当前点左边的第一个点的信号强度值；
[0081]
④
当前点右边的第二个点的信号强度值小于当前点右边的第一个点的信号强度值；
[0082]
⑤
当前点左边的第二个点和右边的第二个点的信噪比估计值均大于信噪比阈值t2；
[0083]
2-4)遍历数据r1中的所有点，按照所述步骤2-3)的方法，找出数据r1中的所有峰顶点，构建得到峰顶点集合p。
[0084]
3)计算预处理后的数据r1中的每个点的信噪比估计值；
[0085]
所述步骤2)和步骤3)中计算信噪比估计值的方法包括以下步骤：
[0086]
s1、将数据r1分为若干个数据单元，对于每一个数据单元，按照以下方法计算数据单元中的所有数据点的信噪比估计值；
[0087]
s2、计算数据单元中所有数据点的期望e(x)和标准差stdev(x)，其中x表示构成数据点的向量，记数据选择阈值为ins
max
，ins
max
＝e(x) ηstdev(x)，η为常数；本实施例中，η＝3；
[0088]
将数据r1的直方图的箱数设置为n
bin
，n
bin
＝30；则直方图被划分为n
bin
段，每一段的范围长度ins
size
为：
[0089][0090]
从而将数据r1中超过ins
max
的数据去除，将数据r1中小于ins
max
的数据点统计到直方图中，直方图的分段区间为：
[0091]
[0,ins
size
)，[ins
size
,2ins
size
)，
……
，[(n
bin-1)ins
size
,n
bin
ins
size
)；
[0092]
将数据r1中小于ins
max
的所有数据点统计到该分段区间中，计算落入到每个分段区间内的数据点的数量；
[0093]
s3、对n
bin
段分段区间按照落入其中的数据点的数量进行排列，筛选出数据点的数量的中位数对应的分段区间[(n
m-1)ins
size
,nmins
size
)，则噪声的估计初始值n0为：
[0094]
n0＝(n
m-0.5)ins
size
；
[0095]
噪声需大于等于1，所以修正之后的噪声估计值n为：
[0096]
n＝max{1,(n
m-0.5)ins
size
}；
[0097]
s4、计算该数据单元中每个点的信噪比估计值ynj：
[0098]
ynj＝yj/n；
[0099]
其中，yj表示该数据单元中的点j的信号值；
[0100]
s4、计算数据r1中的所有数据单元中的数据点的信噪比估计值。
[0101]
4)寻找峰顶点集合p中每个峰的峰起点：
[0102]
对于峰顶点集合p中的峰顶点pi，在数据r1中寻找与峰顶点pi最近的点，作为数据r1中选择的峰顶点，记为pi
current
；
[0103]
遍历数据r1中点pi
current
左边所有的点，当某一个点pis同时满足以下条件a1、b1、c1时，则将点pis作为峰顶点pi对应的峰起点：
[0104]
a1、点pis的信号强度比峰顶点pi右边1/2个峰宽w范围内的所有点的信号强度均低；
[0105]
b1、点pis到点p
current
的保留时间的差的绝对值小于预先设定的峰宽w；
[0106]
c1、点pis的信噪比估计值大于输入的信噪比阈值t1；
[0107]
遍历峰顶点集合p中所有的峰顶点，按照上述方法寻找得到p中每个峰顶点对应的峰起点；
[0108]
5)寻找每个峰的峰终点：
[0109]
遍历数据r1中点p
current
右边所有的点，当某一个点pid同时满足以下条件a2、b2、c2时，则将点pid作为峰顶点pi对应的峰终点：
[0110]
a2、点pid的信号强度比峰顶点pi左边所有点的信号强度均低；
[0111]
b2、点p
current
到点pid的保留时间的差的绝对值小于预先设定的峰宽w；
[0112]
c2、点pid的信噪比估计值大于输入的信噪比阈值t1；
[0113]
遍历峰顶点集合p中所有的峰顶点，按照上述方法寻找得到p中每个峰顶点对应的
峰终点；
[0114]
6)根据每个峰顶点p
current
以及对应的峰起点、峰终点得到数据r1中所有的色谱峰。
[0115]
在一种实施例中，t1＝t2。
[0116]
本发明将基于直方图统计的信噪比估计用到峰识别算法中，并改进了直方图统计，采用中值直方图方式，减少了计算复杂度，效果好于最大似然估计算法。
[0117]
最大似然(ml)估计是经典的基于数据辅助信噪比估计算法，其方法可简要概括为根据噪声的概率密度函数求的联合概率密度函数，求解获得有用信号功率和噪声方差估计值，从而获得信噪比估计值。ml估计需要知道信号信息或其对应的估计值，特点是估计准确，但运算量很大。以本发明所涉及仪器信号的噪声为例，仪器噪声概率密度函数属于未知，对未知量的计算则将增加计算量、降低算法效率。本发明使用中值直方图方式，输入相关参数即可对信号做噪声估计，减少了计算复杂度，在lc-ms/ms谱峰识别的应用上效果好于最大似然估计算法。
[0118]
本发明采用基于噪声估计算法相比较于传统的小波变换算法，具有算法时间复杂度低、峰识别速度更快且识别精度相当的优点，但更适用于商业应用。
[0119]
小波变换算法普遍应用于数据的滤波处理，传统的小波变换算法公式为：
[0120][0121]
从变换公式上可发现，小波变换计算的时间复杂度为o(nlog(n))。本发明使用s-g平滑算法用于数据的滤波处理，s-g平滑算法的时间复杂度为o(n)。从算法时间复杂度上看o(n)＜o(nlog(n))，因此本发明算法时间复杂度与小波变换算法相比较低。
[0122]
常见小波变换算法在峰识别应用过程中，应提前对目标峰所在频率做设定才能实现峰识别，获取目标峰的频率范围需对数据做频谱分析。本发明采用信噪比估计算法对峰识别筛选的相关信息做简单规定，峰识别过程所使用参数简单，在识别过程中的速度更快。所以说小波变换算法在峰识别应用时速度与本发明相比较低。
[0123]
本发明基于噪声估计的谱峰识别结果在下图中灰色区域标出，如图1所示：图内识别结果仅包含有效峰信息且未包含无关的噪声信息，说明本发明可对谱峰内无用噪声信息做较好的排除，同时对有效峰信息做保留和识别，本发明谱峰识别精度与小波变换算法相当。
[0124]
实施例2
[0125]
本实施例提供一种基于噪声估计的谱峰识别系统，其采用实施例1的方法进行lc-ms/ms谱峰识别。
[0126]
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节。