Abstract
A blind deconvolution algorithm with modified Tikhonov regularization is introduced. To improve the spectral resolution, spectral structure information is incorporated into regularization by using the adaptive term to distinguish the spectral structure from other regions. The proposed algorithm can effectively suppress Poisson noise as well as preserve the spectral structure and detailed information. Moreover, it becomes more robust with the change of the regularization parameter. Comparative results on simulated and real degraded Raman spectra are reported. The recovered Raman spectra can easily extract the spectral features and interpret the unknown chemical mixture.
© 2014 Chinese Laser Press
Raman spectra often suffer from common problems of bands overlapping and Poisson noise [1,2]. Raman spectra measured by a spectrophotometer can be mathematically modeled as
where and are the measured Raman spectrum and actual spectrum, and stands for the instrument function, which mainly collects the instrumental broadening. denotes the Poisson noise. The goal of blind spectral deconvolution is to seek the best estimates of and based on measured spectrum and prior information about the actual spectrum scene. Over the years, many deconvolution methods have been developed, such as high-order statistic [3], Fourier self-deconvolution (FSD) [4], the Wiener filtering method [5], maximum entropy deconvolution [6], and the semi-blind deconvolution method [7,8]. In [9], Katrasnik et al. extended the Richardson–Lucy (RL) method to the spectroscopy field. Nevertheless, owing to the ill-posed nature of inverse problems, this algorithm often lacks stability and uniqueness. The amount of noise will increase with the iteration process. And the algorithm must be stopped before the noise rises above a certain level. Tikhonov regularization (TR) was initially introduced as a regularizer for spectral processing in [10,11]. Since then it has been used extensively and with great success for inverse problems because of its ability to suppress noise.The maximum a posteriori (MAP) technique is a commonly used approach to estimate the actual spectrum given the measured spectrum . This technique maximizes the conditional probability of an actual spectrum given a certain measured spectrum. Based on Bayes’s rule, MAP can be converted to one likelihood probability multiplied by two priori probabilities. In this paper, the MAP framework is employed to construct the spectral deconvolution model. Considering the cases in which the data are contaminated by Poisson noise, the intensity of each pixel in the measured spectrum is a random variable that obeys an independent Poisson distribution. Hence the likelihood probability can be written as
where denotes the spectral length.To estimate and , iterative deconvolution algorithms can be used. One option is to use the RL method [9], which minimizes the following functional to maximize Eq. (2):
The RL algorithm does not converge to the solution because the noise is amplified after a small number of iterations. In order to get better convergence, the RL method with Tikhonov regularization [11] (TR-RL) was proposed, and the cost functional to be minimized is then
where . Symbol is the regularization parameter, which plays a very important role, controlling the constraint strength. If is too small, the noise will not be well suppressed; conversely, if it is too large, the spectral structure will be removed. This means that the regularization parameter needs to be adaptively adjusted according to spectral structure information. For better structure information preservation, in the TR term, we added a nonnegative monotonically decreasing function to control the regularization about the spectrum [12]: where is a constant between 0.01 and 1. Combining the adaptive term into the TR model, the modified Tikhonov regularization (MTR) model with spectral structure information is defined: which can distinguish three types spectral regions (flat, noise, and structure regions), but TR cannot distinguish those regions. The difference is shown in Fig. 1. Substituting the TR term in Eq. (4) with MTR in Eq. (6), we introduce the total functional:Minimizing Eq. (7) by the expectation maximization (EM) algorithm, the numerical iterative algorithm is then obtained,
where the superscript indexes the number of the iteration. The term is calculated by the previous iteration result .Summarizing, two good advantages of MTR-RL can be drawn: (i) for the flat region because is close to 1, which means a large TR strength is enforced to these points, and then the noise will be suppressed; (ii) for structure regions because is small and almost close to 0 and weakens the TR strength, so the spectral detail and structure will be preserved. In other words, MTR-RL has the ability to adaptively adjust the regularization strength according to the spectral structure information.
To quantitatively evaluate the deconvoluted spectra, three metrics, the normalized mean square error (NMSE) , the full width at half-maximum ratio (FWHMR) [9], and the noise suppression ratio (NSR) , are employed. and are the widths of the bands in the measured spectrum and deconvoluted spectrum , respectively. NSR was defined as the ratio of the total variation of the measured spectrum to the total variation of the deconvoluted spectrum , which is the noise attenuation measure. Among these three metrics, the NMSE requires the existence of a reference spectrum . Therefore, it can only be used in simulated experiments. FWHMR and NSR, which are nonreference measures of spectra, can also be used in experiments performed on real Raman spectra. It has been verified that the two metrics can reflect the width reduction and noise suppression [9]. The larger the values of FWHMR and NSR, the higher the spectral quality.
To demonstrate the effectiveness of the proposed method, we execute simulations with three test spectra under the Poisson noise process. The spectra come from the demo spectral library of Bruker Optics Incorporation. They were degraded by a Gaussian function with standard variance 12 [shown in Fig. 2(b)]. The degraded spectrum becomes much smoother and less resolved, in which the bands become wider and lower. For instance, it is difficult to distinguish the peak at from that at . Then the overlap spectrum was contaminated by Poisson noise [illustrated in Fig. 2(c)]. Three type regions can be defined according to the spectral gradient value, which is shown in Fig. 1. For comparison, the FSD, RL, and TR-RL methods were adopted.
Figures 2(d)–2(f) show a simulation example of restoration by RL, TR-RL, and MTR-RL. We chose and 200 iterations. Comparing Fig. 2(d) with Fig. 2(f), it can be found that MTR-RL obviously reduces the artifacts as well as better preserves the spectral structure. The effects of the regularization parameter were also investigated. For the TR-RL and MTR-RL methods, the NMSE of the deconvoluted spectrum versus varying regularization parameter on the simulated degraded spectrum of Fig. 2(c) is plotted in Fig. 3. The MTR-RL method appears more robust with the change of the regularization parameter at a wide range from to at intervals of 0.005. However, the change of regularization parameter has a great impact on the performance of TR-RL. In particular, the NMSE value increases dramatically when the regularization parameter becomes large.
Moreover, the NMSE of degraded spectra and the best deconvolution spectra by the three methods are compared in Table 1. The best deconvolution spectrum is selected to be the one with the lowest NMSE when the regularization parameter changes. It can be seen that the MTR-RL method has achieved the smallest NMSE among the three methods. The NMSE versus the number of iterations by the three methods on the Raman spectrum [methyl formate ()] is shown in Fig. 4, where the convergence is also highlighted.
Finally, the proposed algorithm was applied to real Raman spectra. Owing to photon-limited detection, Raman spectra are always corrupted by Poisson noise. We tested eight Raman spectra, which are downloaded from [13,14]. We use , starting with a 69 nm Gaussian function with a standard deviation of 2 for the initial instrument function. In order to save space, only two of them are illustrated. Figure 5(a) shows the 600 length spectrum of Cr:LisAF crystal from 300 to 900 nm. The two peaks at 631 and 661 nm overlap each other, and only two pinnacles can be distinguished. Figures 5(b) and 5(c) are the deconvoluted results by TR-RL and MTR-RL. The two peaks at 631 and 661 nm are split into four peaks at 580, 604, 631, and 663 nm, respectively. However, the TR-RL result only has three peaks. The instrument function is plotted in Fig. 5(d), whose width is about 69 nm.
Figure 6(a) shows the 690 length Raman spectrum of (D+)-glucopyranose [14] from 10 to . Figure 6(b) shows the deconvolution data. Here we choose . The peak at is split into two peaks at 395 and , respectively, and the peak at is separated into two peaks at 542 and , respectively. The spectral resolution has been improved considerably, and the distortion of the relative intensity can be revised to some degree. The same spectrum has been deconvoluted by [15], but the deconvoluted result seems rather noisy. Table 2 shows the FWHMR and NSR values of all the eight Raman spectra. It is seen that all the blind deconvolution methods raise the FWHMR and NSR, but the proposed method obtains the highest values.
This paper has presented a new spectral blind deconvolution method for Raman spectrum with Poisson noise. The proposed method can automatically balance the regularized strength between different spectral structures. Comparative results on simulated and real spectra show that MTR-RL can effectively reduce the Poisson noise and preserve the spectral structure information. Moreover, it becomes more robust with the regularization parameter, which makes the method more preferable in practical applications. The recovered Raman spectra can easily extract the spectral features and interpret the unknown chemical mixture. Although the application considered here is Raman spectroscopy, the method is more generally applicable to fluorescence data, and so on.
Nevertheless, there may still be room for improvement in our proposed method. In fact, the instrumental broadening usually varies with wavenumber in dispersive Raman spectrometers, because of the monochromator. Thus, the instrument function is not necessarily unique. We will examine this in our future work.
ACKNOWLEDGMENTS
This research was partially funded by the Project of the Program for New Century Excellent Talents in University (NCET-11-0654), the National Key Technology Research and Development Program (2013BAH72B01), the National Key Technology Research and Development Program (2013BAH18F02), the Scientific R & D Project of State Education Ministry and China Mobile (MCM20121061), the National Social Science Fund of China (14BGL131), and the New PhD Researcher Award from the Chinese Ministry of Education. The authors sincerely thank Prof. Lórenzf-Fonfría for helpful discussions and providing the source codes of the FSD method.
REFERENCES
1. J. J. Weber and D. D. Yavuz, “Broadband spectrum generation using continuous-wave Raman scattering,” Opt. Lett. 38, 2449–2451 (2013). [CrossRef]
2. K. Y. Song and H. J. Yoon, “Observation of narrowband intrinsic spectra of Brillouin dynamic gratings,” Opt. Lett. 35, 2958–2960 (2010). [CrossRef]
3. J. Yuan, Z. Hu, and J. Sun, “High-order cumulant-based blind deconvolution of Raman spectra,” Appl. Opt. 44, 7595–7601 (2005). [CrossRef]
4. J. K. Kauppinen, D. J. Moffatt, H. H. Mantsch, and D. G. Cameron, “Fourier self-deconvolution: a method for resolving intrinsically overlapped bands,” Appl. Spectrosc. 35, 271–276 (1981). [CrossRef]
5. C. W. Helstrom, “Image restoration by the method of least squares,” J. Opt. Soc. Am. 57, 297–303 (1967). [CrossRef]
6. V. A. Lórenz-Fonfría and E. Padrós, “Maximum entropy deconvolution of infrared spectra: use of a novel entropy expression without sign restriction,” Appl. Spectrosc. 59, 474–486 (2005). [CrossRef]
7. L. Yan, H. Liu, S. Zhong, and H. Fang, “Semi-blind spectral deconvolution with adaptive Tikhonov regularization,” Appl. Spectrosc. 66, 1334–1346 (2012). [CrossRef]
8. H. Liu, T. Zhang, L. Yan, H. Fang, and Y. Chang, “A MAP-based algorithm for spectroscopic semi-blind deconvolution,” Analyst 137, 3862–3873 (2012). [CrossRef]
9. J. Katrasnik, F. Pernu, and B. T. Likar, “Deconvolution in acousto-optical tunable filter spectrometry,” Appl. Spectrosc. 64, 1265–1273 (2010). [CrossRef]
10. F. Stout, J. H. Kalivas, and K. Héberger, “Wavelength selection for multivariate calibration using Tikhonov regularization,” Appl. Spectrosc. 61, 85–95 (2007). [CrossRef]
11. D. K. Buslov and N. A. Nikonenko, “Regularized method of spectral curve deconvolution,” Appl. Spectrosc. 51, 666–672 (1997). [CrossRef]
12. L. Z. Deng, L. Cao, and H. Zhu, “Spectral semi-blind deconvolution with hybrid regularization,” Infrared Phys. Technol. 64, 91–96 (2014). [CrossRef]
13. NASA, “Absorption spectral data of Cr:LiSAF crystal,” http://aesd.larc.nasa.gov/gl/laser/spectra/spectra.htm.
14. S. B. Engelson, “Raman spectral of (D+)-glucopyranose,” http://www.models.life.ku.dk/specarb.
15. J. Yuan, Z. Hu, G. Wang, and Z. Xu, “Constrained high-order statistical blind deconvolution of spectral data,” Chin. Opt. Lett. 3, 552–555 (2005).