CN103134789A - Spectrum recover method based on Laplacian-Markov field - Google Patents

Abstract

The invention discloses a spectrum recover method based on the Laplacian-Markov field. The a spectrum recover method based on the Laplacian-Markov field comprises the following steps: conducting normalization processing on discretized spectrum to obtain normalized spectrum intensity; calculating a first derivative of the normalized spectrum intensity; calculating a neighborhood standard deviation of each normalized sectrum fi, finding out a maximum standard deviation and a minimum standard deviation, constructing a weighting matrix according to maximum standard deviation and the minimum standard deviation; and using the Laplacian-Markov field, and solving a raman spectrum deconvolution and split Bregman iteration method. The spectrum recover method based on the Laplacian-Markov field suppresses noise in places gentle in intensity change, adopts weak reatrain to store detail of the spectrum in places drastic in intensity change, achieves balance of recover resolution ratio and noise compress capability, introduces the split Bregman iteration method to overcome indifferentiable problem of conventional gradient descenting method, and provides technical support for application of the raman spectrum in quality detection, material analysis and the like.

Description

Spectrum recovery method based on Laplacian-Markov field

Technical Field

The invention belongs to the technical field of Raman spectrum data processing, and particularly relates to a method for recovering a degenerated Raman spectrum by utilizing a Laplace-Markov Laplacian-Markov) field, a deconvolution and splitting Bregman iteration method, which provides technical support for the Raman spectrum in the aspects of quality detection, material analysis and the like.

Background

In recent years, the research on raman spectroscopy technology is active, so that raman spectroscopy is widely applied to the fields of food quality monitoring, chemical rapid detection, material analysis, biomedical science and the like. However, during the measurement process of the spectrometer, the spectral data is often affected by factors such as natural broadening of spectral lines, collisional broadening, low-pass characteristics of detectors and circuits, and the like, which causes strong noise and sub-band overlapping phenomena of the spectral data. These noise interferences cause difficulties in the spectral application and can even lead to serious identification and estimation errors, which greatly limits the application of the spectrum. Therefore, it is essential to preprocess data and improve the resolution of the spectrum.

At present, there are many methods for recovering the degraded spectrum, and these methods can be mainly classified into two main categories: curve fitting and spectral deconvolution. The curve fitting method considers that the spectrum curve is formed by combining a Gaussian function and a Lorentzian function, namely, the measured spectrum can be obtained by fitting the functions. The method can well solve the problem of subband overlapping, namely, overlapping peaks can be fitted, but the method is extremely sensitive to noise. Spectral deconvolution methods, mainly to recover the degraded spectrum by different constraints, are High Order Statistics (HOS), Tikhonov regularization, homomorphic filtering, etc. In these methods, strong constraints can well suppress raman spectral noise (Tikhonov regularization), but too strong constraints can cause that severely degraded spectral details cannot be recovered. The weaker constraint can better recover the degraded details (higher order statistics), but the noise suppression capability is reduced. Therefore, such methods are difficult to balance in recovering spectral resolution and noise suppression capability.

Disclosure of Invention

The invention provides a spectrum recovery method based on a Laplacian-Markov field, which solves the problem that Raman spectrum data measured by a spectrometer is often influenced by random noise, instrument errors and the like to cause low resolution.

A spectrum recovery method based on a Laplacian-Markov field is carried out according to the following steps:

(1) for discretized spectrum

Normalized processing is carried out to obtain normalized spectral intensity f_iI is 1,2,3 … n, n is the total number of spectra;

(2) calculating normalized spectral intensity f_iFirst derivative f of_i′；

(3) Calculating each spectral point f_iNeighborhood of [ f ]_i-3，f_i-2，f_i-1，f_i，f_i+1，f_i+2，f_i+3]Standard deviation of (S)_iFrom which the maximum standard deviation S is found_maxAnd minimum standard deviation S_minConstructing a weighting matrix Q of n x n dimensions, the diagonal elements of the weighting matrix QSetting the other elements to be 0, wherein ln is a natural logarithm, and e is a natural constant;

(4) and (3) solving the Raman spectrum by adopting a splitting iterative method:

(41) solving parameters by an initialization splitting method: initial value f of the spectrum⁰F, the initial value b of the over-noise-removal compensation ⁰0, substitute variable initial value d ⁰0, and 0 is the iteration number k;

(42) updating calculation f^k+1And d^k+1：

$f^{k+1}=f^k+\mathrm{\Δt}(\▿E\left(\stackrel{\‾}{f}\right)),$ $\▿E\left(\stackrel{\‾}{f}\right)=H^T{Q}^{T}Q(H\stackrel{\‾}{f}-g)+\mathrm{\β}(d^k-D\stackrel{\‾}{f}-b^k)$

Where max is the maximum value, H is the instrument response function, H^TIs the transpose of H, Δ t is the time step, n × n dimensional difference matrix of the spectrum

Normalized spectral intensity seti ═ 1,2, …, n, regularization parameters λ ∈ (0,1) and β ∈ (0, 1);

(43) update calculation b^k+1＝b^k+((Df)^k+1-d^k+1)；

(44) Judgment of f^k+1Whether or not an iteration stop condition (| f) is satisfied^k+1-f^k‖/‖f^k|) > epsilon, epsilon is an iteration stop threshold, if not, then return to step (42), otherwise output raman spectrum f^k+1。

The instrument response function H is a Gaussian function or a Lorentzian function.

The technical effects of the invention are as follows: first, a method of spectral deconvolution with data weighting is proposed. For the raman spectrum data curve, the influence of interference on the data itself is not the same for all data, and the influence of noise is larger for the wave number with gentle intensity change and smaller for the wave number with severe intensity change. The data weighting adopts different constraints according to different intensity changes, and adopts stronger constraint at the place with gentle intensity change to restrain good noise; the details of the good spectrum are preserved by taking weaker constraints where the intensity changes are severe. This achieves a balance between the ability to restore resolution and noise suppression. Secondly, a Laplacian-Markov prior is applied in the solution of the spectrum deconvolution, and is a weak smoothness prior which can well keep the details of the spectrum, so that the Laplacian-Markov prior is first applied in the deconvolution of the Raman spectrum data. Thirdly, a splitting Bregman iteration method is introduced to solve the deconvolution model, the splitting Bregman iteration method overcomes the infinitesimal problem of the conventional gradient descent method, and compared with the conventional method, the method has extremely high algorithm solving speed and high practical value.

Drawings

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a flow chart of a Laplacian-Markov spectrum deconvolution algorithm;

FIG. 3 is a diagram of a mathematical model of spectral degradation;

FIG. 4 shows glucose (C)₁₅H₂₉NO₇) Actually measuring a Raman spectrum curve;

FIG. 5 shows glucose (C)₁₅H₂₉NO₇) Actually measuring a spectrum curve after normalization of the Raman spectrum curve;

FIG. 6 is a neighborhood standard deviation S output spectral curve for each spectral point;

FIG. 7 is a spectral curve of each spectral point normalized by the neighborhood standard deviation S;

FIG. 8 is a diagram of kernel functions of different shapes having the same full width at half maximum;

FIG. 9 shows glucose (C)₁₅H₂₉NO₇) Result graphs of different iteration times of the spectrum; wherein,

FIG. 9(a) shows glucose (C)₁₅H₂₉NO₇) Performing spectrum curve original drawing;

FIG. 9(b) shows glucose (C)₁₅H₂₉NO₇) 3 times of iteration result graphs of the spectral curve;

FIG. 9(C) shows glucose (C)₁₅H₂₉NO₇) 5 times of iteration result graphs of the spectrum curve;

FIG. 9(d) shows glucose (C)₁₅H₂₉NO₇) 7 times of iteration result graphs of the spectral curve;

FIG. 9(e) shows glucose (C)₁₅H₂₉NO₇) A spectrum curve 10 times iteration result graph;

FIG. 9(f) shows glucose (C)₁₅H₂₉NO₇) A spectrum curve 20 times iteration result graph;

FIG. 9(g) shows glucose (C)₁₅H₂₉NO₇) A spectrum curve 30 times iteration result graph;

FIG. 9(h) shows glucose (C)₁₅H₂₉NO₇) A spectrum curve 40 times of iteration result graphs;

FIG. 9(i) is glucose (C)₁₅H₂₉NO₇) A spectrum curve 50 times iteration result graph;

FIG. 9(j) is glucose (C)₁₅H₂₉NO₇) A spectrum curve 60 times iteration result graph;

FIG. 9(k) shows glucose (C)₁₅H₂₉NO₇) A spectrum curve 70 times of iteration result graphs;

FIG. 9(l) is glucose (C)₁₅H₂₉NO₇) A spectrum curve 80 times iteration result graph;

FIG. 9(m) shows glucose (C)₁₅H₂₉NO₇) A spectrum curve 90 times iteration result graph;

FIG. 9(n) is glucose (C)₁₅H₂₉NO₇) A spectrum curve 100 times iteration result graph;

FIG. 10 is a comparison of the spectrum recovery results of degraded glucose by different methods, wherein 10(a) is a graph of the actually measured spectrum of a glucose instrument, 10(b) is a graph of the experimental result of the TRSD method, and 10(c) is a graph of the experimental result of the HOS method, and 10(d) is a graph of the experimental result of the LM-QSD method;

Detailed Description

The present invention is described in further detail below with reference to the attached drawings and examples.

For the raman spectrum data curve, the influence of interference on the data itself is not the same for all data, and the influence of noise is larger for the wave number with gentle intensity change and smaller for the wave number with severe intensity change. If the adaptive processing is performed according to the difference of the intensity change, the balance of the resolution restoring capability and the noise suppression capability can be achieved. Therefore, the invention provides a spectrum deconvolution method with data weighting, wherein the data weighting adopts different constraints according to different intensity changes, and adopts stronger constraint at a place with mild intensity change to inhibit noise; the details of the good spectrum are preserved by taking weaker constraints where the intensity changes are severe. On the other hand, a Laplacian-Markov prior is a weak smoothness prior. The smoothness of the prior constrained Raman spectrum is adopted, the Laplacian-Markov prior is applied to deconvolution of Raman spectrum data for the first time, and an experimental result also proves that the method is very effective for recovering the degraded spectrum.

Most raman spectral data measured by raman spectroscopy can be described as a convolution of the true signal with the corresponding function of the measurement instrument, as shown in fig. 2.

g(v)＝Hf(v)+n(v),

Wherein, f (v) is a real signal, g (v) is an observed signal measured by a measuring instrument, and H is a response function of the measuring instrument, also called a fuzzy core, and mainly comprises an internal linear function and an instrument broadening function.

As shown in FIG. 3, FIG. 3 is glucose (C)₆H₁₂O₆) The measured Raman spectrum curve has the abscissa Wavenumber of the curve as the wave number, the ordinate Raman Intensity as the Raman spectrum Intensity, and the resolution ratio of the spectrum curve is 1cm^-1The curve is 2850 to 2900cm^-1The wavenumber range has severe sub-band overlapping phenomena and is affected by noise, and the spectrum can be used for testing the resolution recovery capability of the algorithm on the spectrum curve and the noise suppression capability.

Referring to fig. 1, the method of the present invention specifically comprises the following steps:

(1) degraded spectrum to discretization

Preprocessing is performed to find the values of the points where the spectral intensity is the greatest and least, denoted as f_maxAnd f_minNormalizing the spectra in order to adapt the algorithm to different spectral curves, i.e.

$f_i=\frac{{\hat{f}}_i-f_{\min }}{f_{\max }-f_{\min }}{f}_i\∈\left(\mathrm{0,1}\right),i\∈(\mathrm{0,1,2,3}...n)$

As shown in FIG. 4, the glucose concentration is 2850 to 2900cm^-1The maximum value of the wavelength intensity in the wavenumber range is 15, but the maximum intensity variation range is greatly different for different wavenumber ranges and different substances, so that the normalization process is required, and the result of the normalization process performed on fig. 4 is shown in fig. 5.

(2) Normalized spectral intensity f_iThe first derivative and the second derivative of (a),

f_i′＝(f_i+1-f_i)/2，

(3) and solving the weighting matrix Q. Q has dimension n, for each spectral point f_iFind its neighborhood (f)_i-3，f_i-2，f_i-1，f_i，f_i+1，f_i+2，f_i+3) Standard deviation of S. When i is 1, f_iIs taken as (f)₁，f₁，f₁，f₁，f₂，f₃，f₄) (ii) a When i is 2, f_iIs taken as (f)₁，f₁，f₁，f₂，f₃，f₄，f₅) (ii) a When i is 3, f_iIs taken as (f)₁，f₁，f₂，f₃，f₄，f₅，f₆). Taking the maximum value and the minimum value of S as S respectively_max、S_minThen the diagonal element Q of the weighting matrix Q_iiCan be expressed as:

$q_{\mathrm{ii}}=\ln [\frac{(e-1)(S-S_{\min })}{S_{\max }-S_{\min }}+1]$

wherein ln is a natural logarithm, and e is a natural constant, and the value is 2.7.

Let sigma_i＝q_iiThen Q can be expressed as:

the standard deviation S is solved for the spectrum data in fig. 5, the solving result is shown in fig. 6, and it can be seen from comparison between fig. 5 and fig. 6 that the intensity change of the spectrum data is more drastic, the larger the value of the corresponding point S is, the noise is in an oscillation change rule, so the standard deviation S value is smaller, and different weights can be given to the spectrum recovery and the denoising according to S according to the characteristic of the spectrum. Normalizing S to obtain q_iiThereby obtaining a matrix Q, the result of which is shown in fig. 7.

(4) Setting an instrument response function H, the instrument response function usually has two selectable functions, a Gaussian function and a Lorentz function, as shown in fig. 8, which are a Gaussian function curve (Gaussian shape) and a Lorentz function curve (Lorentz shape), respectively, and it can be seen that although the shapes of the two functions are different, the half-widths and the heights are the same, and glucose (C) is considered₆H₁₂O₆) The symmetry of the molecule, a gaussian function was chosen as the kernel function for this example, with the width of the narrowest peak as the width of this gaussian kernel.

(5) And setting a splitting method to solve parameters. Lambda is a regular parameter, the function of lambda lies in balancing a data item and smoothness prior, the value range is (0,1), beta is the regular parameter in the splitting method, the value range is (0,1), k is the iteration number, and d is used for replacing f_iThe' is a first derivative of the spectrum, b is an over-denoising compensation value, the value range is (0,1), and the actual value is an empirical value. Initialization of the spectral initial value f⁰F, the initial value b of the over-noise-removal compensation ⁰0, substitute variable initial value d⁰＝0。

(6) Solving for f using a split Bregman iteration method^k+1And d^k+1To f for^k+1Linear solution by iterative Gauss-Seidel, H^TIs the transpose of matrix H. The delta t is a step of time,

energy functional, Df, is the derivation of the matrix.

$f^{k+1}=f^k+\mathrm{\Δt}(\▿E\left(\stackrel{\‾}{f}\right)),$ $\▿E\left(\stackrel{\‾}{f}\right)=H^T{Q}^{T}Q(H\stackrel{\‾}{f}-g)+\mathrm{\β}(d^k-D\stackrel{\‾}{f}-b^k)$

Wherein the n × n dimensional difference matrix of the spectrum

Normalized spectral intensity set

i＝1,2,…,n。

If the difference is (Df)^k+1+b^kIf the | is smaller than the threshold lambda/beta, the noise point is considered to be the noise point, and the intensity of the noise point is replaced by 0; if it is larger than the threshold lambda/beta, it is considered as a minutia point of the spectrum, it is smoothness-constrained with a Laplacian-Markov prior. Then-update b, b^k+1＝b^k+((Df)^k+1-d^k+1) Judgment f^k+1Whether or not an iteration stop condition [ (iif) is satisfied^k+1-f^k‖/‖f^kII > epsilon, epsilon is a threshold value for stopping iteration, and the value range is (10)^-9~10^-7) If the condition is not met, the loop step (6) is carried out to continue solving, and if the condition is met, the Raman spectrum f is output^k+1。

Setting glucose (C)₆H₁₂O₆) Process parameters, initialize f first⁰＝g，b⁰＝0，d⁰When the value is 0, two regular parameters in the splitting method take values respectively as follows: α =0.05 and β = 0.1. And (3) carrying out iterative solution, wherein the intermediate output result of the solution is shown in fig. 9, and as the iteration times are increased, the curve changes less and less, and finally converges when 100 times are carried out.

In order to verify the effectiveness of the algorithm of the present invention, the present invention has been compared with a Tikhonov Regularized Spectrum Deconvolution (TRSD) (see J.Ottaway, J.H.Kalivas, and E.Andries, "Spectral Multivariate Calibration with wavelet reconstruction Using Variants of" applied. Spectroscopy.64, 1388-1395 (2010)), a method of Higher Order Statistics (HOS) (see J.Yuan, Z.Hu, G.Wang, and Z.Xu, "structured high-order statistical breakdown of Spectral data," Chin.Opt.Lett.3,552-555(2005), in short, the name of the method of the present invention is-ordered Spectral data), and the length of the spectrum can be calculated as a ratio of the length of the spectrum to the length of the spectrum, as shown in the longitudinal length chart, and the length of the LM can be calculated as a ratio of the length of the longitudinal spectrum, as shown in the chart of the longitudinal length of the curve, and the length of the LM can be calculated as a ratio of the length of the longitudinal spectrum, as shown in the chart 1; compared with the TRSD algorithm and the HOS algorithm, the LM-QSD algorithm is faster than the TRSD algorithm and the HOS algorithm, and is faster than the TRSD algorithm and the HOS algorithm by more than 4 times, which is mainly benefited by the splitting method introduced by the invention.

TABLE 1 CPU time consumption comparison units(s) for three algorithms

Spectral curve	TRSD	HOS	LM-QSD
				D(+)-Glucopyranose	15.365s	17.524s	4.365s
methyl	25.625s	20.152s	5.574s

It will be understood by those skilled in the art that the foregoing is only a preferred embodiment of the present invention, and is not intended to limit the invention, and that any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the present invention.

Claims (3)

1. A spectrum recovery method based on a Laplacian-Markov field is carried out according to the following steps:

(1) for discretized spectrum

Normalized processing is carried out to obtain normalized spectral intensity f_iI is 1,2,3 … n, n is the total number of spectra;

(2) calculating normalized spectral intensity f_iFirst derivative f of_i′；

(3) Calculating each spectral point f_iNeighborhood of f_i-3，f_i-2，f_i-1，f_i，f_i+1，f_i+2，f_i+3]Standard deviation of (S)_iFrom which the maximum standard deviation S is found_maxAnd minimum standard deviation S_minConstructing a weighting matrix Q of n x n dimensions, the diagonal elements of the weighting matrix Q

Setting the other elements to be 0, wherein ln is a natural logarithm, and e is a natural constant;

(4) and (3) solving the Raman spectrum by adopting a splitting iterative method:

(41) solving parameters by an initialization splitting method: initial value f of the spectrum⁰F, the initial value b of the over-noise-removal compensation⁰0, substitute variable initial value d⁰0, and 0 is the iteration number k;

(42) updating calculation f^k+1And d^k+1：

$f^{k+1}=f^k+\mathrm{\Δt}(\▿E\left(\stackrel{\‾}{f}\right)),$ $\▿E\left(\stackrel{\‾}{f}\right)=H^T{Q}^{T}Q(H\stackrel{\‾}{f}-g)+\mathrm{\β}(d^k-D\stackrel{\‾}{f}-b^k)$

Where max is the maximum value, H is the instrument response function, H^TIs the transpose of H, Δ t is the time step, n × n dimensional difference matrix of the spectrum

Normalized spectral intensity set

i ═ 1,2, …, n, regularization parameters λ ∈ (0,1) and β ∈ (0, 1);

(43) update calculation b^k+1＝b^k+((Df)^k+1-d^k+1)；

(44) Judgment of f^k+1Whether or not an iteration stop condition (| f) is satisfied^k+1-f^k‖/‖f^k|) > epsilon, epsilon is an iteration stop threshold, if not, then return to step (42), otherwise output raman spectrum f^k+1。

2. The method for spectrum recovery according to claim 1, wherein the instrument response function H is a gaussian function or a lorentzian function.

3. A spectrum recovery system based on Laplacian-Markov field comprises

A first module for discretizing the spectrum

Normalized processing is carried out to obtain normalized spectral intensity f_iI is 1,2,3 … n, n is the total number of spectra;

a second module for calculating the normalized spectral intensity f_iFirst derivative f of_i′；

A third module for calculating each spectral point f_iNeighborhood of [ f ]_i-3，f_i-2，f_i-1，f_i，f_i+1，f_i+2，f_i+3]Standard deviation of (S)_iFrom which the maximum standard deviation S is found_maxAnd minimum standard deviation S_minConstructing a weighting matrix Q of n x n dimensions, the diagonal elements of the weighting matrix Q

Setting the other elements to be 0, wherein ln is a natural logarithm, and e is a natural constant;

a fourth module for solving the raman spectrum using a split iteration method, comprising:

the first submodule is used for initializing the splitting method to solve parameters: initial value f of the spectrum⁰F, the initial value b of the over-noise-removal compensation⁰0, substitute variable initial value d⁰0, and 0 is the iteration number k;

a second submodule for updating the calculation f^k+1And d^k+1：

$f^{k+1}=f^k+\mathrm{\Δt}(\▿E\left(\stackrel{\‾}{f}\right)),$ $\▿E\left(\stackrel{\‾}{f}\right)=H^T{Q}^{T}Q(H\stackrel{\‾}{f}-g)+\mathrm{\β}(d^k-D\stackrel{\‾}{f}-b^k)$

Where max is the maximum value, H is the instrument response function, H^TIs the transpose of H, Δ t is the time step, n × n dimensional difference matrix of the spectrum

Normalized spectral intensity set

i ═ 1,2, …, n, regularization parameters λ ∈ (0,1) and β ∈ (0, 1);

a third submodule for updating the calculation b^k+1＝b^k+((Df)^k+1-d^k+1)；

A fourth sub-module for judging f^k+1Whether or not iteration is satisfiedStop condition (| f)^k+1-f^k‖/‖f^k|) > epsilon, epsilon is an iteration stop threshold, if not, then return to step (42), otherwise output raman spectrum f^k+1。

Images