Adaptive Operator-Based Spectral Deconvolution With the Levenberg-Marquardt Algorithm

Spectral distortion often occurs in spectral data due to the influence of the bandpass function of the spectrometer. Spectral deconvolution is an effective restoration method to solve this problem. Based on the theory of the maximum posteriori estimation, this paper transforms the spectral deconvolution problem into a multi-parameter optimization problem, and a novel spectral deconvolution method is proposed on the basis of Levenberg-Marquardt algorithm. Furthermore, a spectral adaptive operator is added to the method, which improves the effect of the regularization term. The proposed methods, Richardson-Lucy (R-L) method and Huber-Markov spectroscopic semi-blind deconvolution (HMSBD) method, are employed to deconvolute the white light-emitting diode (LED) spectra with two different color temperatures, respectively. The correction errors, root mean square errors, noise suppression ability, and the computation speed of above methods are compared. The experimental results prove the superiority of the proposed algorithm.


Introduction
Charge coupled device (CCD) spectrometer is a basic equipment for measuring and analyzing the material structure and composition by using the optical principle. It is widely applied in agriculture, astronomy, environmental detection, film industry, food safety, color measurement, and other fields [1][2][3][4][5]. However, since influenced by the bandpass function of the spectrometer, the measured spectrum usually produces spectral distortion. The bandpass function is mainly related to the incident slit of the spectrometer. In order to obtain high signal-to-noise ratio (SNR) signals, it is necessary to increase the luminous flux, thus the width of the slit must be enlarged to ensure sufficient light energy. However, this will also lead to the augment of the bandpass function of the spectrometer which decreases the spectral resolution. To overcome the shortcoming, the spectral deconvolution algorithm is an effective and commonly used correction method to obtain the accurate measurement result.
Spectral deconvolution is a classical problem in the field of the spectroscopic instrument. The development of the spectral deconvolution algorithm can be divided into two stages. In the first stage, only the effect of the bandpass function on the original spectrum is considered. The original spectrum is solved by using the mathematical relationship among the original spectrum, the measured spectrum, and the bandpass function. Typical methods include the Stearns and Sterns (S-S) method [6], improved S-S method [7], maximum entropy deconvolution [8], higher-order statistical method [9], and differential operator method [10]. However, there are some defects in these methods, the problem itself is ill-posed [11,12]. When dealing with such a problem, these methods will amplify the noise and other errors in the spectrum and generate the artificial noise. Since the CCD contains fewer pixels previously, the data interval of the measured spectrum is larger, and it is insensitive to noise, these methods can also be well applied. In recent years, due to the improvement of the CCD technology, the spectral data interval has been reduced and the resolution is increased which makes the above methods no longer applicable. Different from the first stage, the methods of the second stage take into account the influence of the error term on the measurement spectrum by utilizing regularization methods to suppress the spectrum, avoiding the amplification of the error term. For example, Eichstädt et al. extended Richardson-Lucy (R-L) method [13] applied in image correction to spectroscopy. The R-L method introduces a new stopping function, which stops the iteration before the function converges completely, thus playing the role of regularization and avoiding over-correction. However, the stopping condition will make the R-L algorithm stop iteration earlier, resulting in undercorrection. Considering spectral deconvolution as a maximum a posteriori (MAP) estimation problem, Liu et al. proposed a Huber-Markov spectroscopic semi-blind deconvolution (HMSBD) method [14] based on Huber-Markov priori [15], which has better protection for the high frequency effective signal. However, since this method is essentially a gradient descent method, the convergence of the algorithm is slow. Jin et al. proposed a bandwidth correction method [16] for the light-emitting diode (LED) spectrum based on the Levenberg-Marquart (L-M) algorithm [17,18]. By optimizing the parameters of the He-Zheng LED model [19], the noise [20][21][22] of the LED spectrum can be eliminated and the spectrum itself can be smoothed. However, this method can only be used for deconvolution of LED spectra. When the measured spectra are quite different from the He-Zheng model, the effect of this method is poor.
In this paper, a spectral deconvolution method based on the Levenberg-Marquardt algorithm is proposed, and an adaptive operator is introduced as the coefficient of the regularization term. To verify the validity and applicability of this method, the white LED spectra affected by the bandpass function are constructed. The proposed methods, R-L method, and HMSBD method, are used to deconvolute the measured spectra. The correction error and the root mean square error, the noise suppression ability and the computation speed of all three methods, are compared. The experimental results show the superiority of the proposed method.

Proposed model
Since the broadening effect of spectrometers, spectral distortion often occurs in spectral data. Most spectral data measured by spectrometers can be modeled mathematically by convolution of the real spectrum and the bandpass function, as well as the errors superimposed on the spectra which are the mainly noise produced during the spectral measurement. This process can generally be described as where ( ) M λ and ( ) S λ denote the real spectrum and measured spectrum, b represents the bandwidth function (also called line spread function or blur kernel), and ( ) N λ represents the noise in the measured spectrum. It is assumed that the bandpass functions at any wavelength position are identical in (1), but in practice, they are similar, but not identical. This is due to the different aberrations at different wavelength positions. Therefore, it is more appropriate to use the following matrix form to describe the process. = + M BS N (2) where N denotes the noise of spectrometer which mainly includes the dark noise, readout noise, optoelectronic noise, and fixed pattern noise. B is the bandwidth function matrix, which contains the bandwidth function of all wavelength positions and can be expressed as 11 12 1 where n is the number of pixels. Hence, M , S , and N in (2) can be represented as Applying Bayes formula and logarithmic function, (4) can be converted to (5): In (5), ( | ) p M S denotes the probability of the measured spectrum behaving as M which entirely depends on the distribution of the noise while the real spectrum S and the bandpass function matrix B are known. For spectrometers, many literatures [20][21][22] have proved that the noise in instruments is mainly Gaussian noise. Thus, ( ) p N obeys the Gaussian distribution and can be expressed as where σ denotes the standard deviation of the noise. Equation (6) can be expressed in the form of norms as where 1 C is the constant coefficient. In Bayesian inference, the last item ( ) p S in (5) expresses the prior knowledge about the spectrum S, which is the spectral constraint item. Hence, its choice doesn't need to be solely artificial and can be selected by the prior knowledge. It still plays the role of a regularization term which is introduced to solve the ill-posed problem and make the spectrum smooth as long as it is a proper distribution. Some conventional models such as Laplacian prior constrain, total variation prior, and Gauss-Markov prior regularize the corresponding ill-posed problem by forcing spatial smoothness on the spectrum. Gauss-Markov priori provides a statistical description of the spectrum. It focuses on the distribution of the pixels in the spectrum relative to their adjacent pixels, so it can correctly reflect the spectral correlation. Gauss-Markov priori can also be incorporated with the local statistical features of the spectrum as prior information into the iterative restoration process of the spectrum, maximizing the use of the prior information of the spectrum. Due to these advantages, in this paper, Gauss-Markov priori is selected as the prior probability, and the expression of ( ) p S can be written as where 2 C is the constant coefficient. Equations (7) and (8) are substituted into (5) and simplified, the deconvolution problem is transformed into a multi-parameter optimization problem, and MAP estimation is equivalent to minimizing the objective function: where 2 α S acts as the regularization term [23][24][25][26][27][28][29][30][31] and resembles the well-known Tikhonov regularization, and α is the regularization parameter. The regularization parameter α is usually unknown and needs to be determined by special methods according to the practical problem.
There are three methods to solve regularization parameter α most commonly, which include the general cross-validation (GCV) method [32][33], L-curve method [13,34], and mean square of estimation (MSE) minimum method [35,36]. The three methods determine the regularization parameter based on different specifications, and the value of the regularization parameter is also different. However, although the introduction of the regularization term can effectively avoid the amplification of error factors such as noise in the solving process of the inverse problem, there are still some problems. When the spectral intensity is weak, which means the value of S is small, the regularization term 2 α S is also small, which has little effect on spectral constraints and cannot suppress noise well. On the contrary, when the spectral intensity is strong, which means the value of S is lager, the regularization term 2 α S is also lager, which has a strong inhibition effect on the spectrum and often results in undercorrection.
For better noise suppression and preventing undercorrection, adding an adaptive operator that can control the regularization term is a good solution [23][24][25]. The adaptive operator should depend on noise and spectral properties and should be a non-negative function. When the spectral intensity is weak, the adaptive operator needs to increase the value of regularization term appropriately to suppress the noise, and when the spectral intensity is strong, it needs to decrease the value of regularization term appropriately to avoid undercorrection. To satisfy the requirements of the above analysis, an adaptive operator for this problem can be formulated as where 1 K and 2 K are constants, and the adaptive operator ( ) S β is a function of S. Thus, the regularization term with an adaptive operator is rewritten to 2 αβ S . After this improvement, 2 αβ S is larger than 2 α S when the spectral intensity is lower, that is, the value of S is smaller.  As can be seen from Fig. 1, the regularization term with an adaptive operator satisfies the above analysis requirements and can be adapted to different spectra by adjusting the parameters 1 K and 2 K . Thus, the objective equation is updated to

Optimization
The optimization method refers to the method of determining the value of some optional variables under certain constraints, so that the selected objective function can be optimized. Mathematically speaking, it is a method of finding the extremum, that is, under a set of constraints of equality or inequality, the objective function reaches the extremum, i.e., the maximum or the minimum.
The L-M algorithm is the most widely used non-linear least squares algorithm. It is an algorithm to get the maximum (minimum) value by gradient, which belongs to the hill climbing method. It has the advantages of both the gradient method and Newton's method, which means fast convergence speed and strong stability of the algorithm. When the damping coefficient is very small, the step size is equal to that of Newton's method, and when the damping factor is very large, the step size is equal to that of the gradient descent method. This makes the L-M algorithm deal with redundant parameter problems effectively, so that the probability of the objective function falling into the local minimum is greatly reduced.
Because of the superiority of this algorithm, in this paper, the L-M algorithm is selected to solve (11). Assuming that δ is the optimal increment of the S, the updated S can be expressed as follows through the Taylor series expansion: where J is the Jacobian matrix, and T ( ) δ δ O denotes the remainder term. By substituting (12) into (11), the (13) can be obtained: Since (13) only has the minimum value without the maximum value, the minimum value is obtained when the derivative of (13) to δ is zero. To obtain the value of δ , let the derivative of (13) be zero: To obtain the optimal increment δ through (14), the inverse matrix of T ( ) αβ + J J I needs to be calculated. However, the necessary and sufficient condition for a matrix to be an invertible matrix is that all its eigenvalues are not zero. Therefore,  (15) where λ is the damping factor. Therefore, a new estimate of S can be obtained, as shown below: where k represents the number of iterations. Through multiple iterations of (16), the estimated spectrum can achieve the required accuracy. It is particularly important for the iterative process to select the initial values of parameters and stopping rules. For such non-constrained least squares problems, a good initial value selection can accelerate the convergence process as well as reduce the running time of the program, and the suitable stopping rules can prevent over-fitting, thus improving the correction effect of the algorithm.
Generally speaking, the measured spectrum is chosen as the initial value of the estimated spectrum for iteration. This applies a certain prior knowledge, because the measured spectrum is closest to the true spectrum compared with other values. Setting it to the initial value can prevent the algorithm from falling into the local optimal solution and converge quickly. The iteration stopping rules of the algorithm generally include two kinds. The first one is that the iteration of the algorithm has reached a certain number of times, and the second one is that the increment of the iteration has been reduced to a certain extent. In general, iteration increment δ decreases with the number of iterations. Although continuing iteration may make the corrected spectrum closer to the real spectrum, it also takes more time and may result in over-fitting. Therefore, after reaching a certain accuracy, it is a general choice to stop iteration and output the corrected spectrum. In this paper, the iteration is stopped and the result is output only if 500 k > or 8 10 k δ − < is satisfied.
According to the aforementioned analysis, the detail of the proposed method is as follows: (1) Select the values of α , λ , 1 (3) Output the estimated spectrum k S .

Experiments and discussion
In order to verify the effectiveness of the spectral deconvolution method based on the L-M algorithm and adaptive regularization term (LMAR) and whether the LMAR method is better than the spectral deconvolution method based on the L-M algorithm and original regularization term (LMOR), a series of experiments were carried out. Two kinds of white LED spectra with different color temperatures were simulated as the original spectra. The original spectra were convolved with a bandpass function and the noise was added to simulate the measured spectra. It is noted that in order to simulate the situation as realistically as possible, the relevant parameters of USB4000 spectrometer (Ocean Optics) were used here. The normalized Hg-Ar light spectrum measured by USB4000 was used as the bandpass function. The level of simulated noise was the same as that of USB4000 under 25 ℃ ambient temperature and 3.8 ms integration time. The original spectra and measured spectra of white LEDs are shown in Figs. 2 and 3, respectively.
As can be seen from Figs. 2 and 3, the original spectra are changed to varying degrees due to the influence of the bandwidth function. The narrowband part of the original spectrum has larger deformation, while the deformation of the broadband part is smaller.

Parameter analysis
The selection of the regularization parameter and the adaptive operator parameter is an important part of the algorithm, which balances the original objective function and the regularization term. In order to better explain the role of the three parameters, the experiment of optimal parameters selection for the LED spectrum was designed. The root mean square errors (RMSEs) and its type A uncertainties of the original spectrum and the corrected spectra were calculated by changing the values of the parameters in turn, and the relationship between the parameters and RMSE was obtained. The RMSE was used to represent the deviation between the corrected spectrum and the original spectrum. The smaller the numerical value is, the better the correction effect of this method is. It can be calculated by (17): where n represents the total number of the spectrum, y denotes the original spectrum and ˆi y represents the measured spectrum.
Uncertainty indicates the reliability of the results. The smaller the uncertainty is, the higher the quality of the data is and the higher the credibility of the results is. The greater the uncertainty is, the lower the quality of the data is and the lower the credibility of the results is. The type A uncertainty at confidence level P=0.95 is derived from (18): where n represents the repetition times of experiments and when n=10, and t=2.26. x denotes the corrected spectrum, and x is the average of the corrected spectrum.
The parameters corresponding to the minimum RMSE are chosen as the optimal parameters.  From these figures, it can be seen that with an increase in the parameters, RMSE curves first decrease and then rise, all of which have the lowest points. A set of parameters of RMSE which can make the correction spectrum reach the minimum are 0.005 α = , 1 1200 K = , and 2 65000 K = , respectively. The proposed correction algorithm uses only the optimal α , the optimal α and 1 K , and the optimal α , 1 K , and 2 K for correction, respectively. The correction results of LED 1# are shown in Fig. 7. The same operations are done on LED 2# and the optimal parameters are 0.008 α = It can be seen from Figs. 7 and 8 that the proposed algorithm has a good correction effect. Compared with the measured spectra, the corrected spectra are closer to the original spectra, which proves the effectiveness of the algorithm. However, if only the regularization parameter α is used, it can be clearly seen that the spectrum is undercorrected in the position of strong spectral intensity, and the phenomenon of oscillation appears in the position of low spectral intensity, due to the weak effect of regularization on noise suppression. When the adaptive operator is introduced and suitable K 1 and K 2 are selected, the oscillation at low spectral intensity is weakened and the correction effect of the position with strong spectral intensity is improved, which proves the superiority of the adaptive operator.

Comparison with other algorithms
In the next step of the experiment, the R-L method, the HMSBD method, and the LMAR method were used to deconvolute the measured spectra. Each method selected the number of iterations based on the individual stopping rule and the correction results are shown in Figs As can be seen from Figs. 9 and 10, all three methods show excellent improvement of the white LEDs, which can correct deformation and suppress noise well. However, in restoring the spectrum at the position with larger deformation and suppressing noise at the flat position, the LMAR method is better than the R-L method and HMSBD method, whose corrected spectra are closer to the original spectrum. To further analyze the correction effect of several methods, the correction errors, RMSE, and noise suppression ability of these methods were calculated. The correction error was obtained by subtracting the correction spectrum from the original spectrum, which could most intuitively reflect the correction effect of various correction algorithms at different wavelength positions. The RMSE was used to represent the deviation between the corrected spectrum and the original spectrum, and the noise suppression ability denotes the noise suppression effect of the algorithm. Since the measured spectra contain both noise and deformation, RMSEs in flat regions (350 nm -400 nm) of the two spectra were calculated to characterize the noise suppression degree of the algorithm. Figures 11 and 12 show the correction errors of all three methods for two different white LED spectra. The RMSE and noise suppression ability of the three methods were calculated, as shown in Tables 1 and 2   As can be seen from Figs. 11 and 12, in the position of larger deformation, the LMAR method has better correction effect and its correction error is closer to zero. Compared with the R-L method and the HMSBD method, the correction error of the LMAR method is smoother, which means that its correction spectrum is smoother. As can be seen from Table 1, the RMSEs of the LMAR method and the HMSBD method are much smaller than that of the R-L method. When dealing with the spectrum of white LED 1#, the RMSE of the LMAR method is smaller, and while dealing with the spectrum of white LED 2#, the RMSE of the HMSBD method is relatively smaller. This is because the spectrum of white LED 1# is mainly the narrowband spectrum while that of white LED 2# is the broadband spectrum. The LMAR method has better correction effect for the narrowband spectrum, while the HMSBD method has better correction effect for the broadband spectrum, which can also be seen in Figs. 11 and 12. As can be seen from Table 2, the RMSE of LMAR method in flat regions of the two spectra is the smallest of the three methods, which means that its noise suppression effect is the best among the three methods.
In addition to the correction effect, the computation speed of the correction methods is also an important consideration particularly if a large number of spectra are to be processed. In this paper, various methods were implemented by using MATLAB 2016a. The central processing unit (CPU) of the computer was Intel (R) Core (TM) i5-7300HQ, and the memory was 8GB. The computation speed of each method is shown in Table 3. It can be seen that the computation speed of the HMSBD method is much slower than those of the other methods, since the essence of the HMSBD method is a gradient descent method, of which convergence rate is fairly slow. The LMAR method and R-L method have similar computation speed, but the computation speed of the LMAR method is faster than that of the R-L method, which suggests that the convergence speed of the LMAR method is faster. According to the above discussion, one can conclude that the LMAR method has the excellent correction effect and fast computation speed, which is suitable for current real-time spectral processing.

Conclusions
Because of the broadening effect of the spectrometer, spectral distortion often occurs in the spectral data. Spectral deconvolution is an effective method to obtain accurate spectral data. In this paper, a spectral deconvolution method based on the L-M algorithm is proposed, and an adaptive operator is introduced as the coefficient of the regularization term. In order to verify the effectiveness of the proposed method and the effect of the adaptive operator, the proposed method, R-L method, and MAP method are used to correct different LED spectra. The correction error, the root mean square error, the noise suppression ability, and the computation speed of all three methods are compared. The results show that the LMAR method has the excellent correction effect and fast computation speed, and the adaptive operator can effectively enhance the correction effect of the method.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.