CN102495383B - 一种自动相位校正方法 - Google Patents

Abstract

本发明涉及一种自动相位校正方法，首先针对校正好谱图的基线最为平滑的特性，提出了一种新的基于基线连续性的自动校正方法，这种方法可以适应于所有谱图，其适应性非常强但有时不能达到准确的校正效果。为此，在第一步校正后的谱图基础上，利用正峰范围内不存在负数和负峰范围内不存在正数的特性，提出了一种新的基于负值抑制函数的相位校正方法。因此，本发明具有如下优点：可以适应于没有两个独立的单峰信号的谱图；可以适应于数字分辨率和信噪比较差的谱图；可以适用于存在畸变谱峰的压水实验数据和既存在正峰又存在负峰的DEPT-135谱图；相位校正结果精确，自动相位校正的结果和手动相位校正结果只有0.5o以内的偏差。

Description

一种自动相位校正方法

技术领域

本发明涉及一种图像处理方法，尤其是涉及一种自动相位校正方法。适用于适应于没有两个独立的单峰信号、数字分辨率和信噪比较差的谱图的校正。

背景技术

理想情况下，分析仪器如核磁共振波谱仪采集到的时域数据经过Fourier变换得到的会是纯吸收线形谱和纯色散线形谱，但实际上一般得到都是具有相位偏差的谱图数据。造成这种相位偏差的原因有很多，主要有发射机相位和接收机相位不一致引起的零阶相位，发射和接收之间死时间引起的一阶相位。由于纯的吸收线形谱相比于绝对值谱和具有混合相位的谱图有更高的分辨率，所以目前谱图分析都是基于吸收线形谱。而为了得到纯吸收线形谱，必须经过相位校正，所以相位校正是谱图后期处理的一个必须步骤。

相位校正过程是指将Fourier变换后的数据逐点利用如下的公式进行调整，以达到纯吸收线形和纯色散线形：

$R_i=R_i^0*\cos \left({\mathrm{\φ}}_i\right)-I_i^0*\sin \left({\mathrm{\φ}}_i\right)$

$I_i=R_i^0*\sin \left({\mathrm{\φ}}_i\right)+I_i^0*\cos \left({\mathrm{\φ}}_i\right)---\left(1\right)$

其中 ${\mathrm{\φ}}_i=\mathrm{Ph}0+\mathrm{Ph}1*\frac{i}{n}$

上式中R_i和I_i分别表示相位校正后谱图第i个点的实部和虚部，而

和

分别代表校正前谱图第i点的实部和虚部。φ_i代表加在第i个点的相位。n表示谱图数据点数。通常来说，φ为频率的一次函数，包含了与频率无关的零阶相位Ph0和与频率成正比的一阶相位Ph1。

虽然手动相位校正一般都能取得很好的相位校正效果，但其比较耗时，尤其在代谢组学中往往需要批量处理大量的数据，所以自动相位校正在各个谱图数据处理软件中都有实现以方便进行数据处理。早在1969年，Ernst就提出了一种基于色散信号积分为零和吸收信号积分最大的自动相位校正方法。1988年，Marshall提出了一种基于DISPA圆的自动相位校正方法[1]。Heuer提出了利用谱峰对称性进行自动相位校正的方法[2]。Chen提出了基于谱图熵最小和谱峰负值抑制的的一种自动相位校正方法[3]。Brouwer提出了一种结合基线校正，峰面积最小和负值抑制的自动相位校正方法[4]。

虽然上述的自动相位校正都能取得一定的自动相位校正效果，但其都存在一定的不足：(1)Ernst提出的基于色散信号积分为零和吸收信号积分最大的方法对信噪比，数字分辨率以及基线失真较为敏感；(2)DISPA在校正的过程中需要两个独立的单峰信号，而在很多复杂的谱图中这个条件并不能满足；(3)利用谱峰对称性的方法在有些聚合物和复杂的代谢组学的谱图中并不适用，因为其谱峰非常重叠而且并不一定对称。(4)Chen提出的基于谱图熵最小化和Brouwer提出的谱峰面积最小化的方法，由于没有考虑谱图中可能包含负峰或者畸变峰的情况，所以都不能适应DEPT135或者压水峰等谱图。

发明内容

本发明针对现有技术存在的上述技术问题，在于提供了一种自动相位校正方法，可以适应于没有两个独立的单峰信号的谱图；可以适应于数字分辨率和信噪比较差的谱图；可以适用于存在畸变谱峰的压水实验谱图和既存在正峰又存在负峰的DEPT-135谱图；并且其校正结果非常精确，自动相位校正的结果和手动相位校正结果只有0.1％以内的偏差。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种自动相位校正方法，包括以下步骤：

步骤1，输入需要进行相位校正的谱图数据；

步骤2，利用连续小波变换和滑动窗口法识别出步骤1中谱图数据的信号区间，利用结构体Peak存取各个信号区间，结构体Peak中包括信号区间起始点和信号区间终止点两个信息；

步骤3，将步骤1输入的整个谱图数据的谱宽区间减去步骤2中识别出的信号区间，得到各个基线区间，并利用结构体base存取各个基线区间，结构体base中主要包括基线区间起始点和基线区间终止点两个信息；

步骤4，根据步骤3中得到的各个基线区间，利用基线区间的连续性进行第一次自动相位校正，得到所需的校正相位Ph0和Ph1以及经过第一次自动相位校正的谱图；

步骤5，根据步骤4得到的经过第一次自动相位校正的谱图和步骤2识别出的信号区间，将谱图中的谱峰信号区分为正峰信号、负峰信号和畸变峰信号；

步骤6，根据谱图中谱峰中正峰信号中不含有负数点，负峰信号中不含有正数点的特性，对步骤4得到的经过第一次自动相位校正的谱图进行第二次自动相位校正，得到所需改进的相位Ph0_add和Ph1_add；

步骤7，将步骤4得到的相位Ph0及Ph1和步骤6得到的相位Ph0_add及Ph1_add相加，得到最终所需的校正相位Ph0+Ph0_add和Ph1+Ph1_add。

如上所述的步骤2中，获取谱图数据的信号区间的方法为：

步骤2.1，利用连续小波变换对谱图数据进行导数计算，得到谱图数据的数值导数；

步骤2.2，利用滑动窗口法和步骤2.1中计算的谱图的数值导数，识别出谱图数据的信号区间。

如上所述的步骤2.2中，滑动窗口法的步骤如下：

步骤2.2a、根据步骤2.1中得到的谱图数值的数值导数计算出噪声水平σ_noise，采用如下方式计算：首先将数值导数数据分成32份，分别计算每一份的均方差，取其中最小的均方差作为噪声水平σ_noise；

步骤2.2b、设定阈值为n×σ_noise，其中n为参数值，n取值为3～6；

步骤2.3c、比较滑动窗口的高度与步骤2.2b中确定的阈值的大小，滑动窗口的高度大于阈值，则滑动窗口的中心点位于信号区间；滑动窗口的高度小于等于阈值，则滑动窗口的中心点位于基线区间；遇到滑动窗口中心点位于信号区间时，记录这个中心点位置，将其作为该谱峰的起始点，继续滑动窗口直到滑动窗口中心点位于基线区间，并将此时的滑动窗口中心点作为该谱峰的结束点，由此确定一个谱峰区间。

如上所述的步骤4中，进行第一次自动相位校正是基于自定义的衡量基线区间连续性的函数BC，具体方法是遍历每个基线区间，计算当前基线区间的终止点和下一个基线区间起始点的偏差的绝对值，然后将所有的偏差绝对值求和作为基线连续性的量化值，所述函数BC表述如下：

$\mathrm{BC}=\underset{N}{\mathrm{\Σ}}\mathrm{abs}(S(\mathrm{base}(i+1).\mathrm{Start})-S(\mathrm{base}\left(i\right).\mathrm{End}));$

其中：N表示识别出来的基线区间的个数，S代表谱图数据，base(i)和base(i+1)分别表示两个相邻的基线区间，Start为基线区间起始点，End为基线区间终止点，abs为求绝对值。

如上所述的步骤5中将谱图中的谱峰信号区分为正峰信号、负峰信号和畸变峰信号，包含以下步骤：

步骤5.1，遍历谱峰信号区内的谱峰数据点，找出绝对值最大的点作为谱峰峰高；

步骤5.2，谱峰峰高对应点为正数，则先假定该谱峰为准正峰信号，反之为负数，则假定为准负峰信号；

步骤5.3，遍历谱峰数据点，准正峰信号中存在绝对值大于谱峰峰高10％的负数点，该谱峰发生了畸变，将该准正峰信号归入畸变峰信号，同理的准负峰中间存在大于谱峰峰高10％的正数点，将该准负峰归入畸变峰信号。

如上所述的步骤6中，第二次自动相位校正是基于自定义的谱峰特性抑制函数，方法正峰中出现了负数峰点，则将该负数值平方加入抑制函数；负峰中出现正数峰点，就将该值平方后加入抑制函数，该自定义的衡量谱峰特性的函数表述如下：

$P=\underset{N}{\mathrm{\Σ}}P\left(i\right),P\left(i\right)=\left\{\begin{array}{c}\underset{\mathrm{Peak}.\mathrm{Start}}{\overset{\mathrm{Peak}.\mathrm{End}}{\mathrm{\Σ}}}{(S\left(k\right)-\mathrm{abs}\left(S\left(k\right)\right))}^2,\mathrm{if\; peak\; is\; positive}\\ 0,\mathrm{if\; peak\; is\; distorted}\\ \underset{\mathrm{Peak}.\mathrm{Start}}{\overset{\mathrm{Peak}.\mathrm{End}}{\mathrm{\Σ}}}{(-S\left(k\right)-\mathrm{abs}\left(S\left(k\right)\right))}^2,\mathrm{if\; peak\; is\; negative}\end{array}\right.;$

其中N代表谱峰个数，S表示校正好谱图的实部数据，abs表示求绝对值，Peak.Start为一个信号区间起始点，Peak.End为一个信号区间终止点，abs为求绝对值。

本发明与现有技术相比，具有以下优点：

1.由于本方法对谱图中的谱峰信号没有要求，其适应性更广，可以适应于没有两个独立的单峰信号的谱图；

2.由于本方法基于基线的连续性和谱峰的整体特性，所以其相对于基于谱峰对称性等自动相位校正方法来说可以更适应于数字分辨率和信噪比较差的谱图；

3.由于本方法中的第一次相位校正是基于基线的连续性，所以其可以适用于存在畸变谱峰的压水实验数据和既存在正峰又存在负峰的DEPT-135谱图；

4.校正结果精确，自动相位校正的结果和手动相位校正结果只有0.5°以内的偏差，下表表示各种不同的谱图数据手动相位校正和自动相位校正偏差，表1中给出了手动相位校正和自动相位校正偏差对比。

表1手动相位校正和自动相位校正偏差

附图说明

图1是一种自动相位校正方法的流程图。

图2是利用不同的自动相位校正方法处理压水峰核磁共振谱图的结果示意图；(a)未经过相位校正的压水峰NMR谱图数据；(b)利用商用软件Topspin2.1相位校正结果；(c)利用本发明中提出的自动相位校正方法处理后校正结果。

图3是将这种自动相位校正方法应用于DEPT-135核磁共振谱图的结果示意图；其中(a)未经过相位校正的DEPT-135核磁共振谱图；(b)利用本发明中提出的自动相位校正方法处理后校正结果。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

如图1所示，一种自动相位校正方法，包括以下步骤：

步骤1，输入需要进行相位校正的谱图数据；

步骤2，利用连续小波变换和滑动窗口法识别出步骤1中谱图数据的信号区间，利用结构体Peak存取各个信号区间，结构体Peak中包括信号区间起始点和信号区间终止点两个信息；

步骤3，将步骤1输入的整个谱图数据的谱宽区间减去步骤2中识别出的信号区间，得到各个基线区间，并利用结构体base存取各个基线区间，结构体base中主要包括基线区间起始点和基线区间终止点两个信息；

步骤4，根据步骤3中得到的各个基线区间，利用基线区间的连续性进行第一次自动相位校正，得到所需的校正相位Ph0和Ph1以及经过第一次自动相位校正的谱图；

步骤5，根据步骤4得到的经过第一次自动相位校正的谱图和步骤2识别出的信号区间，将谱图中的谱峰信号区分为正峰信号、负峰信号和畸变峰信号；

步骤6，根据谱图中谱峰中正峰信号中不含有负数点，负峰信号中不含有正数点的特性，对步骤4得到的经过第一次自动相位校正的谱图进行第二次自动相位校正，得到所需改进的相位Ph0_add和Ph1_add；

步骤7，将步骤4得到的相位Ph0及Ph1和步骤6得到的相位Ph0_add及Ph1_add相加，得到最终所需的校正相位Ph0+Ph0_add和Ph1+Ph1_add。

步骤2中，获取谱图数据的信号区间的方法为：

步骤2.1，利用连续小波变换对谱图数据进行导数计算，得到谱图数据的数值导数；

步骤2.2，利用滑动窗口法和步骤2.1中计算的谱图的数值导数，识别出谱图数据的信号区间。

步骤2.2中，滑动窗口法的步骤如下：

步骤2.2a、根据步骤2.1中得到的谱图数值的数值导数计算出噪声水平σ_noise，采用如下方式计算：首先将数值导数数据分成32份，分别计算每一份的均方差，取其中最小的均方差作为噪声水平σ_noise；

步骤2.2b、设定阈值为n×σ_noise，其中n为参数值，n取值为3～6；

步骤2.3c、比较滑动窗口的高度与步骤2.2b中确定的阈值的大小，滑动窗口的高度大于阈值，则滑动窗口的中心点位于信号区间；滑动窗口的高度小于等于阈值，则滑动窗口的中心点位于基线区间；遇到滑动窗口中心点位于信号区间时，记录这个中心点位置，将其作为该谱峰的起始点，继续滑动窗口直到滑动窗口中心点位于基线区间，并将此时的滑动窗口中心点作为该谱峰的结束点，由此确定一个谱峰区间。

步骤4中，进行第一次自动相位校正是基于自定义的衡量基线区间连续性的函数BC，具体方法是遍历每个基线区间，计算当前基线区间的终止点和下一个基线区间起始点的偏差的绝对值，然后将所有的偏差绝对值求和作为基线连续性的量化值，所述函数BC表述如下：

$\mathrm{BC}=\underset{N}{\mathrm{\Σ}}\mathrm{abs}(S(\mathrm{base}(i+1).\mathrm{Start})-S(\mathrm{base}\left(i\right).\mathrm{End}));$

其中：N表示识别出来的基线区间的个数，S代表谱图数据，base(i)和base(i+1)分别表示两个相邻的基线区间，Start为基线区间起始点，End为基线区间终止点，abs为求绝对值。

步骤5中将谱图中的谱峰信号区分为正峰信号、负峰信号和畸变峰信号，包含以下步骤：

步骤5.1，遍历谱峰信号区内的谱峰数据点，找出绝对值最大的点作为谱峰峰高；

步骤5.2，谱峰峰高对应点为正数，则先假定该谱峰为准正峰信号，反之为负数，则假定为准负峰信号；

步骤5.3，遍历谱峰数据点，准正峰信号中存在绝对值大于谱峰峰高10％的负数点，该谱峰发生了畸变，将该准正峰信号归入畸变峰信号，同理的准负峰中间存在大于谱峰峰高10％的正数点，将该准负峰归入畸变峰信号。

步骤6中，第二次自动相位校正是基于自定义的谱峰特性抑制函数，方法正峰中出现了负数峰点，则将该负数值平方加入抑制函数；负峰中出现正数峰点，就将该值平方后加入抑制函数，该自定义的衡量谱峰特性的函数表述如下：

$P=\underset{N}{\mathrm{\Σ}}P\left(i\right),P\left(i\right)=\left\{\begin{array}{c}\underset{\mathrm{Peak}.\mathrm{Start}}{\overset{\mathrm{Peak}.\mathrm{End}}{\mathrm{\Σ}}}{(S\left(k\right)-\mathrm{abs}\left(S\left(k\right)\right))}^2,\mathrm{if\; peak\; is\; positive}\\ 0,\mathrm{if\; peak\; is\; distorted}\\ \underset{\mathrm{Peak}.\mathrm{Start}}{\overset{\mathrm{Peak}.\mathrm{End}}{\mathrm{\Σ}}}{(-S\left(k\right)-\mathrm{abs}\left(S\left(k\right)\right))}^2,\mathrm{if\; peak\; is\; negative}\end{array}\right.;$

其中N代表谱峰个数，S表示校正好谱图的实部数据，abs表示求绝对值，Peak.Start为一个信号区间起始点，Peak.End为一个信号区间终止点，abs为求绝对值。

如图2所示，其中(a)未经过相位校正的代谢组学NMR谱图数据；(b)利用商用软件Topspin2.1相位校正结果；(c)利用本发明中提出的相位校正方法的校正结果。图中X轴为化学定标所用的ppm，Y轴为谱图数值，可以看出采用商用软件自动相位后的谱图在0ppm附近的谱峰仍存在相位失真，而采用我们自动相位校正结果非常接近手动相位校正结果。

如图3所示，其中(a)未经过相位校正的DEPT-135核磁共振谱图数据；(b)利用本发明中提出的自动相位校正方法处理后校正结果。图中X轴为化学定标所用的ppm，Y轴为谱图数值，从图中可以看出我们提出的自动相位校正方法也可适用于既存在正峰又存在负峰的DEPT-135核磁共振谱图数据。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种自动相位校正方法，其特征在于，包括以下步骤：

步骤1，输入需要进行相位校正的谱图数据；

步骤2，利用连续小波变换和滑动窗口法识别出步骤1中谱图数据的信号区间，利用结构体Peak存取各个信号区间，结构体Peak中包括信号区间起始点和信号区间终止点两个信息；

步骤3，将步骤1输入的整个谱图数据的谱宽区间减去步骤2中识别出的信号区间，得到各个基线区间，并利用结构体base存取各个基线区间，结构体base中主要包括基线区间起始点和基线区间终止点两个信息；

步骤4，根据步骤3中得到的各个基线区间，利用基线区间的连续性进行第一次自动相位校正，得到所需的校正相位Ph0和Ph1以及经过第一次自动相位校正的谱图；

步骤5，根据步骤4得到的经过第一次自动相位校正的谱图和步骤2识别出的信号区间，将谱图中的谱峰信号区分为正峰信号、负峰信号和畸变峰信号；

步骤6，根据谱图中谱峰中正峰信号中不含有负数点，负峰信号中不含有正数点的特性，对步骤4得到的经过第一次自动相位校正的谱图进行第二次自动相位校正，得到所需改进的相位Ph0_add和Ph1_add；

步骤7，将步骤4得到的相位Ph0及Ph1和步骤6得到的相位Ph0_add及Ph1_add相加，得到最终所需的校正相位Ph0+Ph0_add和Ph1+Ph1_add；

其中，步骤2中获取谱图数据的信号区间的方法为：

步骤2.1，利用连续小波变换对谱图数据进行导数计算，得到谱图数据的数值导数；

步骤2.2，利用滑动窗口法和步骤2.1中计算的谱图的数值导数，识别出谱图数据的信号区间；

其中，步骤2.2中滑动窗口法的步骤如下：

步骤2.2a、根据步骤2.1中得到的谱图数值的数值导数计算出噪声水平σ_noise，采用如下方式计算：首先将数值导数数据分成32份，分别计算每一份的均方差，取其中最小的均方差作为噪声水平σ_noise；

步骤2.2b、设定阈值为n×σ_noise，其中n为参数值，n取值为3～6；

步骤2.3c、比较滑动窗口的高度与步骤2.2b中确定的阈值的大小，滑动窗口的高度大于阈值，则滑动窗口的中心点位于信号区间；滑动窗口的高度小于等于阈值，则滑动窗口的中心点位于基线区间；遇到滑动窗口中心点位于信号区间时，记录这个中心点位置，将其作为该谱峰的起始点，继续滑动窗口直到滑动窗口中心点位于基线区间，并将此时的滑动窗口中心点作为该谱峰的结束点，由此确定一个谱峰区间。

2.根据权利要求1所述的一种自动相位校正方法，其特征在于，所述的步骤4中，进行第一次自动相位校正是基于自定义的衡量基线区间连续性的函数BC，具体方法是遍历每个基线区间，计算当前基线区间的终止点和下一个基线区间起始点的偏差的绝对值，然后将所有的偏差绝对值求和作为基线连续性的量化值，所述函数BC表述如下：

$\mathrm{BC}=\underset{N}{\mathrm{\Σ}}\mathrm{abs}(S(\mathrm{base}(i+1).\mathrm{Start})-S(\mathrm{base}\left(i\right).\mathrm{End}));$

其中：N表示识别出来的基线区间的个数，S代表谱图数据，base(i)和base(i+1)分别表示两个相邻的基线区间，Start为基线区间起始点，End为基线区间终止点，abs为求绝对值。

3.根据权利要求1所述的一种自动相位校正方法，其特征在于，所述的步骤5中将谱图中的谱峰信号区分为正峰信号、负峰信号和畸变峰信号，包含以下步骤：

步骤5.1，遍历谱峰信号区内的谱峰数据点，找出绝对值最大的点作为谱峰峰高；

步骤5.2，谱峰峰高对应点为正数，则先假定该谱峰为准正峰信号，反之为负数，则假定为准负峰信号；

步骤5.3，遍历谱峰数据点，准正峰信号中存在绝对值大于谱峰峰高10%的负数点，该谱峰发生了畸变，将该准正峰信号归入畸变峰信号，同理的准负峰中间存在大于谱峰峰高10%的正数点，将该准负峰归入畸变峰信号。

4.根据权利要求1所述的一种自动相位校正方法，其特征在于，所述的步骤6中，第二次自动相位校正是基于自定义的谱峰特性抑制函数，方法正峰中出现了负数峰点，则将该负数值平方加入抑制函数；负峰中出现正数峰点，就将该值平方后加入抑制函数，该自定义的衡量谱峰特性的函数表述如下：

$P=\underset{N}{\mathrm{\Σ}}P\left(i\right)P\left(i\right)=\left\{\begin{array}{c}\underset{\mathrm{Peak}.\mathrm{Start}}{\overset{\mathrm{Peak}.\mathrm{End}}{\mathrm{\Σ}}}{(S\left(k\right)-\mathrm{abs}\left(S\left(k\right)\right))}^2,\mathrm{ifpeakispositive}\\ 0,\mathrm{ifpeakisdistorted}\\ \underset{\mathrm{Peak}.\mathrm{Start}}{\overset{\mathrm{Peak}.\mathrm{End}}{\mathrm{\Σ}}}{(-S\left(k\right)-\mathrm{dbs}\left(S\left(k\right)\right))}^2,\mathrm{ifpeakisnegative}\end{array}\right.;$

其中N代表谱峰个数，S表示校正好谱图的实部数据，abs表示求绝对值,Peak.Start为一个信号区间起始点，Peak.End为一个信号区间终止点，abs为求绝对值。

Images