CN103018194A - Asymmetric least square baseline correction method based on background estimation - Google Patents

Abstract

The invention discloses an asymmetric least square baseline correction method based on background estimation. The asymmetric least square baseline correction method comprises the following steps of: calculating a background value by utilizing a histogram background estimation method to be taken as an initial value of a baseline; establishing an improved asymmetric least square baseline correction model; solving a model in an iteration manner; and judging iteration spectral background values at two times, if the background value is not changed, stopping the whole algorithm, and outputting the final baseline. The asymmetric least square baseline correction method based on the background estimation provided by the invention has the advantages that the operation is simple, the correction method is convenient and applicable, high in reliability and high in pertinence, an infrared spectrogram peak position and peak shape are kept to be invariant, meanwhile, the baseline of a spectrogram is effectively eliminated, and a good effect is obtained.

Description

Asymmetric least square baseline correction method based on background estimating

Technical field

The present invention relates to infrared spectrum Chemical Measurement qualitative/quantitative analysis technical field, particularly relate to a kind of asymmetric least square baseline correction method based on background estimating.

Background technology

Can be subject to the impact of the external environments such as ambient temperature, humidity and instrument self working condition in the spectrogram gatherer process, the spectrogram that causes finally obtaining presents baseline wander in various degree.Baseline wander can affect the precision of qualitative analysis on the one hand, can cause on the other hand spectrum peak height and peak area to measure distortion, thereby affects the precision of quantitative test.In addition, even similar spectrum samples, its baseline also may be very different, thereby the mathematical model that causes setting up based on these spectrum samples lacks robustness.In order to set up stable, a reliable qualitative analysis or Quantitative Analysis Model, must carry out baseline correction.

Baseline correction has two kinds of manual synchronizing and automatic calibrations.Manual synchronizing generally is to choose unique point on the spectrogram by man-machine two mutual modes, then these unique points is fitted to a curve.Time-consuming, the effort of this method, and because the unique point randomness of artificial selection is larger, its reappearance is relatively poor.The people such as Schulze have carried out comprehensive summary [1] to baseline automatic calibration algorithm, mainly contain at present differential and filtering method [2], morphological method [3], interpolation fitting method [4], background estimating method [5 ~ 8], etc.But these methods have weak point.Differential method can be eliminated constant and linear drift, but has also amplified the high frequency noise in the spectrogram simultaneously; Morphological method can effectively be eliminated the baseline baseline wander, but owing to the peak width of each absorption peak in the infrared spectrum is inconsistent, the selection of its structural element is very thorny; In the interpolation fitting method, automatically choosing of interpolation knot needs to rely on artificial experience, and the baseline that different interpolating functions simulates also is not quite similar; Background estimating method such as SNIP[8], its method is fast and convenient, but the selection of iterations is a difficult problem, excessive too small all improper.

Asymmetric least square baseline correction algorithm is that the people such as Eilers put forward [9] in 2003, and it has strict mathematical theory basis, when processing spectrogram preferably effect is arranged, thereby has been subject to paying close attention to widely.The method is estimated baseline [10 ~ 11] based on Whittaker smoother, requires the given spectrum of asymmetricly match of baseline to be asked, and treats simultaneously and asks baseline to apply certain slickness constraint.But this algorithm has only been considered the smoothness constraint of second derivative, and in fact under the level and smooth as far as possible prerequisite of baseline, we not only require the numerical value and the error between the raw data that simulate very little, but also require their first order derivative very approaching.Therefore, it is to have circumscribedly that asymmetric least-squares algorithm is directly used in baseline correction, must take into full account the physical characteristics of infrared spectrum self.

[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

Summary of the invention

Technical matters to be solved by this invention provides a kind of simple to operate, convenient applicable, reliability is high, and is with strong points, when keeping infrared spectrum peak position and peak shape constant, effectively eliminated the asymmetric least square baseline correction method based on background estimating of the baseline of spectrogram.

In order to solve above-mentioned technical matters: the present invention has designed a kind of asymmetric least square baseline correction method based on background estimating, comprises following concrete steps:

Step (1): utilize background value that histogram background estimating method calculates as the initial value c of baseline _i

Step (2): set up improved asymmetric least square baseline correction model;

Step (3): the improved asymmetric least square baseline correction model in the solution procedure (2) iteratively obtains the background value c of the spectrum after the iteration ⁽ⁱ⁺¹⁾

Step (4): compare the background value of iteration front and back spectrum, if change hardly, then whole algorithm stops, and exports final baseline.

As a kind of optimization method of the present invention: described step (1) comprises following concrete processing:

Step (11): initialization i ← 0, X ₀← all samples in the window;

Step (12): calculate X _iStandard deviation sigma _i

Step (13): calculate pillar number n ← round ((max (X _i)-min (X _i))/σ _i);

Step (14): calculate pillar size s ← (max (X _i)-min (X _i))/n;

Step (15): compute histograms H ← hist (X _i, n, s);

Step (16): the center c that finds out maximum pillar in the histogram _i

Step (17): upgrade window, X _I+1← { x ∈ X _i|-2 σ _i≤ x-c _i≤ 2 σ _i;

Step (18): judge stop condition, if | X _I+1| 0.95|X _i|, turn step (110);

Step (19): i=i+1 returns the step (2);

Step (110): output c _i

As a kind of optimization method of the present invention: described step (2) comprises following concrete processing:

Step (21): in conjunction with the physical characteristics of infrared spectrum self, find the variation characteristic of spectrogram, set up majorized function:

$\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]$

Wherein, w _iBe defined as weight factor, Δ is defined as the second order difference operator, and lambda definition is the first regularization parameter;

Step (22): add first order derivative bound term in the function of step (21), it is as follows to minimize objective function:

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

D wherein ₁Be defined as the single order differential matrix, D is defined as Second differential matrix, and Q is defined as a diagonal angle weight matrix, λ ₁Be defined as the second regularization parameter.

As a kind of optimization method of the present invention: described step (3) comprises following concrete processing:

Step (31): the objective function L (z) that minimizes in the step (22) finally is converted into the iterative linear equation system:

$(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;$

Step (32): use histogram background estimating method to estimate the background value c of original spectrum y ⁽⁰⁾, based on background value c ⁽⁰⁾Obtain the initial estimate z of baseline ⁽⁰⁾, utilize the initial estimate z of baseline ⁽⁰⁾Determine to proofread and correct spectrum y-z ⁽⁰⁾Be the position of negative value, for the first time during iteration, be constructed as follows matrix:

With the iterative linear equation system in the above-mentioned matrix substitution step (31), find the solution the baseline z that makes new advances ⁽¹⁾

Step (33): the negative loop of proofreading and correct spectrum constantly changes, and the baseline derivation algorithm of this iteration can be expressed as:

$\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.;$

Step (34): repeating step (31) reaches convergence N time to step (33), and wherein, N is defined as natural number.

As a kind of optimization method of the present invention: described step (4) comprises following concrete processing:

Step (41): utilize histogram background estimating method, calculation correction spectrum y-z _(i)Background value c ⁽ⁱ⁺¹⁾,

If | c ⁽ⁱ⁺¹⁾-c ⁽ⁱ⁾|＜ε, then algorithm stops, and obtains baseline, and the spectrogram y after the baseline correction.

The present invention compared with prior art has following advantage:

1. the designed asymmetric least square baseline correction method gained baseline flatness based on background estimating of the present invention is good, and can guarantee below original spectrogram;

2. the designed asymmetric least square baseline correction method based on background estimating of the present invention is simple to operate, and convenient applicable, algorithm is in solution procedure, and the data transformation by suitable has reduced computing time effectively;

3. the designed asymmetric least square baseline correction method algorithm based on background estimating of the present invention is with strong points, in the algorithm design process, fully excavate and utilize the physical characteristics of spectrogram self, make baseline to be estimated more near true baseline, thereby improve the precision of follow-up qualitative/quantitative analytical model;

4. the designed asymmetric least square baseline correction method algorithm based on background estimating of the present invention

Widely applicable, for the spectral signal that comprises narrow peak and broad peak, can both obtain good calibration result.

Description of drawings

Fig. 1 is the designed asymmetric least square baseline correction method process flow diagram based on background estimating of the present invention.

Embodiment

The present invention is described in further detail below in conjunction with accompanying drawing:

As shown in Figure 1, a kind of asymmetric least square baseline correction method based on background estimating provided by the invention, basis at asymmetric least square adds first order derivative bound term, and utilize the calculated value of histogram background estimating method as the initial value of baseline, find the solution iteratively, obtain final baseline, and the spectrogram after the baseline correction, thereby baseline eliminated to the impact of follow-up qualitative/quantitative analytical model precision.

The improved asymmetric least square baseline correction model of following given first.

Asymmetric least square (Asymmetric Least Squares, AsLS) match thinking derives from the Whittaker smoother, and its optimization aim is as follows:

$\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)$

Wherein, w _iBe weight factor, Δ is the second order difference operator, and λ is regularization parameter.First asymmetric fitting degree that represents fitting function and raw data in the formula (1); Second is in order to guarantee the slickness of fitting function, and the compromise factor lambda plays the effect of balance asymmetric approximation degree and slickness.

Minimize formula (1), can derive following equation system:

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

In the formula, the diagonal matrix that W is comprised of vectorial w, i.e. W=diag (w), D is the matrix of second derivatives of z, i.e. the Dz=Δ ²Z.Find the solution the baseline that formula (2) obtains estimating:

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

Weight coefficient w in the formula _iSelect according to asymmetrical mode, work as y _iZ _iThe time, w _i=p, and y _i≤ z _iThe time, w _i=1-p.General p gets very little value, and its span is 0.001 ~ 0.1; Like this for y _iZ _iPoint, its weight is very little, and for y _i≤ z _iPoint, its weight is very large, the cause of saying that Here it is " asymmetric least square ".λ generally gets very large value, and its scope is 10 ²~ 10 ⁹Because optimization aim is convex function, iterative process is rapid convergence very.In actual the finding the solution, general 5 to 10 iteration can restrain.

For asymmetric least-square fitting approach, only considered the smoothness constraint of baseline in the bound term of formula (1).This paper proposes improved asymmetric least square baseline correction algorithm, minimizes following objective function:

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

D wherein ₁Be the single order differential matrix, D is Second differential matrix, and Q is a diagonal angle weight matrix.The baseline fitting error of first reflected infrared light spectrum of following formula, the first order derivative error of fitting of second reflection spectrum, the 3rd is the slickness constraint of spectrum baseline.λ and λ ₁Be regularization parameter, it plays the compromise effect between the baseline flatness minimizing error of fitting and guarantee.Usually, the span of parameter lambda is 10 ~ 10 ⁴, parameter lambda ₁Value is less than 10 ^-2Get final product.As can be seen from the above equation, not only require numerical value and the error between the raw data simulate very little, but also require their first order derivative very approaching.

Weight matrix Q is made of weight w, works as y _iZ _iThe time, w _i=0, and y _i＜z _iThe time, w _i=1.At this moment, first the positivity constraint that also can regard infrared spectrum as of formula (4), the energy of the negative loop of the spectrum y-z after it is used for guaranteeing to proofread and correct is as far as possible little, and will no longer consider on the occasion of part.

Minimize objective function L (z) and finally be converted into the iterative linear equation system:

$(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)$

At first use histogram background estimating method to estimate the background value c of original spectrum y ⁽⁰⁾Based on background value c ⁽⁰⁾Obtain the initial estimate z of baseline ⁽⁰⁾, and utilize this value to determine to proofread and correct spectrum y-z ⁽⁰⁾Be the position of negative value, for the first time during iteration, structural matrix Q carries out Regularization to these negative loops:

In matrix substitution (5) formula with structure, can find the solution the estimated value z of the baseline that makes new advances ⁽¹⁾

The negative loop of proofreading and correct spectrum constantly changes, and the baseline derivation algorithm of this iteration can be expressed as:

$\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)$

Given initial baseline, above-mentioned algorithm iteration just reaches convergence 5 ~ 10 times.But above-mentioned algorithm Chang Wufa estimates the baseline of spectrum fully, and then the spectrum that causes proofreading and correct still exists baseline wander to a certain degree residual, therefore needs to use iteratively above-mentioned algorithm in the experiment.

During each iteration, therefore corresponding residual being embodied on the spectral background of baseline wander in next iteration, is all utilized the new baseline renewal spectrum of estimating, and is again utilized histogram background estimating method to calculate its background value, then the process asked for of repetitive baseline.If during double iteration, the background value of spectrum changes hardly, and then whole algorithm stops.Algorithm reaches the purpose of eliminating baseline by the mode of this context update just.In the practical application, whole algorithm only needs iteration can reach convergence 3 ~ 5 times.

Based on above-described specific embodiment; purpose of the present invention, technical scheme and beneficial effect are further described; institute is understood that; the above only is specific embodiments of the invention; be not limited to the present invention; within the spirit and principles in the present invention all, any modification of making, be equal to replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (5)

1. the asymmetric least square baseline correction method based on background estimating is characterized in that, comprises following concrete steps:

Step (1): utilize background value that histogram background estimating method calculates as the initial value c of baseline _i

Step (2): set up improved asymmetric least square baseline correction model;

Step (3): the improved asymmetric least square baseline correction model in the solution procedure (2) iteratively obtains the background value c of the spectrum after the iteration ⁽ⁱ⁺¹⁾

Step (4): compare the background value of iteration front and back spectrum, if change hardly, then whole algorithm stops, and exports final baseline.

2. the asymmetric least square baseline correction method based on background estimating according to claim 1 is characterized in that described step (1) comprises following concrete processing:

Step (11): initialization i ← 0, X ₀← all samples in the window;

Step (12): calculate X _iStandard deviation sigma _i

Step (13): calculate pillar number n ← round ((max (X _i)-min (X _i))/σ _i);

Step (14): calculate pillar size s ← (max (X _i)-min (X _i))/n;

Step (15): compute histograms H ← hist (X _i, n, s);

Step (16): the center c that finds out maximum pillar in the histogram _i

Step (17): upgrade window, X _I+1← { x ∈ X _i|-2 σ _i≤ x-c _i≤ 2 σ _i;

Step (18): judge stop condition, if | X _I+1| 0.95|X _i|, turn step (110);

Step (19): i=i+1 returns the step (2);

Step (110): output c _i

3. the asymmetric least square baseline correction method based on background estimating according to claim 1 is characterized in that described step (2) comprises following concrete processing:

Step (21): in conjunction with the physical characteristics of infrared spectrum self, find the variation characteristic of spectrogram, set up majorized function:

$\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]$

Wherein, w _iBe defined as weight factor, Δ is defined as the second order difference operator, and lambda definition is the first regularization parameter;

Step (22): add first order derivative bound term in the function of step (21), it is as follows to minimize objective function:

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

D wherein ₁Be defined as the single order differential matrix, D is defined as Second differential matrix, and Q is defined as a diagonal angle weight matrix, λ ₁Be defined as the second regularization parameter.

4. the asymmetric least square baseline correction method based on background estimating according to claim 3 is characterized in that described step (3) comprises following concrete processing:

Step (31): the objective function L (z) that minimizes in the step (22) finally is converted into the iterative linear equation system:

$(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;$

Step (32): use histogram background estimating method to estimate the background value c of original spectrum y ⁽⁰⁾, based on background value c ⁽⁰⁾Obtain the initial estimate z of baseline ⁽⁰⁾, utilize the initial estimate z of baseline ⁽⁰⁾Determine to proofread and correct spectrum y-z ⁽⁰⁾Be the position of negative value, for the first time during iteration, be constructed as follows matrix:

With the iterative linear equation system in the above-mentioned matrix substitution step (31), find the solution the baseline z that makes new advances ⁽¹⁾

Step (33): the negative loop of proofreading and correct spectrum constantly changes, and the baseline derivation algorithm of this iteration can be expressed as:

Step (34): repeating step (31) reaches convergence N time to step (33), and wherein, N is defined as natural number.

5. a kind of improved asymmetric least square baseline correction method according to claim 4 is characterized in that described step (4) comprises following concrete processing:

Step (41): utilize histogram background estimating method, calculation correction spectrum y-z _(i)Background value c ⁽ⁱ⁺¹⁾If, | c ⁽ⁱ⁺¹⁾-c ⁽ⁱ⁾|＜ε, then algorithm stops, and obtains baseline, and the spectrogram y after the baseline correction.

Images