CN103018194A - 基于背景估计的非对称最小二乘基线校正方法 - Google Patents

Abstract

本发明公开了一种基于背景估计的非对称最小二乘基线校正方法，包括:利用直方图背景估计方法计算出的背景值做为基线的初始值；建立改进的非对称最小二乘基线校正模型；迭代地求解该模型；判断前后两次迭代光谱的背景值，如果几乎不改变，则整个算法终止，输出最终的基线。本发明所设计的基于背景估计的非对称最小二乘基线校正方法操作简单，方便适用，可靠性高，针对性强，在保持红外谱图峰位和峰形不变的同时，有效地消除了谱图的基线，取得了较好的效果。

Description

基于背景估计的非对称最小二乘基线校正方法

技术领域

本发明涉及红外光谱化学计量学定性/定量分析技术领域，特别是涉及一种基于背景估计的非对称最小二乘基线校正方法。

背景技术

谱图采集过程中会受到背景温度、湿度等外界环境以及仪器自身工作状况的影响，导致最终获取的光谱图呈现不同程度的基线漂移。基线漂移一方面会影响定性分析的精度，另一方面会导致光谱峰高和峰面积测量失真，从而影响定量分析的精度。另外，即使是相似的光谱样本，其基线也可能有很大的不同，从而导致基于这些光谱样本建立的数学模型缺乏鲁棒性。为了建立一个稳定、可靠的定性分析或定量分析模型，必须进行基线校正。

基线校正有人工校正和自动校正两种。人工校正一般是通过人机两互的方式选取谱图上的特征点，然后将这些特征点拟合成一条曲线。这种方法费时、费力，且由于人工选择的特征点随意性较大，其重现性较差。Schulze等人对基线自动校正算法进行了全面的综述[1]，目前主要有微分和滤波方法[2]、形态学方法[3]、插值拟合方法[4]、背景估计方法[5~8]，等等。但是这些方法都有不足之处。微分方法可以消除常数和线性漂移，但同时也放大了谱图中的高频噪声；形态学方法可以有效消除基线基线漂移，但是由于红外谱图中各吸收峰的峰宽不一致，其结构元素的选择十分棘手；插值拟合方法中，插值节点的自动选取需要依靠人工经验，并且不同的插值函数拟合出的基线也不尽相同；背景估计方法如SNIP[8]，其方法快速简便，但迭代次数的选择是个难题，过大过小都不合适。

非对称最小二乘基线校正算法是Eilers等人2003年提出来的[9]，它具有严格的数学理论基础，在处理谱图时有较好的效果，因而受到了广泛的关注。该方法基于Whittaker smoother估计基线[10~11]，要求待求基线非对称地拟合给定光谱，同时对待求基线施加某种光滑性约束。但是该算法只考虑了二阶导数的平滑性约束，实际上在基线尽可能平滑的前提下，我们不仅要求拟合出的数值与原始数据之间的误差很小，而且还要求它们的一阶导数很接近。因此，基线校正直接应用非对称最小二乘算法是有局限性的，必须充分考虑红外光谱自身的物理特性。

[1]Schulze G,Jirasek A,Yu M M L,et al.Investigation of SelectedBaseline Removal Techniques as Candidates for Automated Implementation.Applied spectroscopy,2005,59(5):545-574

[2]Leger M N and Ryder a G.Comparison of Derivative Preprocessingand Automated Polynomial Baseline Correction Method for Classificationand Quantification of Narcotics in Solid Mixtures.Applied spectroscopy,2006,60(2):182-193

[3]P.Maragos,R.W.Schafer.Morphological Filters.Part 1.TheirSet-Theoretic Analysis and Relations to Linear Shift-Invariant Filters.IEEE Trans.on Acoustics,Speech,and Signal Processing,VOL.35,NO.8,1987

[4]Brown D E.Fully Automated Baseline Correction of 1d and 2d NmrSpectra Using Bernstein Polynomials.Journal of Magnetic Resonance,Series A,1995,114(2):268-270

[5]Andrade L and Manolakos E S.Signal Background Estimation andBaseline Correction Algorithms for Accurate DNA Sequencing.The Journalof VLSI Signal Processing,2003,35(3):229-243

[6]Andrade L and Manolakos E S.Accurate Estimation of the SignalBaseline in DNA Chromatograms.2002.35-44

[7]Marion D and Bax A.Baseline Correction of 2d Ft Nmr Spectra Usinga Simple Linear Prediction Extrapolation of the Time-Domain Data.J.Magn.Reson,1989,83(1):205-211

[8]Ryan C G,Clayton E J,Griffin W L,et al.Snip,aStatistics-Sensitive Background Treatment for the Quantitative Analysisof Pixe Spectra in Geoscience Applications.Nuclear Instruments andMethods in Physics Research Section B,1988,34(3)

[9]P.H.C.Eilers,H.F.M.Boelens.Baseline Correction withAsymmetric Least Squares Smoothing,2005,10

[10]Newey W K and Powell J L.Asymmetric Least Squares Estimation andTesting.Econometrica:Journal of the Econometric Society,1987:819-847

[11]Eilers P H C.A Perfect Smoother.Anal.Chem,2003,75(14):3631-3636

发明内容

本发明所要解决的技术问题是提供一种操作简单，方便适用，可靠性高，针对性强，在保持红外谱图峰位和峰形不变的同时，有效地消除了谱图的基线的基于背景估计的非对称最小二乘基线校正方法。

为了解决上述的技术问题：本发明设计了一种基于背景估计的非对称最小二乘基线校正方法，包括如下具体步骤：

步骤（1）：利用直方图背景估计方法计算出的背景值做为基线的初始值c_i；

步骤（2）：建立改进的非对称最小二乘基线校正模型；

步骤（3）：迭代地求解步骤（2）中的改进的非对称最小二乘基线校正模型，得到迭代后的光谱的背景值c⁽ⁱ⁺¹⁾；

步骤（4）：比较迭代前后光谱的背景值，如果几乎不改变，则整个算法终止，输出最终的基线。

作为本发明的一种优化方法：所述步骤（1）包括如下具体处理：

步骤（11）：初始化i←0,X₀←{窗口中所有样本}；

步骤（12）：计算X_i的标准差σ_i；

步骤（13）：计算柱子数目n←round((max(X_i)-min(X_i))/σ_i)；

步骤（14）：计算柱子大小s←(max(X_i)-min(X_i))/n；

步骤（15）：计算直方图H←hist(X_i,n,s)；

步骤（16）：找出直方图中最大柱子的中心c_i；

步骤（17）：更新窗口，X_i+1←{x∈X_i|-2σ_i≤x-c_i≤2σ_i}；

步骤（18）：判断停止条件，如果|X_i+1|>0.95|X_i|，转步骤（110）；

步骤（19）：i＝i+1，返回步（2）；

步骤（110）：输出c_i。

作为本发明的一种优化方法：所述步骤（2）包括如下具体处理：

步骤（21）：结合红外光谱自身的物理特性，找到谱图的变化特征，建立优化函数：

$\arg \underset{z}{\min }[\underset{i}{\mathrm{\Σ}}{w}_i{(y_i-z_i)}^1+\mathrm{\λ}\underset{i}{\mathrm{\Σ}}{\left({\mathrm{\Δ}}^2{z}_i\right)}^2]$

其中，w_i定义为权重因子，Δ定义为二阶差分算子，λ定义为第一规整化参数；

步骤（22）：在步骤（21）的函数中加入一阶导数约束项，最小化目标函数如下：

L(z)=||Q(y-z)||²+λ₁||D₁(y-z)||²+λ||Dz||²

其中D₁定义为一阶微分矩阵，D定义为二阶微分矩阵，Q定义为一个对角权重矩阵，λ₁定义为第二规整化参数。

作为本发明的一种优化方法：所述步骤（3）包括如下具体处理：

步骤（31）：将步骤（22）中的最小化目标函数L(z)最终转化为迭代求解线性方程***：

$(Q^{T}Q+{\mathrm{\λ}}_1{D}_1^T{D}_1+{\mathrm{\λD}}^{T}D)z=(Q^{T}Q+{\mathrm{\λ}}_1{D}_1^T{D}_1)y;$

步骤（32）：运用直方图背景估计方法估计出原始光谱y的背景值c⁽⁰⁾，基于背景值c⁽⁰⁾得到基线的初始估计值z⁽⁰⁾，利用基线的初始估计值z⁽⁰⁾来确定校正光谱y-z⁽⁰⁾为负值的位置，第一次迭代时，构造如下矩阵：

将上述矩阵代入步骤（31）中的迭代求解线性方程***，求解出新的基线z⁽¹⁾；

步骤（33）：校正光谱的负值部分不断改变，这种迭代的基线求解算法可以表示为：

$\left\{\begin{array}{c}{Q}^{\left(k\right)}=\mathrm{diag}(y-z^{(k-1)}<0)\\ z^{\left(k\right)}={(Q^{\left(k\right)T}{Q}^{\left(k\right)}+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1+{\mathrm{\&lambda;D}}^{T}D)}^{-1}(Q^{\left(k\right)T}{Q}^{\left(k\right)}+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1)y\end{array}\right.;$

步骤（34）：重复步骤（31）至步骤（33）N次达到收敛，其中，N定义为自然数。

作为本发明的一种优化方法：所述步骤（4）包括如下具体处理：

步骤（41）：利用直方图背景估计法，计算校正光谱y-z_(i)的背景值c⁽ⁱ⁺¹⁾，

如果|c⁽ⁱ⁺¹⁾-c⁽ⁱ⁾|<ε，则算法终止，得到基线，以及基线校正后的谱图y。

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

1.本发明所设计的基于背景估计的非对称最小二乘基线校正方法所得基线平滑性好，并能保证在原始谱图下方；

2.本发明所设计的基于背景估计的非对称最小二乘基线校正方法操作简单，方便适用，算法在求解过程中，通过适当的数据变换，有效地减少了计算时间；

3.本发明所设计的基于背景估计的非对称最小二乘基线校正方法算法针对性强，在算法设计过程中充分挖掘和利用谱图自身的物理特性，使待估基线更接近真实基线，从而提高后续定性/定量分析模型的精度；

4.本发明所设计的基于背景估计的非对称最小二乘基线校正方法算法

适用面广，对于包含窄峰和宽峰的光谱信号，都能取得不错的校正效果。

附图说明

图1为本发明所设计的基于背景估计的非对称最小二乘基线校正方法流程图。

具体实施方式

下面结合附图对本发明做进一步的详细说明：

如图1所示，本发明提供的一种基于背景估计的非对称最小二乘基线校正方法，在非对称最小二乘的基础上加入一阶导数约束项，并利用直方图背景估计方法的计算值做为基线的初始值，通过迭代地求解，得到最终的基线，以及基线校正后的谱图，从而消除基线对后续定性/定量分析模型精度的影响。

下面首先给出改进的非对称最小二乘基线校正模型。

非对称最小二乘（Asymmetric Least Squares，AsLS）拟合思路来源于Whittaker平滑器，其优化目标如下：

$\arg \underset{z}{\min }[\underset{i}{\mathrm{\&Sigma;}}{w}_i{(y_i-z_i)}^1+\mathrm{\&lambda;}\underset{i}{\mathrm{\&Sigma;}}{\left({\mathrm{\&Delta;}}^2{z}_i\right)}^2]---\left(1\right)$

其中，w_i为权重因子，Δ为二阶差分算子，λ是正则化参数。式(1)中第一项表示拟合函数与原始数据的非对称拟合程度；第二项是为了保证拟合函数的光滑性，折中因子λ起到平衡非对称逼近程度和光滑性的作用。

最小化式(1)，可以导出下面的方程***：

(W+λD^TD)z=Wy                  (2)

式中，W是由向量w组成的对角阵，即W=diag(w)，D是z的二阶导数矩阵，即Dz=Δ²z。求解式(2)得到估计的基线：

z=(W+λD^TD)^-1Wy                (3)

式中权重系数w_i根据非对称的方式选择，当y_i>z_i时，w_i=p，而y_i≤z_i时，w_i=1-p。一般p取很小的值，其取值范围为0.001~0.1；这样对于y_i>z_i的点，其权重很小，而对于y_i≤z_i的点，其权重很大，这就是“非对称最小二乘”说法的来由。λ一般取很大的值，其范围为10²~10⁹。由于优化目标是凸函数，迭代过程会很快收敛。实际求解中，一般5到10次迭代即可收敛。

对于非对称最小二乘拟合方法，式(1)的约束项中只考虑了基线的平滑性约束。本文提出改进的非对称最小二乘基线校正算法，最小化下面的目标函数：

L(z)=||Q(y-z)||²+λ₁||D₁(y-z)²+λ||Dz||²               (4)

其中D₁是一阶微分矩阵，D是二阶微分矩阵，Q是一个对角权重矩阵。上式第一项反映红外光谱的基线拟合误差，第二项反映光谱的一阶导数拟合误差，第三项是光谱基线的光滑性约束。λ和λ₁是规整化参数，它在最小化拟合误差和保证基线平滑性之间起到折中作用。一般地，参数λ的取值范围为10~10⁴，参数λ₁取值小于10^-2即可。从上式可以看出，不仅要求拟合出的数值与原始数据之间的误差很小，而且还要求它们的一阶导数很接近。

权重矩阵Q由权值w构成，当y_i>z_i时，w_i=0，而y_i<z_i时，w_i=1。此时，式(4)的第一项也可以看成红外光谱的正性约束，它用来保证校正后的光谱y-z的负值部分的能量尽量小，而正值部分将不再考虑。

最小化目标函数L(z)最终转化为迭代求解线性方程***：

$(Q^{T}Q+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1+{\mathrm{\&lambda;D}}^{T}D)z=(Q^{T}Q+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1)y---\left(5\right)$

首先运用直方图背景估计方法估计出原始光谱y的背景值c⁽⁰⁾。基于背景值c⁽⁰⁾得到基线的初始估计值z⁽⁰⁾，并利用这个值来确定校正光谱y-z⁽⁰⁾为负值的位置，第一次迭代时，构造矩阵Q对这些负值部分进行规整化：

将构造的矩阵代入(5)式中，可以求解出新的基线的估计值z⁽¹⁾。

校正光谱的负值部分不断改变，这种迭代的基线求解算法可以表示为：

$\left\{\begin{array}{c}{Q}^{\left(k\right)}=\mathrm{diag}(y-z^{(k-1)}<0)\\ z^{\left(k\right)}={(Q^{\left(k\right)T}{Q}^{\left(k\right)}+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1+{\mathrm{\&lambda;D}}^{T}D)}^{-1}(Q^{\left(k\right)T}{Q}^{\left(k\right)}+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1)y\end{array}\right.---\left(7\right)$

给定初始基线，上述算法迭代5~10次就达到收敛。但是上述算法常无法充分地估计出光谱的基线，进而导致校正的光谱仍存在一定程度的基线漂移残留，因此实验中需要迭代地运用上述算法。

每次迭代时，对应的基线漂移残留体现在光谱背景上，因此在下一次迭代中，都利用新估计的基线更新光谱，并再次利用直方图背景估计法计算其背景值，然后重复基线求取的过程。如果连续两次迭代时，光谱的背景值几乎不改变，则整个算法终止。算法正是通过这种背景更新的方式来达到消除基线的目的。实际应用中，整个算法只需要迭代3~5次即可达到收敛。

基于以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于背景估计的非对称最小二乘基线校正方法，其特征在于，包括如下具体步骤：

步骤（1）：利用直方图背景估计方法计算出的背景值做为基线的初始值c_i；

步骤（2）：建立改进的非对称最小二乘基线校正模型；

步骤（3）：迭代地求解步骤（2）中的改进的非对称最小二乘基线校正模型，得到迭代后的光谱的背景值c⁽ⁱ⁺¹⁾；

步骤（4）：比较迭代前后光谱的背景值，如果几乎不改变，则整个算法终止，输出最终的基线。

2.根据权利要求1所述的基于背景估计的非对称最小二乘基线校正方法，其特征在于，所述步骤（1）包括如下具体处理：

步骤（11）：初始化i←0,X₀←{窗口中所有样本}；

步骤（12）：计算X_i的标准差σ_i；

步骤（13）：计算柱子数目n←round((max(X_i)-min(X_i))/σ_i)；

步骤（14）：计算柱子大小s←(max(X_i)-min(X_i))/n；

步骤（15）：计算直方图H←hist(X_i,n,s)；

步骤（16）：找出直方图中最大柱子的中心c_i；

步骤（17）：更新窗口，X_i+1←{x∈X_i|-2σ_i≤x-c_i≤2σ_i}；

步骤（18）：判断停止条件，如果|X_i+1|>0.95|X_i|，转步骤（110）；

步骤（19）：i＝i+1，返回步（2）；

步骤（110）：输出c_i。

3.根据权利要求1所述的基于背景估计的非对称最小二乘基线校正方法，其特征在于，所述步骤（2）包括如下具体处理：

步骤（21）：结合红外光谱自身的物理特性，找到谱图的变化特征，建立优化函数：

$\arg \underset{z}{\min }[\underset{i}{\mathrm{\&Sigma;}}{w}_i{(y_i-z_i)}^1+\mathrm{\&lambda;}\underset{i}{\mathrm{\&Sigma;}}{\left({\mathrm{\&Delta;}}^2{z}_i\right)}^2]$

其中，w_i定义为权重因子，Δ定义为二阶差分算子，λ定义为第一规整化参数；

步骤（22）：在步骤（21）的函数中加入一阶导数约束项，最小化目标函数如下：

L(z)=||Q(y-z)||²+λ₁||D₁(y-z)||²+λ||Dz||²

其中D₁定义为一阶微分矩阵，D定义为二阶微分矩阵，Q定义为一个对角权重矩阵，λ₁定义为第二规整化参数。

4.根据权利要求3所述的基于背景估计的非对称最小二乘基线校正方法，其特征在于，所述步骤（3）包括如下具体处理：

步骤（31）：将步骤（22）中的最小化目标函数L(z)最终转化为迭代求解线性方程***：

$(Q^{T}Q+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1+{\mathrm{\&lambda;D}}^{T}D)z=(Q^{T}Q+{\mathrm{\&lambda;}}_1{D}_1^T{D}_1)y;$

步骤（32）：运用直方图背景估计方法估计出原始光谱y的背景值c⁽⁰⁾，基于背景值c⁽⁰⁾得到基线的初始估计值z⁽⁰⁾，利用基线的初始估计值z⁽⁰⁾来确定校正光谱y-z⁽⁰⁾为负值的位置，第一次迭代时，构造如下矩阵：

将上述矩阵代入步骤（31）中的迭代求解线性方程***，求解出新的基线z⁽¹⁾；

步骤（33）：校正光谱的负值部分不断改变，这种迭代的基线求解算法可以表示为：

步骤（34）：重复步骤（31）至步骤（33）N次达到收敛，其中，N定义为自然数。

5.根据权利要求4所述的一种改进的非对称最小二乘基线校正方法，其特征在于，所述步骤（4）包括如下具体处理：

步骤（41）：利用直方图背景估计法，计算校正光谱y-z_(i)的背景值c⁽ⁱ⁺¹⁾，如果|c⁽ⁱ⁺¹⁾-c⁽ⁱ⁾|<ε，则算法终止，得到基线，以及基线校正后的谱图y。

Images