Raman spectroscopy-based method for identifying illegal food additives

Abstract

The application provides a food illegal additive identification method based on Raman spectrum enhancement, comprising collecting Raman spectrum of standard samples containing illegal additives and normal food matrix to construct a data set; a Raman characteristic dictionary is built, and an adaptive mask weighted residual matrix is generated according to the characteristic peak significance during iteration, the dictionary is iteratively optimized through singular value decomposition, and high-identification high-quality characteristic peaks are preferentially extracted; an optimization function is solved according to an alternating direction multiplier method, the spectrum to be measured is divided into three parts, i.e. a substance signal represented by the dictionary, an L1 constraint sparse noise and a low-rank fluorescent background, and initial sparse coefficients are obtained; the background component is matched with the matrix type, a corresponding quadratic penalty type interference model and a judgment threshold are called, and the coefficients are twice modified and optimized; if the atomic coefficients of the illegal additive corresponding to the dictionary exceed the threshold, it can be determined that the sample contains the illegal additive.

Description

Technical Field

[0001] This application belongs to the field of Raman spectroscopy enhancement, and in particular relates to a method for identifying illegal food additives based on Raman spectroscopy enhancement. Background Technology

[0002] Raman spectroscopy can collect the scattering spectral signals produced by specific molecules, reflecting the molecular characteristics of different substances and serving as a kind of molecular identification. However, actual food samples are usually not single pure substances, but complex mixtures of various organic and inorganic substances. When the target illegal additive exists in trace amounts in the food system, its Raman scattering signal is weak and easily masked by the strong signal of the food matrix itself, making it difficult to effectively extract and accurately identify the characteristic peaks of the target substance. In the actual detection and signal acquisition process, the acquired Raman spectral data not only contains the effective molecular signals of the mixture, but is also affected by a variety of complex interference factors. Because the Raman scattering effect itself is relatively weak, the spectral acquisition process is easily interfered with by external environmental factors and instrument detectors, resulting in randomly distributed spike noise in the spectrum. At the same time, organic macromolecules such as proteins and lipids commonly found in food produce strong fluorescence backgrounds under laser excitation. This type of low-frequency broadband fluorescence background usually causes the spectral baseline to rise or even be distorted, further obscuring the weak Raman characteristic peaks. Different types of food matrices exhibit varying degrees of interference and correlation effects on the spectral characteristics of illegal additives, further increasing the difficulty of signal analysis. Therefore, how to separate and extract highly identifiable feature information of target illegal additives under the combined interference of strong fluorescence background, random spike noise, and matrix signals, while reducing the impact of matrix interference and achieving qualitative and identification of illegal additives, is a pressing issue that needs to be addressed in the application of Raman spectroscopy in the field of food safety testing. Summary of the Invention

[0003] This invention proposes a method for identifying illegal food additives based on Raman spectroscopy enhancement. This method addresses the problem that existing technologies cannot separate and extract highly identifiable feature information of target illegal additives under the combined interference of strong fluorescence background, random spike noise, and matrix signals. The method includes:

[0004] The pure Raman spectra of the standard substance containing the target illegal additive and the normal food matrix spectra were obtained as training sample sets, as well as the Raman spectra data to be tested. A Raman feature dictionary is constructed. In the iteration of dictionary atoms, a feature saliency mask that evolves adaptively with the iteration is constructed. The mask weight is determined by the saliency of the feature peaks represented by the dictionary atoms. The residual matrix is ​​weighted using the feature saliency mask. The dictionary atoms are updated to convergence by performing singular value decomposition on the weighted residual matrix. This enables the dictionary to prioritize the identification of feature peaks with high signal-to-noise ratio and strong discriminative power. The Raman spectral data to be measured is sparsely decomposed, and the joint optimization objective function is solved. The Raman spectral data to be measured is decomposed into a linear superposition of three components, namely the mixed substance signal component linearly represented by the Raman feature dictionary, the sparse peak noise component subject to L1 norm constraint, and the low-rank fluorescence background component. The joint optimization objective function is solved by the alternating direction multiplier method to obtain the initial sparse representation coefficients and retain the low-rank fluorescence background component. The initial sparse representation coefficients are corrected, and the low-rank fluorescence background component is matched with a preset matrix background feature library to identify the matrix type. A specific interference model and decision threshold are called according to the matrix type. The interference model is a quadratic penalty matrix that describes the intensity of interference between spectra. The initial sparse representation coefficients are optimized and corrected by a quadratic penalty using the interference model. When the atomic coefficient of the illegal additive in the corrected coefficients exceeds the decision threshold, it is determined that the sample contains the illegal additive.

[0005] Furthermore, this invention also relates to a food additive identification system based on Raman spectroscopy enhancement, comprising the following modules: The acquisition module is used to acquire the pure Raman spectra of standard substances containing the target illegal additives and the spectra of normal food matrices as training sample sets, as well as the Raman spectra data to be tested; The decomposition module is used to construct a Raman feature dictionary. In the dictionary atom iteration, a feature saliency mask that evolves adaptively with the iteration is constructed. The mask weight is determined by the saliency of the feature peaks represented by the dictionary atom. The residual matrix is ​​weighted using the feature saliency mask. The weighted residual matrix is ​​then subjected to singular value decomposition to update the dictionary atom until convergence, so that the dictionary prioritizes the identification of feature peaks with high signal-to-noise ratio and strong discriminative power. The solution module is used to perform sparse decomposition on the Raman spectral data to be measured, solve the joint optimization objective function, decompose the Raman spectral data to be measured into a linear superposition of three components, namely the mixed substance signal component linearly represented by the Raman feature dictionary, the sparse spike noise component subject to L1 norm constraint, and the low-rank fluorescence background component. The joint optimization objective function is solved by the alternating direction multiplier method to obtain the initial sparse representation coefficients and retain the low-rank fluorescence background component. The identification module is used to correct the initial sparse representation coefficients, match the low-rank fluorescence background component with the preset matrix background feature library to identify the matrix type, and call a specific interference model and decision threshold according to the matrix type. The interference model is a quadratic penalty matrix that describes the intensity of interference between spectra. The interference model is used to perform a quadratic penalty optimization correction on the initial sparse representation coefficients. When the atomic coefficient of the illegal additive in the corrected coefficient exceeds the decision threshold, it is determined that the sample contains the illegal additive.

[0006] In the iterative process of constructing the Raman feature dictionary, this invention utilizes a feature saliency mask determined by the saliency of feature peaks to weight the residual matrix, making the dictionary update process more inclined to retain feature peak information with high signal-to-noise ratio and strong discriminative power. Through joint optimization decomposition, the spectral data to be measured is separated into three components: mixed substance signal, sparse peak noise, and low-rank fluorescence background, reducing the influence of noise and background baseline on the extraction of target feature peaks. Furthermore, the separated fluorescence background is used to match and identify the food matrix type, and based on the identification result, the corresponding quadratic penalty matrix interference model is invoked to perform secondary optimization and correction on the initial sparse representation coefficients. The above processing can weaken the spectral cross-interference between different food matrices, reduce the risk of false positives or false negatives, and, combined with a decision threshold, achieve the identification and judgment of illegal food additives. Attached Figure Description

[0007] Figure 1 A flowchart of the first embodiment; Figure 2 This is a schematic diagram of local residual singular value decomposition. Figure 3 This is a schematic diagram of matrix background matching. Detailed Implementation

[0008] To make the purpose, features, and advantages of this application more apparent and understandable, specific embodiments and accompanying drawings will be used to clearly and completely describe the technical solution protected by this application. Obviously, the embodiments described below are only some embodiments of this application, and not all embodiments. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0009] In the first embodiment, the present invention proposes a method for identifying illegal food additives based on Raman spectroscopy enhancement, such as... Figure 1 As shown, it includes: S1. Obtain the pure Raman spectra of the standard substance containing the target illegal additive and the normal food matrix spectra as training sample sets, and the Raman spectra data to be tested.

[0010] Raman spectra of the target illegal additive, a standard pure product, and normal food matrix samples without additives were collected using a portable Raman spectrometer. The acquired spectral data was converted into a processable spectral data matrix, and basic preprocessing was performed on the raw spectral data. During preprocessing, the spectral data were divided into feature training copies for dictionary learning and background retention copies for subsequent fluorescence background separation. For the feature training copies, the Asymmetric Least Squares Smoothing Baseline Correction (ALS) algorithm was used to reduce instrument baseline drift. For the Raman spectral data to be tested and the background retention copies used for matrix background matching, strong baseline subtraction processing that would remove low-frequency fluorescence background was not performed; only dark current correction, smoothing and denoising, and intensity normalization were performed, mapping the spectral intensity values ​​to the 0-1 range to retain the data foundation required for subsequent low-rank fluorescence background component separation and matrix matching. The combined processed standard pure product Raman spectra and normal food matrix Raman spectra were used to construct a training sample matrix, and the sample spectra collected on-site and subjected to the same preprocessing were used to construct a column vector of the Raman spectral data to be tested. Among them, the column vector of the Raman spectrum data to be measured retains the low-frequency fluorescence background components, which are used for subsequent low-rank fluorescence background separation and matrix type matching.

[0011] S2. Construct a Raman feature dictionary. In the iteration of dictionary atoms, construct a feature saliency mask that evolves adaptively with the iteration. The mask weight is determined by the saliency of the feature peaks represented by the dictionary atoms. Use the feature saliency mask to weight the residual matrix. Perform singular value decomposition on the weighted residual matrix to update the dictionary atoms until convergence, so that the dictionary can preferentially identify feature peaks with high signal-to-noise ratio and strong discriminative power.

[0012] An initial dictionary is constructed using the training sample matrix, and iterative optimization is performed based on the K-SVD framework. In each iteration, the sparse representation coefficient matrix under the current dictionary is calculated. For the k-th dictionary atom that needs updating, the residual matrix after removing the contributions of other dictionary atoms (excluding the k-th atom) is first calculated; that is, the contribution of the k-th dictionary atom and its corresponding sparse coefficients to be updated is retained in the current overall reconstruction result. Then, the training sample column with non-zero sparse coefficients in the k-th row is extracted to form a local residual matrix for updating the k-th dictionary atom. For the k-th dictionary atom to be updated, the saliency of the characteristic peak at each wavenumber position is calculated, generating a feature dimension weight vector with a length consistent with the number of spectral sampling points. This feature dimension weight vector is used as the main diagonal element to construct a feature saliency mask matrix; when the number of spectral sampling points is 1024, the feature saliency mask matrix is ​​a 1024×1024 dimensional diagonal matrix. Multiply the 1024×1024-dimensional feature saliency mask matrix by a local residual matrix of dimension 1024×M to achieve weighting of the wavenumber feature dimension, where M is the number of sample columns containing the non-zero sparse coefficients corresponding to the k-th dictionary atom.

[0013] The weighted local residual matrix is ​​subjected to truncated singular value decomposition (SVD). The first left singular vector obtained from the decomposition is extracted as the updated dictionary atom, and the product of the first right singular vector and the largest singular value is extracted to update the corresponding sparse representation coefficients. The feature saliency mask is only used to weight the local residual matrix in terms of feature dimension before SVD; it does not require inverting the mask matrix or performing an inverse mask mapping step. The above process of calculating the local residual and SVD is repeated, updating all dictionary atoms column by column, until the Frobenius norm of the residual matrix is ​​less than the preset convergence tolerance. The Raman feature dictionary that identifies high signal-to-noise ratio and highly discriminative feature peaks is then output. The singular value distribution obtained by SVD of the local residual matrix is ​​as follows: Figure 2 As shown. Since the singular vectors obtained by singular value decomposition have non-unique overall signs, after each dictionary atom update, the updated dictionary atom is correlated with the preset standard Raman spectrum set or the main feature peak template set, and the polarity of the dictionary atom is unified according to the correlation direction of the maximum correlation term; when the correlation coefficient of the maximum correlation term is negative, the dictionary atom and its corresponding sparse coefficient vector are multiplied by -1 to keep the activation direction corresponding to the target feature peak consistent.

[0014] After the Raman feature dictionary is output, the dictionary atoms are calibrated based on the pure Raman spectrum of the target illegal additive standard material. Specifically, the characteristic peak positions, peak shape templates, or normalized standard spectra in the pure Raman spectrum of the target illegal additive standard material are extracted, and the correlation between each dictionary atom and the characteristic peak template or normalized standard spectrum is compared. When the correlation between the dictionary atom and the target illegal additive standard spectrum reaches a preset atom calibration threshold, or when the sparsity coefficient of the dictionary atom in the target illegal additive standard sample is stably greater than a preset activation threshold, the dictionary atom is marked as a feature dictionary atom of the corresponding target illegal additive, and a mapping table between the target illegal additive and the dictionary atom index is formed. The atomic coefficients of the corresponding illegal additive mentioned in subsequent determination steps come from the dictionary atom index determined by this mapping table.

[0015] In an optional embodiment, constructing a feature saliency mask that adaptively evolves with iteration includes: For the k-th dictionary atom to be updated, the magnitude of the second derivative of the spectral intensity at the corresponding wavenumber after smoothing filtering is calculated as the significance of the characteristic peak. The normalized vector is obtained by using the range normalization method on the magnitude of the second derivative; A diagonal matrix is ​​constructed using the elements of the normalized vector as diagonal elements to serve as a feature saliency mask for the atom.

[0016] In the construction of the adaptive feature saliency mask, the k-th dictionary atom vector of length 1024 is received and smoothed using a filter with a window size of 11 and a polynomial order of 3 to reduce the interference of high-frequency white noise on spectral details. After filtering, the spectral intensity is calculated in the range of 200 to 3200. The second derivative within the wavenumber range is used to generate a 1024-dimensional amplitude vector by taking its absolute value. This second derivative amplitude can represent the sharpness and significance of Raman characteristic peaks. The larger the value, the more characteristic peaks of the illegal additives are in the region.

[0017] The calculated 1024-dimensional second-derivative magnitude vector is normalized using the Min-Max range, and a base bias of 0.1 is added to map each element to the interval [0.1, 1.0], resulting in a normalized vector. When the maximum and minimum values ​​are equal, the normalized vector is set to an all-1 vector or normalized after adding a preset minimum positive number to avoid a zero denominator. Using all elements of this 1024-dimensional normalized vector as the main diagonal elements and setting the remaining off-diagonal positions to zero, a 1024×1024-dimensional diagonal matrix is ​​constructed. This matrix serves as the feature saliency mask controlling the local iteration weights.

[0018] In an optional embodiment, the step of weighting the residual matrix using a feature saliency mask and then performing singular value decomposition on the weighted residual matrix to update the dictionary atoms includes: Calculate the residual matrix to be updated for the k-th dictionary atom, and extract the sample column containing the non-zero sparse coefficients corresponding to the k-th dictionary atom to form a local residual matrix. Use the corresponding feature saliency mask to weight the matrix to obtain a weighted local residual matrix. Singular value decomposition of the weighted local residual matrix yields the left singular matrix, the singular value matrix, and the right singular matrix. Update the first column vector of the left singular matrix with the dictionary atoms of the current iteration; The product of the first diagonal elements of the singular value matrix and the transpose of the first column vector of the right singular matrix is ​​updated to the sparse coefficient vector corresponding to the sample column containing the non-zero sparse coefficients in the current iteration. Based on the maximum correlation direction between the updated dictionary atoms and the preset standard Raman spectrum set or characteristic peak template set, synchronous polarity unification is performed on the updated dictionary atoms and their corresponding sparse coefficient vectors.

[0019] In the singular value decomposition update step of the dictionary atoms, the spectral residual dataset containing 500 samples is traversed. The global residual matrix after removing the contributions of all dictionary atoms except the k-th atom is calculated. That is, the k-th atom residual matrix to be updated is obtained by subtracting the sum of the products of the dictionary atoms other than the k-th atom and their sparse coefficients from the input data matrix. The column index set of all non-zero elements in the k-th row of the sparse coefficient matrix is ​​located. Based on this index set, the corresponding sample columns are extracted from the k-th atom residual matrix to be updated, and combined to obtain a local residual matrix of dimension 1024×M, where M is the number of non-zero sparse coefficients, for example, M=150. The previously generated 1024×1024-dimensional feature saliency mask diagonal matrix is ​​multiplied with the local residual matrix, that is, the feature saliency mask matrix is ​​left-multiplied by the local residual matrix, thereby achieving the relative prominence of the feature peak frequency band residual and the suppression of the flat region, resulting in a weighted local residual matrix of dimension 1024×M.

[0020] Singular value decomposition (SVD) is performed on the weighted local residual matrix, yielding a 1024×1024 orthogonal left singular matrix, a 1024×M diagonal singular value matrix, and an M×M orthogonal right singular matrix. When extracting and updating features, the first 1024-dimensional column vector of the left singular matrix is ​​truncated as the principal component, replacing the k-th dictionary atom in the current iteration to identify the Raman feature peak with the highest signal-to-noise ratio. The largest singular value scalar in the first row and first column of the singular value matrix is ​​extracted and multiplied by the transpose of the first M-dimensional column vector of the right singular matrix. The result updates the M-dimensional sparse coefficient vector corresponding to the sample column containing the non-zero sparse coefficients in the current iteration. Based on the maximum correlation direction between the updated dictionary atom and the preset standard Raman spectrum set or feature peak template set, synchronous polarity unification is performed on the dictionary atom and its corresponding sparse coefficient vector, thus completing the optimal update of the dictionary atom in a single iteration.

[0021] S3. The Raman spectral data to be measured is sparsely decomposed, and the joint optimization objective function is solved. The Raman spectral data to be measured is decomposed into a linear superposition of three components, namely the mixed substance signal component linearly represented by the Raman feature dictionary, the sparse peak noise component subject to L1 norm constraint, and the low-rank fluorescence background component. The joint optimization objective function is solved by the alternating direction multiplier method to obtain the initial sparse representation coefficients and retain the low-rank fluorescence background component.

[0022] The Raman spectral data to be measured is modeled as the product of the Raman feature dictionary and sparse representation coefficients, with equality constraints on the sparse peak noise component and the low-rank fluorescence background component, thus constructing a joint optimization objective function. The objective function includes L1 norm regularization terms for the sparse representation coefficients, L1 norm regularization terms for the sparse peak noise component, and nuclear norm regularization terms for the low-rank fluorescence background component. An augmented Lagrangian function is constructed using dual variables and a penalty parameter, and each component is alternately minimized. With other variables fixed, the sparse representation coefficients are updated using the Fast Iterative Thresholding (FISTA) algorithm combined with a soft thresholding operator.

[0023] The algorithm uses a soft-threshold shrinkage operator to process the current data residuals and update the sparse spike noise components. Singular value decomposition is performed on the residual matrix containing background information, and the singular values ​​are truncated using a soft threshold. The matrix is ​​then reconstructed to update the low-rank fluorescence background components. The dual variable is updated using the gradient ascent method. The above update steps are executed iteratively. At the end of each iteration, the changes in the original residual, dual residual, and the scalar value of the joint optimization objective function relative to the previous iteration are calculated. When both the original residual and dual residual are less than the preset residual tolerance, or when the change in the scalar value of the joint optimization objective function is less than 0.0001, the algorithm is considered converged. The converged sparse representation coefficients are extracted as the initial sparse representation coefficients, and the separated low-rank fluorescence background components are saved simultaneously.

[0024] In an optional embodiment, the step of decomposing the Raman spectral data to be measured into a linear superposition of three components includes: Establish an equality constraint equation that ensures the sum of the observed Raman spectral data to be measured is equal to the sum of the mixed material signal component, the sparse peak noise component, and the low-rank fluorescence background component, represented by the product of the Raman feature dictionary and the sparse representation coefficients. The L1 norm of the sparse representation coefficients corresponding to the mixed material signal components, the L1 norm of the sparse spike noise components, and the nuclear norm of the Hankel matrix mapped to the low-rank fluorescence background components are extracted to form a regularization penalty term. The augmented Lagrange multipliers corresponding to the equality constraint equation are introduced using a fixed penalty parameter to construct the augmented Lagrange joint optimization objective function.

[0025] An equality constraint equation of dimension 1024×1 is constructed in the vector space, where the measured Raman spectrum column vector is equal to the sum of the mixed material signal component, the sparse peak noise component, and the low-rank fluorescence background component. The mixed material signal component is obtained by matrix multiplication of a 1024×256 Raman feature dictionary matrix and a 256×1 dimensional sparse representation coefficient column vector. Using structured prior constraints, the L1 norm of the sparse representation coefficients is extracted to promote the sparse representation of the mixed material components, and the L1 norm of the sparse peak noise component is extracted to isolate sudden high-energy outliers such as cosmic rays. Simultaneously, for the smooth and continuous fluorescence background component, its one-dimensional data is reconstructed into a 512×513 Hankel matrix, and the nuclear norm of this matrix is ​​calculated to constrain the low-rank smoothing characteristics of the fluorescence background.

[0026] The 512×513 Hankel matrix is ​​constructed from a 1024-dimensional fluorescent background vector using an overlapping sliding window method, ensuring that elements on the same anti-diagonal line correspond to the same index in the original one-dimensional vector. During the inverse Hankel transformation, the 1024-dimensional fluorescent background column vector is restored by averaging the anti-diagonal elements. These three prior regularization terms are weighted and summed using fixed penalty parameters with preferred values ​​of 0.01, 0.1, and 1.5, respectively, to form the regularization penalty term. The penalty factor is then used... The quadratic penalty function is set to 0.5. Combined with a column vector of Lagrange multipliers of equal dimension 1024×1, the above equality constraints are transformed into an unconstrained augmented Lagrange joint optimization objective function through the augmented Lagrange method.

[0027] In an optional embodiment, the step of solving the joint optimization objective function using the alternating direction multiplier method to obtain the initial sparse representation coefficients while retaining the low-rank fluorescence background component includes: In each iteration step, the sparse representation coefficients, sparse spike noise components, low-rank fluorescence background components, and Lagrange multipliers are updated sequentially using an alternating optimization approach. The sparse spike noise component is updated using a soft threshold shrinkage operator; The singular value shrinkage operator is used to shrink the singular value decomposition matrix corresponding to the Hankel matrix mapped from the intermediate variables containing the observation spectral residuals, and the update is completed by the inverse Hankel transformation; Determine whether the change in the joint optimization objective function between the current iteration result and the previous iteration result satisfies the preset convergence condition. If the preset convergence condition is satisfied, terminate the iteration, output the sparse representation coefficients retained in the last iteration as the initial sparse representation coefficients, and retain the low-rank fluorescence background component.

[0028] In solving the augmented Lagrange joint optimization objective function, an iterative loop based on the alternating direction multiplier method was used, with a maximum of 500 iterations. In each iteration, the 256-dimensional sparse representation coefficients were updated using an iterative soft thresholding algorithm while keeping other variables fixed. A soft thresholding shrinkage operator with a threshold of 0.2 was used to truncate and update the spike noise component, removing transient spike noise. For updating the low-rank fluorescence background, the calculated observed spectral residuals were reorganized into a 512×513-dimensional intermediate Hankel matrix. After performing singular value decomposition on the matrix, a singular value shrinkage operator with a parameter of 3.0 was applied to truncate each singular value to non-negative values ​​by subtracting the threshold, and singular values ​​below the threshold were set to zero.

[0029] After the singular value matrix is ​​shrunk, it is multiplied again to generate a low-rank matrix. An inverse Hankel transform is used to calculate the arithmetic mean of all elements along the anti-diagonal direction, reconstructing it into a 1024-dimensional fluorescence background column vector. The Lagrange multiplier vector is updated synchronously with a step size of 0.5. After each parameter update, the iteration is terminated according to the aforementioned convergence condition. If the condition is met, the current sparse representation coefficients and the low-rank fluorescence background column vector are output. The 256-dimensional sparse representation coefficient column vector corresponding to the last iteration station is output as the initial sparse representation coefficients, while the 1024-dimensional low-rank fluorescence background column vector is retained.

[0030] S4. Correct the initial sparse representation coefficients, match the low-rank fluorescence background component with the preset matrix background feature library to identify the matrix type, and call a specific interference model and decision threshold according to the matrix type. The interference model is a quadratic penalty matrix that describes the intensity of interference between spectra. Use the interference model to perform a quadratic penalty optimization correction on the initial sparse representation coefficients. When the atomic coefficient of the illegal additive in the corrected coefficients exceeds the decision threshold, it is determined that the sample contains the illegal additive.

[0031] Raman spectra of various normal food matrices were pre-acquired, and low-rank fluorescence background curves corresponding to each normal food matrix were extracted using a low-rank background separation method. These low-rank fluorescence background curves were normalized and calibrated to form a preset matrix background feature library. The Pearson correlation coefficients between the separated low-rank fluorescence background components and the standard background curves of each matrix in the preset matrix background feature library were calculated. The category with the largest Pearson correlation coefficient value greater than a set confidence matching threshold was selected as the food matrix type. If the largest Pearson correlation coefficient did not reach the confidence matching threshold, the sample to be tested was classified into a preset general matrix type, and a general interference model and a general decision threshold were used for subsequent correction and judgment. Taking dairy product matrix as an example, the matching comparison effect between the fluorescence background curve extracted from the sample and the reference spectrum of the standard dairy product matrix is ​​as follows: Figure 3 As shown.

[0032] During the matrix model calibration phase, for each preset food matrix, the spectra of normal samples and calibration samples with known concentrations of the target illegal additive were collected under that matrix. This yielded the sparse representation coefficients, reconstruction residuals, and dictionary atom co-occurrence relationships for that matrix. Based on the peak overlap, co-occurrence false trigger frequency, and correlation statistics in the standard samples among different dictionary atoms within the same matrix, a crosstalk matrix corresponding to that matrix was constructed. The matrix type, crosstalk matrix, and decision threshold determined through calibration samples were bound and stored, forming a correspondence table between matrix type, interference model, and decision threshold. A corresponding quadratic penalty matrix and decision threshold were then loaded based on the matrix type. This quadratic penalty matrix describes the spectral crosstalk intensity between different feature dictionary atoms under the current matrix environment. During construction, the Pearson correlation coefficient between the characteristic peaks of each atomic spectrum in the Raman feature dictionary is first calculated, and the absolute value of the correlation coefficient is taken or non-negative truncation is performed to obtain a non-negative crosstalk matrix. This avoids the opposite effect of negative correlation coefficients producing reward co-occurrence under the constraint of non-negative coefficients. Then, a semi-positive definite quadratic penalty matrix is ​​constructed based on the non-negative crosstalk matrix. For example, the non-negative crosstalk matrix is ​​multiplied by its transpose and superimposed with a preset minimum positive multiple of the identity matrix to ensure that the quadratic penalty term is a non-negative penalty term.

[0033] A quadratic penalty optimization problem is constructed with the initial sparse representation coefficients as the optimization variable. The objective function includes a reconstruction fidelity term, a coefficient fidelity term, and a quadratic penalty term. The reconstruction fidelity term ensures that the corrected sparse representation coefficients can still reconstruct the mixed material signal components after subtracting sparse peak noise components and low-rank fluorescence background components from the measured Raman spectrum, together with the Raman feature dictionary. The coefficient fidelity term ensures that the corrected sparse representation coefficients do not deviate excessively from the initial sparse representation coefficients. The quadratic penalty term is obtained by multiplying the transpose of the coefficient vector by the positive semi-definite quadratic penalty matrix corresponding to the current matrix, and then multiplying by the coefficient vector. By adding the reconstruction fidelity and coefficient fidelity terms, the optimization result is prevented from degenerating into a zero vector. The constraint condition includes non-negativity of the coefficients. After solving, the sparse representation coefficient vector corrected by the quadratic penalty optimization is obtained. Based on the mapping table between the target illegal additive and the dictionary atom index, one or more feature dictionary atom coefficients corresponding to the target illegal additive are extracted from the corrected sparse representation coefficient vector. When multiple feature dictionary atoms are involved, their maximum value, weighted sum, or average value can be taken as the comprehensive atom coefficient of the target illegal additive. This coefficient is then compared with the decision threshold under the corresponding matrix type. If the coefficient value is greater than the decision threshold, the sample is determined to contain the illegal additive; otherwise, the sample is determined not to contain the illegal additive.

[0034] In an optional embodiment, the step of matching the low-rank fluorescence background component with a preset matrix background feature library to identify the matrix type includes: Extract all matrix standard background curve data contained in the preset matrix background feature library; The Pearson correlation coefficient between the retained low-rank fluorescence background component and each matrix standard background curve was calculated sequentially. The maximum value among all Pearson correlation coefficient values ​​is selected. When the maximum value reaches the preset information matching threshold, the matrix category set by the matrix standard background curve corresponding to the maximum value is taken as the matrix type of the sample to be tested. When the maximum value does not reach the preset matching threshold, the sample to be tested is set to a preset general matrix type.

[0035] In the matrix interference classification and identification stage, all matrix standard background curve data contained in the preset matrix background feature library are extracted. This feature library contains 20 typical food matrices, such as milk, juice, and edible oil, and their corresponding matrix standard background curve data. Each standard curve is a 1024-dimensional column vector that has been extracted and normalized from low-rank fluorescence background. The Pearson correlation coefficient values ​​of the currently retained 1024-dimensional low-rank fluorescence background components and the 20 matrix standard background curves in the feature library are calculated sequentially, and the ratio of the product of covariance and standard deviation is calculated. After the calculation process is completed, a correlation coefficient array containing 20 continuous floating-point values ​​is obtained, with its value range distributed between [-1.0, 1.0]. The maximum value of the Pearson correlation coefficient is selected. If the maximum value reaches the preset confidence matching threshold, such as 0.85, the category label of the matrix standard background curve corresponding to the maximum correlation coefficient, such as "dairy matrix", is used as the matrix type of the sample to be tested. If the maximum value does not reach the preset confidence matching threshold, the sample to be tested is set as a general matrix type, and the preset general interference model and general decision threshold are called.

[0036] In an optional embodiment, the step of using an interference model to perform secondary penalty optimization correction on the initial sparse representation coefficients includes: Based on the initial sparse representation coefficients, construct the independent variable vector and establish a quadratic penalty function that includes a fidelity term and a quadratic penalty matrix with semi-positive definite quadratic penalty matrix in the called interference model as weights; By solving for the first-order gradient vector of the quadratic penalty function, the independent variable vector is iteratively updated using the gradient descent method. The independent variable vector update step is executed 50 times in a loop, and the result of the last loop calculation is used as the corrected sparse representation coefficient.

[0037] In the quadratic penalty optimization stage for reducing specific matrix collinearity interference, a matching interference model is invoked based on the previously identified matrix type. This model is represented by a 256×256-dimensional positive semi-definite quadratic penalty matrix, where the elements quantize the spectral crosstalk and overlap intensity of different feature dictionary atoms in this matrix environment. Using the output 256-dimensional initial sparse representation coefficients as the starting point of the independent variable vector, a quadratic penalty function is established. This function consists of a fidelity term based on the mean square error of spectral reconstruction using dictionary multiplication, a coefficient fidelity term used to constrain the correction coefficients to be close to the initial sparse representation coefficients, and a quadratic penalty term with the aforementioned positive semi-definite quadratic penalty matrix combined with a preferred weight coefficient of 0.05 as the weight, linearly superimposed.

[0038] In the optimization process, by solving the first-order gradient vector of the quadratic penalty function, the first-order gradient vector equation composed of the reconstruction error gradient, coefficient deviation gradient, and penalty matrix penalty gradient is derived. Using gradient descent with a step size of 0.001, the variables are iteratively updated by subtracting the product of the step size and the first-order gradient vector from the current independent variable vector. After each iteration, coefficients less than zero are truncated to zero to satisfy the non-negativity constraint of the sparse representation coefficients. This independent variable vector update step is executed 50 times to improve the convergence stability of the coefficient correction process. The 256-dimensional result vector output from the last iteration is used as the corrected sparse representation coefficients to determine whether the atomic coefficients corresponding to illegal additives exceed the judgment threshold.

[0039] In a second embodiment, the present invention also relates to a food additive identification system based on Raman spectroscopy enhancement, comprising the following modules: The acquisition module is used to acquire the pure Raman spectra of standard substances containing the target illegal additives and the spectra of normal food matrices as training sample sets, as well as the Raman spectra data to be tested; The decomposition module is used to construct a Raman feature dictionary. In the dictionary atom iteration, a feature saliency mask that evolves adaptively with the iteration is constructed. The mask weight is determined by the saliency of the feature peaks represented by the dictionary atom. The residual matrix is ​​weighted using the feature saliency mask. The weighted residual matrix is ​​then subjected to singular value decomposition to update the dictionary atom until convergence, so that the dictionary prioritizes the identification of feature peaks with high signal-to-noise ratio and strong discriminative power. The solution module is used to perform sparse decomposition on the Raman spectral data to be measured, solve the joint optimization objective function, decompose the Raman spectral data to be measured into a linear superposition of three components, namely the mixed substance signal component linearly represented by the Raman feature dictionary, the sparse spike noise component subject to L1 norm constraint, and the low-rank fluorescence background component. The joint optimization objective function is solved by the alternating direction multiplier method to obtain the initial sparse representation coefficients and retain the low-rank fluorescence background component. The identification module is used to correct the initial sparse representation coefficients, match the low-rank fluorescence background component with the preset matrix background feature library to identify the matrix type, and call a specific interference model and decision threshold according to the matrix type. The interference model is a quadratic penalty matrix that describes the intensity of interference between spectra. The interference model is used to perform a quadratic penalty optimization correction on the initial sparse representation coefficients. When the atomic coefficient of the illegal additive in the corrected coefficient exceeds the decision threshold, it is determined that the sample contains the illegal additive.

[0040] In an optional embodiment, constructing a feature saliency mask that adaptively evolves with iteration includes: For the k-th dictionary atom to be updated, the magnitude of the second derivative of the spectral intensity at the corresponding wavenumber after smoothing filtering is calculated as the significance of the characteristic peak. The normalized vector is obtained by using the range normalization method on the magnitude of the second derivative; A diagonal matrix is ​​constructed using the elements of the normalized vector as diagonal elements to serve as a feature saliency mask for the atom.

[0041] In an optional embodiment, the step of weighting the residual matrix using a feature saliency mask and then performing singular value decomposition on the weighted residual matrix to update the dictionary atoms includes: Calculate the residual matrix to be updated for the k-th dictionary atom, and extract the sample column containing the non-zero sparse coefficients corresponding to the k-th dictionary atom to form a local residual matrix. Use the corresponding feature saliency mask to weight the matrix to obtain a weighted local residual matrix. Singular value decomposition of the weighted local residual matrix yields the left singular matrix, the singular value matrix, and the right singular matrix. Update the first column vector of the left singular matrix with the dictionary atoms of the current iteration; The product of the first diagonal elements of the singular value matrix and the transpose of the first column vector of the right singular matrix is ​​updated to the sparse coefficient vector corresponding to the sample column containing the non-zero sparse coefficients in the current iteration. Based on the maximum correlation direction between the updated dictionary atoms and the preset standard Raman spectrum set or characteristic peak template set, synchronous polarity unification is performed on the updated dictionary atoms and their corresponding sparse coefficient vectors.

[0042] In an optional embodiment, the step of decomposing the Raman spectral data to be measured into a linear superposition of three components includes: Establish an equality constraint equation that ensures the sum of the observed Raman spectral data to be measured is equal to the sum of the mixed material signal component, the sparse peak noise component, and the low-rank fluorescence background component, represented by the product of the Raman feature dictionary and the sparse representation coefficients. The L1 norm of the sparse representation coefficients corresponding to the mixed material signal components, the L1 norm of the sparse spike noise components, and the nuclear norm of the Hankel matrix mapped to the low-rank fluorescence background components are extracted to form a regularization penalty term. The augmented Lagrange multipliers corresponding to the equality constraint equation are introduced using a fixed penalty parameter to construct the augmented Lagrange joint optimization objective function.

[0043] In an optional embodiment, the step of solving the joint optimization objective function using the alternating direction multiplier method to obtain the initial sparse representation coefficients while retaining the low-rank fluorescence background component includes: In each iteration step, the sparse representation coefficients, sparse spike noise components, low-rank fluorescence background components, and Lagrange multipliers are updated sequentially using an alternating optimization approach. The sparse spike noise component is updated using a soft threshold shrinkage operator; The singular value shrinkage operator is used to shrink the singular value decomposition matrix corresponding to the Hankel matrix mapped from the intermediate variables containing the observation spectral residuals, and the update is completed by the inverse Hankel transformation; Determine whether the change in the joint optimization objective function between the current iteration result and the previous iteration result satisfies the preset convergence condition. If the preset convergence condition is satisfied, terminate the iteration, output the sparse representation coefficients retained in the last iteration as the initial sparse representation coefficients, and retain the low-rank fluorescence background component.

[0044] In an optional embodiment, the step of matching the low-rank fluorescence background component with a preset matrix background feature library to identify the matrix type includes: Extract all matrix standard background curve data contained in the preset matrix background feature library; The Pearson correlation coefficient between the retained low-rank fluorescence background component and each matrix standard background curve was calculated sequentially. The maximum value among all Pearson correlation coefficient values ​​is selected. When the maximum value reaches the preset information matching threshold, the matrix category set by the matrix standard background curve corresponding to the maximum value is taken as the matrix type of the sample to be tested. When the maximum value does not reach the preset matching threshold, the sample to be tested is set to a preset general matrix type.

[0045] In an optional embodiment, the step of using an interference model to perform secondary penalty optimization correction on the initial sparse representation coefficients includes: Based on the initial sparse representation coefficients, construct the independent variable vector and establish a quadratic penalty function that includes a fidelity term and a quadratic penalty matrix with semi-positive definite quadratic penalty matrix in the called interference model as weights; By solving for the first-order gradient vector of the quadratic penalty function, the independent variable vector is iteratively updated using the gradient descent method. The independent variable vector update step is executed 50 times in a loop, and the result of the last loop calculation is used as the corrected sparse representation coefficient.

[0046] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for identifying illegal food additives based on Raman spectroscopy enhancement, characterized in that, Includes the following steps: The pure Raman spectra of the standard substance containing the target illegal additive and the normal food matrix spectra were obtained as training sample sets, as well as the Raman spectra data to be tested. A Raman feature dictionary is constructed. In the iteration of dictionary atoms, a feature saliency mask that evolves adaptively with the iteration is constructed. The mask weight is determined by the saliency of the feature peaks represented by the dictionary atoms. The residual matrix is ​​weighted using the feature saliency mask. The dictionary atoms are updated to convergence by performing singular value decomposition on the weighted residual matrix. This enables the dictionary to prioritize the identification of feature peaks with high signal-to-noise ratio and strong discriminative power. The Raman spectral data to be measured is sparsely decomposed, and the joint optimization objective function is solved. The Raman spectral data to be measured is decomposed into a linear superposition of three components, namely the mixed substance signal component linearly represented by the Raman feature dictionary, the sparse peak noise component subject to L1 norm constraint, and the low-rank fluorescence background component. The joint optimization objective function is solved by the alternating direction multiplier method to obtain the initial sparse representation coefficients and retain the low-rank fluorescence background component. The initial sparse representation coefficients are corrected, and the low-rank fluorescence background component is matched with a preset matrix background feature library to identify the matrix type. A specific interference model and decision threshold are called according to the matrix type. The interference model is a quadratic penalty matrix that describes the intensity of interference between spectra. The initial sparse representation coefficients are optimized and corrected by a quadratic penalty using the interference model. When the atomic coefficient of the illegal additive in the corrected coefficients exceeds the decision threshold, it is determined that the sample contains the illegal additive.

2. The method according to claim 1, characterized in that, The construction of the feature saliency mask that adaptively evolves with iteration includes: For the k-th dictionary atom to be updated, the magnitude of the second derivative of the spectral intensity at the corresponding wavenumber after smoothing filtering is calculated as the significance of the characteristic peak. The normalized vector is obtained by using the range normalization method on the magnitude of the second derivative; A diagonal matrix is ​​constructed using the elements of the normalized vector as diagonal elements to serve as a feature saliency mask for the atom.

3. The method according to claim 1, characterized in that, The step of weighting the residual matrix using a feature saliency mask and then performing singular value decomposition on the weighted residual matrix to update the dictionary atoms includes: Calculate the residual matrix to be updated for the k-th dictionary atom, and extract the sample column containing the non-zero sparse coefficients corresponding to the k-th dictionary atom to form a local residual matrix. Use the corresponding feature saliency mask to weight the matrix to obtain a weighted local residual matrix. Singular value decomposition of the weighted local residual matrix yields the left singular matrix, the singular value matrix, and the right singular matrix. Update the first column vector of the left singular matrix with the dictionary atoms of the current iteration; The product of the first diagonal elements of the singular value matrix and the transpose of the first column vector of the right singular matrix is ​​updated to the sparse coefficient vector corresponding to the sample column containing the non-zero sparse coefficients in the current iteration. Based on the maximum correlation direction between the updated dictionary atoms and the preset standard Raman spectrum set or characteristic peak template set, synchronous polarity unification is performed on the updated dictionary atoms and their corresponding sparse coefficient vectors.

4. The method according to claim 2, characterized in that, The process of decomposing the Raman spectral data to be measured into a linear superposition of three components includes: Establish an equality constraint equation that ensures the sum of the observed Raman spectral data to be measured is equal to the sum of the mixed material signal component, the sparse peak noise component, and the low-rank fluorescence background component, represented by the product of the Raman feature dictionary and the sparse representation coefficients. The L1 norm of the sparse representation coefficients corresponding to the mixed material signal components, the L1 norm of the sparse spike noise components, and the nuclear norm of the Hankel matrix mapped to the low-rank fluorescence background components are extracted to form a regularization penalty term. The augmented Lagrange multipliers corresponding to the equality constraint equation are introduced using a fixed penalty parameter to construct the augmented Lagrange joint optimization objective function.

5. The method according to claim 1, characterized in that, The method of using alternating direction multipliers to solve the joint optimization objective function, obtaining the initial sparse representation coefficients while retaining the low-rank fluorescence background component, includes: In each iteration step, the sparse representation coefficients, sparse spike noise components, low-rank fluorescence background components, and Lagrange multipliers are updated sequentially using an alternating optimization approach. The sparse spike noise component is updated using a soft threshold shrinkage operator; The singular value shrinkage operator is used to shrink the singular value decomposition matrix corresponding to the Hankel matrix mapped from the intermediate variables containing the observation spectral residuals, and the update is completed by the inverse Hankel transformation. Determine whether the change in the joint optimization objective function between the current iteration result and the previous iteration result satisfies the preset convergence condition. If the preset convergence condition is satisfied, terminate the iteration, output the sparse representation coefficients retained in the last iteration as the initial sparse representation coefficients, and retain the low-rank fluorescence background component.

6. The method according to claim 1, characterized in that, The step of matching low-rank fluorescence background components with a preset matrix background feature library to identify matrix types includes: Extract all matrix standard background curve data contained in the preset matrix background feature library; The Pearson correlation coefficient between the retained low-rank fluorescence background component and each matrix standard background curve was calculated sequentially. The maximum value among all Pearson correlation coefficient values ​​is selected. When the maximum value reaches the preset information matching threshold, the matrix category set by the matrix standard background curve corresponding to the maximum value is taken as the matrix type of the sample to be tested. When the maximum value does not reach the preset matching threshold, the sample to be tested is set to a preset general matrix type.

7. The method according to claim 1 or 6, characterized in that, The method of using an interference model to perform secondary penalty optimization and correction on the initial sparse representation coefficients includes: Based on the initial sparse representation coefficients, construct the independent variable vector and establish a quadratic penalty function that includes a fidelity term and a quadratic penalty matrix with semi-positive definite quadratic penalty matrix in the called interference model as weights; By solving for the first-order gradient vector of the quadratic penalty function, the independent variable vector is iteratively updated using the gradient descent method. The independent variable vector update step is executed 50 times in a loop, and the result of the last loop calculation is used as the corrected sparse representation coefficient.

8. A food additive identification system based on Raman spectroscopy enhancement, characterized in that, Includes the following modules: The acquisition module is used to acquire the pure Raman spectra of standard substances containing the target illegal additives and the spectra of normal food matrices as training sample sets, as well as the Raman spectra data to be tested; The decomposition module is used to construct a Raman feature dictionary. In the iteration of dictionary atoms, a feature saliency mask that evolves adaptively with the iteration is constructed. The mask weight is determined by the saliency of the feature peaks represented by the dictionary atoms. The residual matrix is ​​weighted using the feature saliency mask. The weighted residual matrix is ​​then subjected to singular value decomposition to update the dictionary atoms until convergence, so that the dictionary prioritizes the identification of feature peaks with high signal-to-noise ratio and strong discriminativeness. The solution module is used to perform sparse decomposition on the Raman spectral data to be measured, solve the joint optimization objective function, decompose the Raman spectral data to be measured into a linear superposition of three components, namely the mixed substance signal component linearly represented by the Raman feature dictionary, the sparse spike noise component subject to L1 norm constraint, and the low-rank fluorescence background component. The joint optimization objective function is solved by the alternating direction multiplier method to obtain the initial sparse representation coefficients and retain the low-rank fluorescence background component. The identification module is used to correct the initial sparse representation coefficients, match the low-rank fluorescence background component with the preset matrix background feature library to identify the matrix type, and call a specific interference model and decision threshold according to the matrix type. The interference model is a quadratic penalty matrix that describes the intensity of interference between spectra. The interference model is used to perform a quadratic penalty optimization correction on the initial sparse representation coefficients. When the atomic coefficient of the illegal additive in the corrected coefficient exceeds the decision threshold, it is determined that the sample contains the illegal additive.

9. The system according to claim 8, characterized in that, The construction of the feature saliency mask that adaptively evolves with iteration includes: For the k-th dictionary atom to be updated, the magnitude of the second derivative of the spectral intensity at the corresponding wavenumber after smoothing filtering is calculated as the significance of the characteristic peak. The normalized vector is obtained by using the range normalization method on the magnitude of the second derivative; A diagonal matrix is ​​constructed using the elements of the normalized vector as diagonal elements to serve as a feature saliency mask for the atom.

10. The system according to claim 8, characterized in that, The step of weighting the residual matrix using a feature saliency mask and then performing singular value decomposition on the weighted residual matrix to update the dictionary atoms includes: Calculate the residual matrix to be updated for the k-th dictionary atom, and extract the sample column containing the non-zero sparse coefficients corresponding to the k-th dictionary atom to form a local residual matrix. Use the corresponding feature saliency mask to weight the matrix to obtain a weighted local residual matrix. Singular value decomposition of the weighted local residual matrix yields the left singular matrix, the singular value matrix, and the right singular matrix. Update the first column vector of the left singular matrix with the dictionary atoms of the current iteration; The product of the first diagonal elements of the singular value matrix and the transpose of the first column vector of the right singular matrix is ​​updated to the sparse coefficient vector corresponding to the sample column containing the non-zero sparse coefficients in the current iteration. Based on the maximum correlation direction between the updated dictionary atoms and the preset standard Raman spectrum set or characteristic peak template set, synchronous polarity unification is performed on the updated dictionary atoms and their corresponding sparse coefficient vectors.

Images