CN105718723B - 一种质谱数据处理中谱峰位置检测方法 - Google Patents

Abstract

本发明涉及一种质谱数据处理中谱峰位置检测方法，对原始谱图进行连续小波变换，形成二维小波系数矩阵；在矩阵中设置多个粒子，使粒子按照固定规则移动，逐渐聚集至局部极值点位置；对矩阵位置进行度量，采用噪声阈值和原始谱峰强度修正度量值形成度量值矩阵；在度量值矩阵中搜索脊点，连接脊线，根据脊线确定谱峰位置。该方法综合利用了原始谱峰强度和小波系数矩阵信息进行谱峰位置检测，克服了传统CWT方法进行弱峰和重叠峰检测时误检率上升的问题，提高了算法的灵敏度和对重叠峰的分辨能力。

Description

一种质谱数据处理中谱峰位置检测方法

技术领域

本发明涉及质谱数据预处理以及信息提取方法，特别涉及一种基于连续小波变换的质谱谱峰位置检测方法。

背景技术

质谱法以其高灵敏度、高分辨率等优点成为一种应用广泛的分析技术。质谱计产生原始谱图数据后，需要利用算法进行精确的谱峰位置和强度的检测，以获取样品的成分和结构信息。谱峰位置的检测结果直接影响着质谱计的灵敏度和精度。目前常用的寻峰方法包括直接比较法、导数法、曲线拟合和连续小波变换法(以下简称:CWT)等。直接比较法和导数法对谱图噪声和背景敏感，而曲线拟合法只能提取若干单峰的信息，不具备全谱寻峰能力。

CWT采用特定小波基函数的伸缩和平移来分解谱图,将谱图由时域转换到小波域，形成一个二维的小波系数矩阵，利用小波系数矩阵进行谱峰检测。该方法具有准确率高和多尺度分辨的优点，并且对噪声和背景有较强抑制能力。目前CWT法寻峰通常的做法是：在小波系数矩阵中采用简单的局部比较方式搜索局部极大值确定脊点，将不同尺度上的脊点连接形成脊线，通过脊线确定谱峰位置。该方法单纯根据小波系数确定谱峰位置，在进行微弱谱峰检测时，易受到噪声干扰，存在误检率上升的问题。谱峰重叠时，强峰会对弱峰的小波系数产生严重影响，导致弱峰处的局部极值减弱，甚至成为负值，此时无法实现对弱峰的检测，严重影响了算法对重叠峰的分辨能力。

本发明从文献[1]：R.A.Carmona,W.L.Hwang,and B.Torresani,"Multiridgedetection and time-frequency reconstruction,"Signal Processing,IEEETransactions on,vol.47,pp.480-492,1999.中引入Crazy-climber寻峰方法的思想，根据质谱脊线特征提出了一种新的脊点搜索方式，结合传统CWT方法形成了一种新的谱峰位置检测方法。该方法综合利用了原始谱峰强度和小波系数矩阵信息进行谱峰位置检测，克服了传统CWT方法进行弱峰和重叠峰检测时误检率上升的问题，提高了算法的灵敏度和对重叠峰的分辨能力。

发明内容

本发明的目的是针对上述现有寻峰方法的不足，提出一种质谱数据处理中谱峰位置检测方法。

本发明的目的是通过以下技术方案实现的：

一种质谱数据处理中谱峰位置检测方法，包括以下步骤：

Ⅰ、选择一定尺度范围对原始质谱图进行连续小波变换，形成二维小波系数矩阵C(a,b)，其中a为尺度参数，b为位移参数；

Ⅱ、在小波系数矩阵C(a,b)中均匀设置N个粒子，使粒子按照固定规则移动，逐渐聚集在局部极值点的位置，粒子移动规则如下：

粒子进行t次移动后，粒子位置X(t)＝(k,j)，下一次移动后的位置X(t+1)＝(k’,j’)由以下规则确定：

A、对于小波系数矩阵C(a,b)，以b为水平方向，a为垂直方向。b的范围为1,2,…,B，a的范围为1,2,…,A；

B、粒子首先进行垂直方向移动，如果2≤k≤A-1，那么k’＝k+1和k’＝k-1各取1/2的概率；如果粒子到达下边界，即k＝1，那么粒子向上移动一格，即k’＝2；如果粒子到达上边界，即k＝A，那么粒子向下移动一格，即k’＝A-1；

C、垂直移动完成后进行水平方向移动，仍先判断粒子是否处于左右边界，即j＝1或j＝B，若处于边界则粒子向内移动一格，即j’＝2或j’＝B-1，否则各按1/2的概率准备向左或向右移动一格，是否移动根据小波系数大小进行判断；粒子移动判断标准如下：

①、如果C(k’,j’)>C(k,j)则粒子发生该移动，X(t+1)＝(k’,j’)；

②、如果C(k’,j’)≤C(k,j)则计算差值ΔC＝|C(k’，j’)-C(k，j)|，设置阈值Temp，若ΔC<Temp，则发生该移动，即X(t+1)＝(k’,j’)，否则不动，即X(t+1)＝(k,j)；

③、随移动次数增加，逐渐减小Temp，使Temp逐渐接近零；该方法类似模拟退火算法思想，随移动次数增加，温度降低，粒子趋于稳定；

Ⅲ、按照粒子占据情况对小波系数矩阵中所有位置进行度量，形成度量值矩阵，并采用噪声阈值和对应位置谱峰强度对度量值矩阵进行修正，矩阵位置度量方式如下：

设t次移动后矩阵位置(k,j)的粒子个数为n_(k,j)(t)，粒子移动总次数为T，则矩阵位置初始度量值为：

度量值矩阵修正方法如下：

设阈值i₀为噪声点所产生的度量值，从上述初始度量值中减去阈值，并与原始谱图谱峰强度S(j)作乘积形成最终的度量值矩阵：

I_(k,j)＝S(j)·(i_(k,j)-i₀) (2)

度量值矩阵中所有度量值均为非负实数，上述度量方式中，初始度量值减去噪声阈值i₀后，若度量值小于零则将该值归零，若噪声阈值选择得当，此时度量值矩阵将成为稀疏矩阵，从而有效突出小波系数矩阵中的局部极值点，降低后续脊点搜索难度。

Ⅳ、在修正后的度量值矩阵中进行局部极大值搜索，并设置强度阈值过滤局部极值得到脊点；

Ⅴ、连接脊点形成脊线，设置脊线长度阈值，去除过短脊线；

Ⅵ、根据脊线确定谱峰位置，完成谱峰位置检测，谱峰位置为脊线中度量值最大的脊点位置。

步骤Ⅰ中所述的一定尺度范围a＝1～Na，尺度间隔为指数间隔，最大尺度Na的确定方法包括以下步骤：

a、在所处理的谱图中截取一个独立的谱峰，选择一个大的尺度范围对谱峰进行小波变换，观察谱峰位置处小波系数随尺度参数的变化曲线；

b、随尺度参数增加，小波系数应表现为先增大后减小的趋势，若未出现该趋势，则继续增大尺度范围进行小波变换；

c、出现先增大后减小的趋势后，选择小波系数达到最大后的2-3个尺度点作为最大尺度Na。

有益效果：本发明与一般采用简单的局部比较方式搜索局部极大值确定脊点进行谱峰位置检测的方法不同。提出一种粒子随机移动，根据粒子占据情况度量矩阵位置，采用噪声阈值和谱峰强度修正度量值，在度量值矩阵中进行脊线搜索，确定谱峰位置的方法。综合利用了原始谱峰强度和小波系数矩阵信息进行谱峰位置检测，有效增加微弱谱峰的检测能力。通过粒子占据情况进行矩阵位置度量，采用噪声阈值修正度量值，去除大部分噪声对应的局部极值点，降低误检率。克服了传统CWT方法进行弱峰和重叠峰检测时误检率上升的问题，提高了算法的灵敏度和对重叠峰的分辨能力。

附图说明

图1一种质谱数据处理中谱峰位置检测方法中寻峰方法流程图

图2锆石TOF-SIMS原始谱图

图3一种质谱数据处理中谱峰位置检测方法流程图

图4尺度参数对小波系数的影响图

图5初始度量值图

图6修正度量值图

图7原始谱图及其脊点图

图8寻峰结果图

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细描述：

一种质谱数据处理中谱峰位置检测方法，包括以下步骤：

Ⅰ、选择一定尺度范围对原始质谱图进行连续小波变换，形成二维小波系数矩阵C(a,b)，其中a为尺度参数，b为位移参数；

Ⅱ、在小波系数矩阵C(a,b)中均匀设置N个粒子，使粒子按照固定规则移动，逐渐聚集在局部极值点的位置，粒子移动规则如下：

粒子进行t次移动后，粒子位置X(t)＝(k,j)，下一次移动后的位置X(t+1)＝(k’,j’)由以下规则确定：

A、对于小波系数矩阵C(a,b)，以b为水平方向，a为垂直方向。b的范围为1,2,…,B，a的范围为1,2,…,A；

B、粒子首先进行垂直方向移动，如果2≤k≤A-1，那么k’＝k+1和k’＝k-1各取1/2的概率；如果粒子到达下边界，即k＝1，那么粒子向上移动一格，即k’＝2；如果粒子到达上边界，即k＝A，那么粒子向下移动一格，即k’＝A-1；

C、垂直移动完成后进行水平方向移动，仍先判断粒子是否处于左右边界，即j＝1或j＝B，若处于边界则粒子向内移动一格，即j’＝2或j’＝B-1，否则各按1/2的概率准备向左或向右移动一格，是否移动根据小波系数大小进行判断；粒子移动判断标准如下：

①、如果C(k’,j’)>C(k,j)则粒子发生该移动，X(t+1)＝(k’,j’)；

②、如果C(k’,j’)≤C(k,j)则计算差值ΔC＝|C(k’，j’)-C(k，j)|，设置阈值Temp，若ΔC<Temp，则发生该移动，即X(t+1)＝(k’,j’)，否则不动，即X(t+1)＝(k,j)；

③、随移动次数增加，逐渐减小Temp，使Temp逐渐接近零；该方法类似模拟退火算法思想，随移动次数增加，温度降低，粒子趋于稳定；

Ⅲ、按照粒子占据情况对小波系数矩阵中所有位置进行度量，形成度量值矩阵，并采用噪声阈值和对应位置谱峰强度对度量值矩阵进行修正，矩阵位置度量方式如下：

设t次移动后矩阵位置(k,j)的粒子个数为n_(k,j)(t)，粒子移动总次数为T，则矩阵位置初始度量值为：

度量值矩阵修正方法如下：

设阈值i₀为噪声点所产生的度量值，从上述初始度量值中减去阈值，并与原始谱图谱峰强度S(j)作乘积形成最终的度量值矩阵：

I_(k,j)＝S(j)·(i_(k,j)-i₀) (2)

度量值矩阵中所有度量值均为非负实数，上述度量方式中，初始度量值减去噪声阈值i₀后，若度量值小于零则将该值归零，若噪声阈值选择得当，此时度量值矩阵将成为稀疏矩阵，从而有效突出小波系数矩阵中的局部极值点，降低后续脊点搜索难度。

Ⅳ、在修正后的度量值矩阵中进行局部极大值搜索，并设置强度阈值过滤局部极值得到脊点；

Ⅴ、连接脊点形成脊线，设置脊线长度阈值，去除过短脊线；

Ⅵ、根据脊线确定谱峰位置，完成谱峰位置检测，谱峰位置为脊线中度量值最大的脊点位置。

步骤Ⅰ中所述的一定尺度范围a＝1～Na，尺度间隔为指数间隔，最大尺度Na的确定方法包括以下步骤：

a、在所处理的谱图中截取一个独立的谱峰，选择一个大的尺度范围对谱峰进行小波变换，观察谱峰位置处小波系数随尺度参数的变化曲线；

b、随尺度参数增加，小波系数应表现为先增大后减小的趋势，若未出现该趋势，则继续增大尺度范围进行小波变换；

c、出现先增大后减小的趋势后，选择小波系数达到最大后的2-3个尺度点作为最大尺度Na。

实施例1

采用同位素地质学专用TOF-SIMS锆石谱图中锆同位素谱峰片段(质量数为88-98)对本发明作进一步详细描述，原始谱图如图2所示。

对上述谱图进行谱峰位置检测，检测方法流程图见图3，包括以下步骤：

(1)选择一定尺度范围对原始质谱图进行连续小波变换，形成二维小波系数矩阵C(a,b)，其中a为尺度参数，b为位移参数；位移参数等于原始谱图长度，尺度范围a＝1～Na，优选的尺度间隔为按1.18指数增长，最大尺度Na的确定方法如下：

截取谱图中一个单峰，选择一个大的尺度范围对其进行小波变换，观察谱峰位置处尺度参数对小波系数的影响，如图4所示。从图中可以看出，随着尺度参数的增大，小波系数逐渐增大。当小波基函数与谱峰宽度相同时，达到最大的小波系数值，之后小波系数逐渐减小。为了提高寻峰精度和寻峰算法稳定性，需要选择多个尺度进行小波变换。小的尺度参数与重叠峰和微弱谱峰识别有关，决定了算法的寻峰灵敏度和重叠峰的分辨能力。尺度参数的最大值Na选择小波系数达到最大之后2-3个尺度即可，尺度参数过大只会增大计算量而不会提高精度。本实施例选择的尺度范围为1～19.67共18个值。

(2)在小波系数矩阵C(a,b)中均匀设置N个粒子，使粒子按照固定规则移动，逐渐聚集在局部极值点的位置。粒子移动规则如下：

粒子进行t次移动后，粒子位置X(t)＝(k,j)，下一次移动后的位置X(t+1)＝(k’,j’)由以下规则确定：

A、对于小波系数矩阵C(a,b)，以b为水平方向，a为垂直方向。b的范围为1,2,…,B，a的范围为1,2,…,A；

B、粒子首先进行垂直方向移动，如果2≤k≤A-1，那么k’＝k+1和k’＝k-1各取1/2的概率；如果粒子到达下边界，即k＝1，那么粒子向上移动一格，即k’＝2；如果粒子到达上边界，即k＝A，那么粒子向下移动一格，即k’＝A-1；

C、垂直移动完成后进行水平方向移动，仍先判断粒子是否处于左右边界，即j＝1或j＝B，若处于边界则粒子向内移动一格，即j’＝2或j’＝B-1，否则各按1/2的概率准备向左或向右移动一格，是否移动根据小波系数大小进行判断；粒子移动判断标准如下：

①、如果C(k’,j’)>C(k,j)则粒子发生该移动，X(t+1)＝(k’,j’)；

②、如果C(k’,j’)≤C(k,j)则计算差值ΔC＝|C(k’，j’)-C(k，j)|，设置阈值Temp，若ΔC<Temp，则发生该移动，即X(t+1)＝(k’,j’)，否则不动，即X(t+1)＝(k,j)；

③、随移动次数增加，逐渐减小Temp，使Temp逐渐接近零；该方法类似模拟退火算法思想，随移动次数增加，温度降低，粒子趋于稳定；

(3)根据粒子占据情况对小波系数矩阵中所有位置进行度量，形成度量值矩阵。矩阵位置度量方式如下：

设t次移动后矩阵位置(k,j)的粒子个数为n_(k,j)(t)，粒子移动总次数为T，则矩阵位置初始度量值为：

设置噪声阈值并结合对应位置谱峰强度对初始度量值矩阵进行修正。设阈值i₀为噪声点所产生的度量值，从上述初始度量值中减去阈值，并与原始谱图谱峰强度S(j)作乘积形成最终的度量值矩阵，即：

I_(k,j)＝S(j)·(i_(k,j)-i₀) (4)

采用上述方法对矩阵位置进行度量，尺度参数a＝2.29时的各位移参数对应的初始度量值如图5所示。设置噪声阈值i₀＝4，根据公式(4)对初始度量值进行校正后得到的最终度量值如图6所示。两图比较可以看出，进行度量值修正后，能够将绝大部分度量值归零，使度量值矩阵成为稀疏矩阵。将原始谱峰强度与度量值作乘积，可以有效突出谱峰位置，降低脊点搜索难度。

(4)根据公式(5)在修正后的度量值矩阵中进行局部极大值搜索，并设置强度阈值过滤局部极值得到脊点。原始谱图和得到的脊点如图7所示。

连接脊点形成脊线，设置脊线长度阈值，去除过短脊线。根据脊线确定谱峰位置，完成谱峰位置检测。

优选的脊线长度阈值为5，将脊点数少于5个的脊线剔除。优选的谱峰位置为脊线中度量值最大的脊点位置，度量值越大，说明该位置处小波基函数与谱峰匹配程度最好，因此将其作为谱峰位置。检测到的谱峰位置及过滤后的脊线如图8所示。从图中可以看出，本发明能够有效检测到谱图中各个谱峰位置。包括质量数89处的微弱谱峰和92处的重叠谱峰，都有很好的寻峰效果。

Claims (2)

1.一种质谱数据处理中谱峰位置检测方法，其特征在于，包括以下步骤：

Ⅰ、选择一定尺度范围对原始质谱图进行连续小波变换，形成二维小波系数矩阵C(a,b)，其中a为尺度参数，b为位移参数；

Ⅱ、在小波系数矩阵C(a,b)中均匀设置N个粒子，使粒子按照固定规则移动，逐渐聚集在局部极值点的位置，粒子移动规则如下：

粒子进行t次移动后，粒子位置X(t)＝(k,j)，下一次移动后的位置X(t+1)＝(k’,j’)由以下规则确定：

A、对于小波系数矩阵C(a,b)，以b为水平方向，a为垂直方向，b的范围为1,2,…,B，a的范围为1,2,…,A；

B、粒子首先进行垂直方向移动，如果2≤k≤A-1，那么k’＝k+1和k’＝k-1各取1/2的概率；如果粒子到达下边界，即k＝1，那么粒子向上移动一格，即k’＝2；如果粒子到达上边界，即k＝A，那么粒子向下移动一格，即k’＝A-1；

C、垂直移动完成后进行水平方向移动，仍先判断粒子是否处于左右边界，即j＝1或j＝B，若处于边界则粒子向内移动一格，即j’＝2或j’＝B-1，否则各按1/2的概率准备向左或向右移动一格，是否移动根据小波系数大小进行判断；粒子移动判断标准如下：

①、如果C(k’,j’)>C(k,j)则粒子发生该移动，X(t+1)＝(k’,j’)；

②、如果C(k’,j’)≤C(k,j)则计算差值ΔC＝|C(k’，j’)-C(k，j)|，设置阈值Temp，若ΔC<Temp，则发生该移动，即X(t+1)＝(k’,j’)，否则不动，即X(t+1)＝(k,j)；

③、随移动次数增加，逐渐减小Temp，使Temp逐渐接近零；该方法类似模拟退火算法思想，随移动次数增加，温度降低，粒子趋于稳定；

Ⅲ、按照粒子占据情况对小波系数矩阵中所有位置进行度量，形成度量值矩阵，并采用噪声阈值和对应位置谱峰强度对度量值矩阵进行修正，矩阵位置度量方式如下：

设t次移动后矩阵位置(k,j)的粒子个数为n_(k,j)(t)，粒子移动总次数为T，则矩阵位置 初始度量值为： <mrow> <msub> <mi>i</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>T</mi> </mfrac> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>t</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>T</mi> </msubsup> <msub> <mi>n</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

度量值矩阵修正方法如下：

设阈值i₀为噪声点所产生的度量值，从上述初始度量值中减去阈值，并与原始谱图谱峰强度S(j)作乘积形成最终的度量值矩阵：

I_(k,j)＝S(j)·(i_(k,j)-i₀) (2)

度量值矩阵中所有度量值均为非负实数，上述度量方式中，初始度量值减去噪声阈值i₀后，若度量值小于零则将该值归零，若噪声阈值选择得当，此时度量值矩阵将成为稀疏矩阵，从而有效突出小波系数矩阵中的局部极值点，降低后续脊点搜索难度；

Ⅳ、在修正后的度量值矩阵中进行局部极大值搜索，并设置强度阈值过滤局部极值得到脊点；

Ⅴ、连接脊点形成脊线，设置脊线长度阈值，去除过短脊线；

Ⅵ、根据脊线确定谱峰位置，完成谱峰位置检测，谱峰位置为脊线中度量值最大的脊点位置。

2.按照权利要求1所述的一种质谱数据处理中谱峰位置检测方法，其特征在于，步骤Ⅰ中所述的一定尺度范围a＝1～Na，尺度间隔为指数间隔，最大尺度Na的确定方法包括以下步骤：

a、在所处理的谱图中截取一个独立的谱峰，选择一个大的尺度范围对谱峰进行小波变换，观察谱峰位置处小波系数随尺度参数的变化曲线；

b、随尺度参数增加，小波系数应表现为先增大后减小的趋势，若未出现该趋势，则继续增大尺度范围进行小波变换；

c、出现先增大后减小的趋势后，选择小波系数达到最大后的2-3个尺度点作为最大尺度Na。