CN106960091B - THz spectrum rapid nondestructive testing method for fresh meat freshness K value - Google Patents

THz spectrum rapid nondestructive testing method for fresh meat freshness K value Download PDF

Info

Publication number
CN106960091B
CN106960091B CN201710167209.0A CN201710167209A CN106960091B CN 106960091 B CN106960091 B CN 106960091B CN 201710167209 A CN201710167209 A CN 201710167209A CN 106960091 B CN106960091 B CN 106960091B
Authority
CN
China
Prior art keywords
value
thz
prediction
bpann
sample
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201710167209.0A
Other languages
Chinese (zh)
Other versions
CN106960091A (en
Inventor
赵茂程
齐亮
唐于维一
李忠
居荣华
赵婕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Forestry University
Original Assignee
Nanjing Forestry University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing Forestry University filed Critical Nanjing Forestry University
Priority to CN201710167209.0A priority Critical patent/CN106960091B/en
Publication of CN106960091A publication Critical patent/CN106960091A/en
Application granted granted Critical
Publication of CN106960091B publication Critical patent/CN106960091B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G01N21/35Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
    • G01N21/3563Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light for analysing solids; Preparation of samples therefor
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G01N21/35Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
    • G01N21/3581Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light using far infrared light; using Terahertz radiation
    • G01N21/3586Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light using far infrared light; using Terahertz radiation by Terahertz time domain spectroscopy [THz-TDS]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/02Neural networks
    • G06N3/04Architecture, e.g. interconnection topology
    • G06N3/045Combinations of networks
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/02Neural networks
    • G06N3/08Learning methods
    • G06N3/084Backpropagation, e.g. using gradient descent
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G01N21/35Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
    • G01N21/3563Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light for analysing solids; Preparation of samples therefor
    • G01N2021/3572Preparation of samples, e.g. salt matrices

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Biochemistry (AREA)
  • Molecular Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Software Systems (AREA)
  • Mathematical Physics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Toxicology (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

A THz spectrum rapid nondestructive testing method for fresh meat freshness K value comprises the following steps: first) meat collection and processing: uniformly cutting fresh meat into meat slices, and refrigerating; when sampling and cutting, the fat and connective tissue of the meat are avoided; II) measuring the THz spectral data of the meat: a THz attenuated total reflection ATR detection mode is adopted to rapidly detect a sample; when the spectrum is collected, the temperature and the humidity of the environment are kept stable; thirdly), rapidly detecting the freshness K value of the fresh meat by adopting a relation prediction model between the THz spectrum data and the fresh meat K value: the relation model is a principal component regression PCR prediction model or a back propagation neural network regression BPANN prediction model. And substituting the THz spectrum data obtained in the step two) into the relation prediction model, and quickly calculating the K value of the unknown sample. The method for sample pretreatment of the invention, which only needs to cut the sample to form the meat slices to be put into an instrument for measurement, is much simpler and faster than the existing method.

Description

THz spectrum rapid nondestructive testing method for fresh meat freshness K value
Technical Field
The technical scheme belongs to the technical field of fresh meat detection, and particularly relates to a THz spectrum rapid nondestructive detection method for fresh meat freshness K value.
Background
Meat is a major source of protein, fat, etc. intake for people. The degree of freshness of fresh meat is particularly important in food safety. Taking pork as an example, the pork is the main type of meat product consumption in China, the pork consumption accounts for 77% of the total meat consumption amount in China, how to control the pork quality in the production and sale link and ensure the food safety is highly concerned by production enterprises, consumers and national quality inspection departments. The quality of fresh meat such as pig, cattle, sheep, fish and the like comprises the elements such as nutrient components, flavor, tenderness, water retention, freshness and the like, wherein the freshness is an important aspect.
Freshness is usually measured by two methods, sensory evaluation and physicochemical analysis. The sensory evaluation method has large subjective factors and poor repeatability; the physical and chemical analysis method is accurate and reliable and has high repeatability. The freshness indexes measured by physical and chemical analysis comprise flesh color, microorganism indexes, volatile salt-based total nitrogen indexes, K value, pH value, polyamine compounds and the like, wherein the K value based on the ATP decomposition process is proved to be feasible and is increasingly emphasized.
Some substances in fresh meat deteriorate along with the increase of the preservation time, and Adenosine Triphosphate (ATP) in the slaughtered meat is automatically decomposed in the following process: adenosine Triphosphate (ATP) → Adenosine Diphosphate (ADP) → adenosine phosphate (AMP) → inosinic acid (HxR) → hypoxanthine (Hx).
According to the ATP decomposition process, the following index equation for measuring the K value is provided:
Figure BDA0001250160060000011
wherein ATP, ADP, AMP, IMP, HxR, Hx-the content of each substance in the sample
The more fresh the meat, the lower the K value, and vice versa the higher the K value. The K value can be generally obtained by High Performance Liquid Chromatography (HPLC). Although the physical and chemical analysis method is accurate and credible, the process is time-consuming and destructive to the tested sample (the method for detecting the physical and chemical parameters in a damage way is to chop meat into muddy meat and then further process and test), the requirements of quick and extensive measurement in production and circulation links cannot be met, and the determination of the K value of pork by using a quick nondestructive detection method is becoming an important research hotspot of scientific research work.
At present, nondestructive detection methods for freshness indexes of fresh meat include near infrared spectrum analysis, visible-near infrared hyperspectral image analysis, electronic nose technology and the like. The spectral analysis method establishes a mathematical model for nondestructive detection of pork freshness based on the principle that chemical substances related to fresh meat freshness have characteristic expression in visible and near infrared bands, and realizes rapid and accurate determination of the freshness. However, in the experimental process of the near-infrared spectrometer, the sample bin needs to be vacuumized, online rapid detection cannot be realized, the optical fiber type spectrometer does not need to be vacuumized, but optical fiber is used for collecting spectral information of a sample, the coupling optical fiber is fine, the detection area is limited, and one-time comprehensive collection of the surface of the sample cannot be realized. Although the visible-near infrared hyperspectral image analysis technology can realize comprehensive acquisition, a standard reflectivity calibration plate is needed to calibrate a spectrum acquisition system before the start of each test and in the long-time test process, so that the temperature drift of a camera photosensitive device and the influence of the brightness change of a light source on a detection system are prevented, and online continuous detection is not facilitated. The electronic nose technology uses a gas sensor to collect the putrefactive odor of meat, and the freshness of the meat is determined based on the principle that the putrefactive degree is deeper and the odor release degree is stronger.
Disclosure of Invention
The Terahertz waves (THz) are positioned between millimeter waves and infrared rays, belong to far infrared bands, and have the frequency within the range of 0.1 THz-10 THz (the THz waves with the frequency within the range of 0.1 THz-2 THz are selected in the scheme). The THz wave is energetically between electrons and photons. The diverse nature of THz waves allows many chemical molecules to exhibit molecular motion properties at the THz band that are less than those at other bands. Nucleotides (such as ADP and ATP) and related substances (such as IMP) belong to biological small molecules, the absorption in the terahertz waveband is mainly due to rotation and vibration of molecules or integral vibration of molecular groups, the terahertz waveband has a terahertz spectrum characteristic structure, and a plurality of absorption peaks exist in the THz waveband.
Therefore, the invention takes the K value as a research parameter, establishes a relation model between THz spectral data and the fresh meat freshness K value by exploring the THz spectral characteristic of the fresh meat freshness, realizes the quick nondestructive detection of the fresh meat K value, and provides a novel quick nondestructive detection method of the fresh meat freshness, and the specific technical scheme is as follows:
a THz spectrum rapid detection method for fresh meat freshness K value is characterized by comprising the following steps:
first) meat collection and processing: uniformly cutting fresh meat into meat slices, and refrigerating; when sampling and cutting, the fat and connective tissue of the meat are avoided;
II) measuring the THz spectral data of the meat: adopting THz Attenuated total reflection Attenuated TotalReflection, and quickly detecting a sample in an ATR detection mode; when the spectrum is collected, the temperature and the humidity of the environment are kept stable;
thirdly), rapidly detecting the freshness K value of the fresh meat by adopting a relation prediction model between the THz spectrum data and the fresh meat K value: the relation model is a principal component regression PCR prediction model or a back propagation neural network regression BPANN prediction model;
substituting the THz spectrum data obtained in the step two) into a relation prediction model, and quickly calculating the K value of the unknown sample;
in the second step):
sampling the THz spectrum of the fresh meat sample in a THz Attenuated Total Reflection (ATR) detection mode; sampling the THz spectrum of the fresh meat sample to obtain sample spectrum data for processing;
all modeling samples are randomly divided into a correction set and a prediction set; (number of samples in calibration set and PreThe ratio of the number of samples in the measurement set is about 2: 1) the spectral data in the correction set is used for establishing a prediction model; the spectral data of the prediction set is used to verify the accuracy of the prediction model. Xn×mIs particularly directed to the spectral data in the correction set.
The sample spectral data matrix is
Figure BDA0001250160060000021
In the formula, n is the number of samples in the correction set, m is the number of spectral wave points, and x is the ATR reflectivity of the samples;
firstly, compressing the spectrum information: the matrix Xn×mDecomposed into two matrices U and P, denoted as Un×a=Xn×mPm×a(ii) a In the formula, the matrix U represents data xij(i 1,2, …, n; j 1,2, …, m) in a new coordinate system; the column vector of the matrix P represents the linear transformation coefficient of the new coordinate system and the original coordinate system; the dimension a of U is less than the dimension m of X; taking data I of each row of the matrix Ui(i ═ 1,2, …, a) as input to the predictive model of step three);
in the third step):
A. the PCR prediction model comprises the following steps:
performing multiple linear regression on the matrix U to obtain a multiple linear equation: y is UB + E, namely a PCR prediction model, and is used for predicting the K value; in the formula, Y is the actually measured K value of the sample set, and E is a random error variable matrix obeying normal distribution;
the parameter B is obtained by the following method: in the modeling process of the PCR prediction model, randomly dividing a plurality of samples into a plurality of correction sets and a plurality of prediction sets; the spectral data of the correction set is used for establishing a prediction model; the spectral data of the prediction set is used for checking the accuracy of the prediction model; and (3) forming the actually measured K values of the correction set samples into a matrix Y to obtain the solution of the equation: b ═ U '(U' U)-1U 'Y, where U' is the transposed matrix of U. And after the parameter matrix B is obtained, the THz spectrum data of the prediction set samples can be substituted into the multivariate linear equation to obtain the prediction K value of the prediction set samples.
B. The modeling step of the BPANN prediction model comprises the following steps:
1) constructing a multilayer back propagation neural network:
the method comprises the steps that a multi-layer back propagation neural network BPANN is regarded as a nonlinear function, the input and the predicted value of the BPANN are respectively an independent variable and a dependent variable of the nonlinear function, and the BPANN expresses a function mapping relation from the independent variable to the dependent variable;
the topological structure of the BPANN consists of three layers of networks, namely an input layer, a hidden layer and an output layer, wherein each layer comprises a plurality of neurons, and the neurons of the adjacent layers are in one-way connection with weight values;
Iiis the input value of the input layer, i ═ 1,2, …, a; hjIs the output value of the hidden layer, j ═ 1,2, …, b; the output layer being a node
Figure BDA0001250160060000031
The value of (a) is a predicted value of BPANN, namely a predicted K value of fresh meat; hwijIs the connection weight between the neuron of the input layer and the neuron of the hidden layer, owjThe connection weight between the neuron of the hidden layer and the neuron of the output layer;
number of nodes in hidden layer
Figure BDA0001250160060000032
In the formula: a is the number of nodes of an input layer, c is a constant within the range of 1-10, and 1 is the number of nodes of an output layer;
2) training the model of step 1), comprising:
2.1) network initialization: determining a and b according to the input dimension and the output dimension; initializing hwij、owjHidden layer threshold thajThe output layer threshold thb, wherein the weight and the threshold are random numbers in the range of-1, and the given learning rate η is 0.01;
2.2) hidden layer output calculation: output of hidden layer
Figure BDA0001250160060000033
In the formula, the function f (x) is a hidden layer excitation function, and the expression of the function is:
Figure BDA0001250160060000034
2.3) output layer output calculation: predicted output of BPANN
Figure BDA0001250160060000035
2.4) error calculation: prediction error of BPANN
Figure BDA0001250160060000036
In the formula, y is a K value actually measured by an HPLC method;
2.5) weight updating: hwij=hwij+ηHj(1-Hj)Iiowje,i=1,2,…,a;j=1,2,…,b;owj=owj+ηHje,j=1,2,…,b;
2.6) threshold update: thaj=thaj+ηHj(1-Hj)owje,j=1,2,…,b;thb=thb+e;
2.7) judging whether the iteration of the algorithm is finished according to the iteration times and the final iteration error, and returning to the step 2.2) if the iteration is not finished.
The modeling step of the BPANN model further comprises the step 3) of optimizing the BPANN model to obtain a BPANN improved model, the method is to optimize the BPANN prediction model as a weak predictor to obtain a strong predictor, and the steps comprise:
3.1) initialization: n is the number of samples in the correction set, and the distribution weight of training data in the initial correction set is Dt(j) Determining a threshold phi of a prediction error of a weak predictor, determining a BPANN weak predictor structure according to the input and output dimensions of a sample, and initializing network parameters;
3.2) training weak predictors: when the t-th weak predictor is trained, the weak predictor is trained by using training data to obtain a predicted K value output by the weak predictor
Figure BDA0001250160060000041
Then calculating the prediction error
Figure BDA0001250160060000042
T is the number of weak predictors;
3.3) computing the weights of the Weak predictors
Figure BDA0001250160060000043
3.4) adjusting the distribution weight of the training data and normalizing:
Figure BDA0001250160060000044
in the formula (I), the compound is shown in the specification,
Figure BDA0001250160060000045
3.5) circulating the steps of 3.2-3.4) for T times to obtain T weak predictors, wherein the prediction model is ft. Superposing according to the updated weight to obtain a prediction model of the strong predictor
Figure BDA0001250160060000046
The phi value is the phi value when the performance of the model is best selected by adopting a leave-one cross verification method, namely, when the phi value is different in a certain range and is investigated by a fixed step length, the RMSECV value under the leave-one cross verification method in the modeling step is taken as the optimal threshold value when the RMSECV value is minimum;
the operation process of the leave-one-cross verification method is as follows: in n samples of the correction set, each sample is taken as a prediction set independently, and the rest n-1 samples are taken as the correction set; obtaining a predicted value of a sample of the prediction set through a modeling step by taking the correction set as a basis; obtaining the predicted values of all samples through n times of circulation;
obtaining the performance index RMSECV of a leave-one-cross verification method,
Figure BDA0001250160060000047
in the second step), the spectral frequency range is 0.1 THz-2 THz.
In the second step), preprocessing the sample spectrum data, and using the data obtained by preprocessing as the data of the matrix X; the method for preprocessing the sample spectrum data comprises the following steps:
first, preprocessing the spectral data by using a first derivative FD:
Figure BDA0001250160060000051
wherein x is the ATR reflectivity of the sample at the number of wave points i; when i is 1, its first order differential value is
Figure BDA0001250160060000052
… …, respectively; when i is m, its first order differential value is
Figure BDA0001250160060000053
And smoothing the spectrum by adopting an SG polynomial: selecting adjacent data points in the odd number of spectral data sequences, forming a window by the data points, and smoothing the central point of the window; in the smoothing process, fitting a polynomial to the data points in the window, and calculating the smoothed points according to the obtained polynomial; the expression of the data point after smoothing is obtained as follows:
Figure BDA0001250160060000054
wherein the width of the window is 2d +1, s is constant, and λjSmoothing the coefficient sequence by SG;
the fitting of the polynomial uses a least squares method.
The method for sample pretreatment of the invention, which only needs to cut the sample to form the meat slices to be put into an instrument for measurement, is much simpler and faster than the existing method.
Drawings
FIG. 1 is a schematic of an ATR accessory optical path;
FIG. 2-1 is a graph showing the trend of change in K value by HPLC;
FIG. 2-2 is a diagram showing 6 pieces of original THz spectra obtained by continuously measuring a meat sample 6 times;
FIGS. 2-3 are graphs of spectra after SVSRS treatment;
FIG. 3-1 is the original spectral line plot of the THz spectrum for 80 samples;
FIG. 3-2 is a spectrum plot of the THz spectrum after 15-point first derivative and SG smoothing;
FIG. 4-1 is a graph of the predicted performance variation of the PCR model at different differential widths;
FIG. 4-2 is a K-value prediction scatter plot of the PCR model;
FIG. 5-1 is a block diagram of the BPANN model;
FIG. 5-2 is an algorithmic flow chart of the BPANN model;
FIG. 5-3 is a graph of RMSECV versus differential width and number of hidden layer nodes for the BPANN model (the minimum value of RMSECV is marked with red circles);
FIGS. 5-4 are K-value prediction scatter plots of the BPANN model;
FIG. 6-1 is an algorithmic flow chart of a BPANN improved model;
FIG. 6-2 is a graph of RMSECV versus threshold φ in the BPANN improved model (threshold φ is 0.05-0.25);
FIG. 6-3 is a graph of RMSECV versus threshold φ in the BPANN improved model algorithm (threshold φ is 0.13-0.14);
FIGS. 6-4 are K-value prediction scatter plots of the BPANN improved model.
Detailed Description
The following explains the research, development and experimental processes of the technical scheme:
1. the THz spectrum itself is characterized by:
one significant feature of terahertz technology is safety. Compared to X-rays having photon energy of kilo-electron volts, terahertz radiation has photon energy lower than the bond energy of various chemical bonds, and thus it does not cause harmful ionization reactions. In addition, terahertz radiation cannot penetrate the skin of the human body because water strongly absorbs terahertz waves. Therefore, even if the terahertz radiation is intense, the influence on the human body is only localized on the surface layer of the skin, unlike the case where the microwave penetrates into the inside of the human body. The THz safety is beneficial to designing safe and convenient nondestructive testing equipment and protecting testing personnel from radiation hazards.
Another significant feature of terahertz technology is spectral resolvability. THz radiation in the far infrared band has relatively low photon energy compared to near and mid infrared radiation, but this band still contains abundant spectral information. Organic molecules, due to their transitions in rotation and vibration (including collective vibration), exhibit strong absorption and dispersion characteristics in this frequency band.
2. The embodiment of molecular substances related to the freshness K value of fresh meat in the THz wave band:
fresh meat freshness-K value is measured by the content ratio of nucleotides (such as ADP, ATP) and related substances (such as IMP), the nucleotides and the related substances belong to biological small molecules, the absorption in a terahertz wave band is mainly due to the rotation and vibration of molecules or the integral vibration of a molecular group, and the characteristic structure of the terahertz wave spectrum is clear. The purine polycrystal has a plurality of absorption peaks in a 0.2-2.5 THz wave band.
3. Selection of THz detection mode:
the biological molecules contain more atoms and have high density of collective vibration modes; the interaction between the atoms in the molecule causes the vibration of the biomolecule to have large dissonance, so the terahertz spectrum of the biomolecule is often unevenly broadened and overlapped.
This problem can be solved in two ways. One, narrow the broadening of the absorption peak by using a low-temperature cooling mode (such as liquid nitrogen cooling), separate each absorption peak, observe the absorption characteristics; in another method, the THz spectrum is obtained under a Normal temperature measurement method, and various preprocessing methods are used for the spectrum, such as Standard Normal transformation (SNV), Multiple Scattering Correction (MSC), derivatives (including first and second derivatives), Savitzky-golay (sg) polynomial smoothing, wavelet transformation, and the like. Noise information such as high-frequency random noise, baseline drift, light scattering and the like in the acquired original spectrum is reduced as much as possible, a linear or nonlinear mathematical model is used for analyzing the complex relation between the spectral line and the measured substance, and a reliable and stable THz spectral analysis model is established.
The latter method is selected in consideration of the purpose of realizing nondestructive rapid detection on the detected fresh meat sample, and the measurement is carried out in a normal temperature mode.
4. Selection of THz test sample
The epidermis, muscle and fat in fresh meat (e.g. pork) samples have different refractive indices and transmittances in the THz band. Different tissue components, namely fat meat and lean meat, in a slice of meat are placed in a THz spectral analysis system, in the detection range (0.1 to 4THz) of an ATR reflection mode, the absorption coefficient of the lean meat is larger than that of the fat meat, and the difference between the absorption coefficient and the fat meat is larger and larger as the frequency of the THz wave is increased. As can be seen, the spectral differences of THz waves are large for different tissue compositions of meat.
Generally, consumers prefer to purchase meat products with leaner meat, so the quality of lean meat determines the quality of the overall meat. The technical scheme selects a lean meat part in a sample as a detection object of the THz spectrum.
The method is further illustrated below with reference to the fresh pork example:
1. materials and methods
1.1 test materials
The test material was chilled pork fillet, purchased daily from a local supermarket, returned to the laboratory in the morning and evenly cut into 2.5cm x 0.5cm pieces. The fat and connective tissue of pork was avoided during sampling to prevent interference of these components with the THz assay results. And packaging the pork sample with a freshness protection package, numbering and placing in a freezer at 4 ℃ for storage to be tested. Continuously collecting meat samples for 8 days, numbering according to the date, and sequentially taking the number of 7 d-0 d. In the afternoon of 0d, THz spectrum acquisition and K value determination were completed for 8 meat samples. The test is repeated for 10 times, and the spectral data and the physicochemical value of 80 meat samples are obtained.
1.2 THz Spectrum Collection
The model of the THz detection device is TAS7500SP (advanced, japan), the THz detection device operates in a room temperature environment, the frequency resolution is 7.6GHz, the detection frequency range is 0.1THz to 4THz, 498 sampling frequency points (wave points) are provided in total, and the spectral line of the sample is obtained by 2048 times of automatic scanning and averaging.
THz has strong interaction with water, and the transmission depth of a meat sample rich in water is only a few hundred micrometers, so that the transmission mode is not suitable for nondestructive detection of freshness of meat products, and the reflection mode cannot be detected because the reflected wave is absorbed by water. The THz Attenuated Total Reflection (ATR) detection mode can overcome the strong absorption of THz waves by free water and bound water which are rich in a sample, so that the THz characteristic of the chemical substances with micron-sized thickness on the surface of the sample can be reflected in the spectrum of the THz Total reflection waves. The optical path of an ATR detection accessory is shown in fig. 1.
After each sample was removed from the freezer, the package was removed and the sample was placed flat on the surface of the ATR test window as shown in fig. 1. THz spectral data are collected for 3 times on the upper surface and the lower surface of the sample respectively, 6 parts of THz spectral data can be obtained from each sample, and the 6 parts of the THz spectral data are arithmetically averaged to be used as final THz spectral data of the sample. Because the THz spectrometer is sensitive to temperature and humidity, the temperature and the humidity in a laboratory are kept basically consistent during spectrum collection.
1.3K value determination
And immediately after the THz spectral data of the sample is collected, K value measurement is carried out on the same sample. And detecting the K value in the meat product by using an HPLC method. Separating the six chemical components by utilizing the different flow rates of the association products of ATP decomposition (ATP, ADP, AMP, IMP, HxR and Hx) in the mobile phase in the stationary phase, respectively measuring the contents of the components, and calculating the K value of the sample to be measured according to the formula 1.
The sample pretreatment process for K value detection is as follows: the sample to be tested was chopped into a puree, from which (2.00. + -. 0.05) g was taken into a 50mL centrifuge tube, 20mL of a cooled 10% perchloric acid solution was added, vortexed for 1min, centrifuged at 8000r/min for 10min, and the supernatant was removed. The analyte in the precipitate was re-extracted with 20mL of 5% perchloric acid solution, centrifuged at 8000r/min for 10min, and the supernatants were combined. Adjusting the pH value of the extracting solution to be about 6.0 by using 10mol/L NaOH solution, then continuously adjusting the pH value to be 6.0-6.4 by using 1.0mol/L NaOH solution, and then using ultrapure water to fix the volume to 50 mL. Filtering with 0.45 μm microporous membrane, and storing the filtrate at 4 deg.C.
The HPLC conditions were as follows: finnigan Surveyor liquid chromatograph (Sammerfo, USA), AQ-C18 chromatographic column (Sammerfo, USA), mobile phase 0.05mol/L K3PO4Buffer (pH 6.5), buffer was prepared with ultrapure waterThe sample amount was 1. mu.L, the flow rate was 200. mu.L/min, and the detection wavelength was 254 nm. The related products of ATP decomposition are quantified by an external standard method, and the measuring range is 0-0.5 mmol/L.
Details of the reagents used in the assay are as follows: 6 standard substances (purity is more than or equal to 99%, Sigma-Aldrich company, USA) are contained in the ATP-related substances (ATP, ADP, AMP, IMP, HxR and Hx); potassium phosphate, sodium hydroxide, perchloric acid (analytically pure, chemical reagents of national drug group limited, china); the test water was ultrapure water prepared by Millipore Academic.
2. Results and analysis
2.1 analysis of K value
The K values of 80 samples were obtained by HPLC and statistically analyzed, and the results are shown in FIG. 2-1. The meat samples of different test batches have the same storage time, and the K value is different due to individual difference. As can be seen from the average value of the samples, the freshness is continuously reduced with the increase of the storage days, and the K value of the samples is gradually increased, but the increase is not fixed. The increase amounts of 5d to 6d and 6d to 7d are relatively large, which shows that the quality of the meat sample is greatly changed in the storage period; sensory evaluation also confirmed that the meat sample was sticky to the touch and significantly sour in the 6 d-7 d storage period, and therefore it was presumed that the freshness of the meat sample was significantly reduced in this period. The K values of all samples in the test cover different freshness of pork, the coverage range is wide, and the mathematical model for predicting the K value by the THz spectrum has better robustness.
80 samples were randomly divided into 54 correction sets and 26 prediction sets, with a ratio of approximately 2: 1. The correction set is used for establishing a mathematical model of THz spectrum prediction K value; the prediction set is used to verify the accuracy of the model built to predict the unknown sample K values. As can be seen from table 1, the K value ranges of the correction set, the prediction set and the total sample set are substantially the same, and the mean and the standard deviation are not significantly different, so that the sample segmentation of the correction set and the prediction set is suitable.
TABLE 1 statistical information of K-values for correction set and prediction set
Figure BDA0001250160060000081
2.2 spectral pretreatment
The spectral quality evaluation method comprises the following steps: the same surface of the meat sample is continuously measured for 6 times under the same test condition to obtain 6 THz spectrums, theoretically, the 6 spectrums should be completely overlapped, but the 6 spectrums which are continuously and repeatedly measured cannot be completely overlapped due to the influence of noise of an instrument and test errors, and in order to evaluate the quality of the spectrums, n for continuously and repeatedly testing the same surface of the sample for multiple times can be calculatedsStandard Variance of the strip spectra (SVSRS), smaller SVSRS indicates better Spectral quality.
Figure BDA0001250160060000082
In the formula xijThe sample ith measurement of ATR reflectance at wave point j,
Figure BDA0001250160060000083
n at wave point jsThe average of the ATR reflectance was measured.
FIG. 2-2 shows the THz spectrum of a meat sample at 6 times in the range of 0.1THz to 4THz, and FIG. 2-3 shows the SVSRS corresponding to the spectrum, from which it can be seen that the SVSRS value at 2THz to 4THz is significantly increased and the repeatability is poor. Therefore, 0.1 THz-2 THz is selected as the spectral frequency band for modeling.
Fig. 3-1 is the THz spectral line of 80 samples in the 0.1 THz-2 THz region, and it can be seen that there is a great difference in the raw spectral intensity of different pork samples. The spectral differences include not only differences in sample composition, but also measurement errors, baseline drift and background noise. In order to eliminate the interference information, besides keeping the test environmental factors consistent as much as possible, the spectral data must be preprocessed to reduce or remove various interference factors.
The invention adopts the First Order Derivative (FD) to preprocess the spectral data, and the First Order Derivative can reduce the data deviation caused by baseline shift, drift and background interference, so that the spectral characteristic closely related to freshness becomes more obvious.
The first derivative FD used in the invention preprocesses the spectral data:
Figure BDA0001250160060000084
wherein x is the ATR reflectivity of the sample at the number of wave points i; when i is 1, its first order differential value is
Figure BDA0001250160060000085
… …, respectively; when i is m, its first order differential value is
Figure BDA0001250160060000086
Since the derivative calculation adds noise, the derivative pre-processing is followed by a Savitzky-Golay (SG) polynomial smoothing spectrum. In the derivation process, the selection of the differential width is important: if the difference width is too small, the noise is large, and the quality of the established model is influenced; if the differential width is too large, the smoothing is excessive and a large amount of detail information is lost. The SG polynomial smoothing selects an odd number of adjacent data points in the data sequence, forms a window from the data points, and smoothes the center point of the window. In smoothing, a polynomial is fitted to the data points in the window and the smoothed points are calculated from the resulting polynomial. The fitting of the polynomial uses a least squares method. The smoothed data point expression can be obtained by calculation as follows:
Figure BDA0001250160060000091
wherein the width of the window is 2d + 1; s is a constant whose value is related to the window width; lambda [ alpha ]jThe coefficient sequence is smoothed for SG.
Fig. 3-2 is a first derivative spectrum obtained after 15-point first derivative and SG smoothing.
2.3 pork K value prediction model
The content of ATP-related products related to K values changes in the pork decay process, and the biological molecules have sensitive spectral response to THz waves, and the THz spectrum can reflect the change of the content of the biological molecules. There is thus an indirect correlation between THz spectral data and pork freshness. The method uses Principal Component Regression (PCR), nonlinear algorithm-error back propagation Neural Network Regression (BPANN) and BPANN improved models to respectively verify the relevance.
The evaluation method of the detection model comprises the following steps:
the index for evaluating the quality of the model is the correlation coefficient (R) of the correction setC) Correction Set Root Mean Square Error (RMSEC), prediction Set correlation coefficient (R)P) And Prediction Set Root Mean Square Error (RMSEP). These parameters are defined as follows:
Figure BDA0001250160060000092
Figure BDA0001250160060000093
Figure BDA0001250160060000094
Figure BDA0001250160060000095
in the formula
Figure BDA0001250160060000096
-the predicted value of the ith sample in the sample set (including the correction set and the prediction set);
yi-the measured value of the ith sample in the sample set (including the calibration set and the prediction set);
nC-number of samples in the calibration set;
nP-predicting the number of samples of the set;
ym-mean values of the calibration or prediction set samples.
Coefficient of correlation RCAnd RPThe larger and/or the smaller the RMSEC and RMSEP, the better the predictive power of the model.
The spectral data was calculated and modeled using Matlab R2009b (Mathworks, usa) software.
2.3.1 principal component regression PCR prediction model
The THz spectral data of 80 samples can be represented in a matrix form:
Figure BDA0001250160060000101
wherein n is the number of samples, m is the number of spectral points, and x is the ATR reflectivity of the sample; xn×mSpectral data refers to data in the correction set.
Each spectrum in the test contains 250 wave points, the number of the wave points of the spectrum is far larger than that of the sample, and overfitting can occur if the spectrum is directly used for regression analysis, so that the prediction accuracy and the stability of the model are reduced. Meanwhile, the data of the spectral wave points have high collinearity and correlation, and the result of the regression model is also distorted. The spectral information can be compressed by Principal Component Analysis (PCA), and the original spectral data can be approximately reflected by scores of a few Principal components, so that the information redundancy and the correlation of the spectral data points are eliminated.
The PCA algorithm is a common and effective information compression method, and starts from a variable covariance matrix formed by a plurality of samples, and adopts a characteristic decomposition method (namely linear transformation of a coordinate system) to obtain a principal component with the maximum variance to replace an original variable, so as to obtain a new variable, namely a principal component score. The main components are independent and independent, and can eliminate the correlation and information redundancy existing in the original data. Therefore, PCA is widely applied in the fields of multivariate statistical analysis, spectrum information compression, extraction and the like.
Using principal component analysis, the original variable matrix X can be decomposed into two matrices U and P, denoted as Un×a=Xn×mPm×a
In the formula, matrixU is a score matrix (score matrix) representing data xij(i ═ 1,2, …, n;: j ═ 1,2, …, m) the vector position (i.e., coordinate value) in the new coordinate system (principal component coordinate system); the matrix P is a loading matrix (loading matrix), and the column vector of the matrix P represents the linear transformation coefficient of the new coordinate system and the original coordinate system; the dimension a of U is less than the dimension m of X; thus, the high-dimensional spectral data space is reduced to the low-dimensional linearly independent principal component feature space. Taking data I of each row of the matrix Ui(i ═ 1,2, …, a) as an input to a predictive model;
PCR prediction model: performing multiple linear regression on the matrix U to obtain a multiple linear equation: y is UB + E, namely a PCR prediction model, and is used for predicting the K value; in the formula, Y is the actually measured K value of the sample set, and E is a random error variable matrix obeying normal distribution;
the parameter B is obtained by the following method: in the modeling process of the PCR prediction model, randomly dividing a plurality of samples into a plurality of correction sets and a plurality of prediction sets; the spectral data of the correction set is used for establishing a prediction model; the spectral data of the prediction set is used for checking the accuracy of the prediction model; and (3) forming the actually measured K values of the correction set samples into a matrix Y to obtain the solution of the equation: b ═ U '(U' U)-1U 'Y, where U' is the transposed matrix of U. And after the parameter matrix B is obtained, the THz spectrum data of the prediction set samples can be substituted into the multivariate linear equation to obtain the prediction K value of the prediction set samples.
In summary, PCR can be divided into two steps: 1, determining the number of principal components, namely, reducing the dimension of a spectrum matrix X by using a principal component analysis method; 2, linear regression analysis is carried out on the dimension-reduced X matrix.
According to the method, principal component decomposition is carried out on the spectrum matrix after pretreatment by using a PCA method, and a principal component set with the cumulative contribution rate reaching 95% is taken as the input of a prediction model. The difference width of the spectrum preprocessing reaches the highest correlation coefficient R according to a prediction modelPAnd lowest RMSEP. Model calculations show that when the first-order difference width is 15, the number of principal component combinations with the cumulative contribution rate of 95% is 26, and the prediction performance of the PCR is optimal (R is the best (R is the first-order difference width)p0.63, RMSEP 16.78%), as shown in fig. 4-1. FIG. 4-2 shows the K value prediction powder of the PCR model under the condition of the parameterDot plots.
2.3.2 error back propagation neural network regression BPANN prediction model
BPANN is a powerful analytical tool that can reveal complex relationships between input and output signals. BPANN is a multi-layer feedforward neural network that can be understood as a non-linear function with the network inputs and predictors being the independent and dependent variables of the function, respectively. BPANN expresses a functional mapping from independent variables to dependent variables. The topology of the BPANN used in the present invention is shown in fig. 5-1.
The topological structure of the BPANN consists of three layers of networks, namely an input layer, a hidden layer and an output layer, wherein each layer comprises a plurality of neurons, and the neurons of the adjacent layers are in one-way connection with weight values;
Iiis the input value of the input layer, i ═ 1,2, …, a; hjIs the output value of the hidden layer, j ═ 1,2, …, b; the output layer being a node
Figure BDA0001250160060000111
The value of (a) is a predicted value of BPANN, namely a predicted K value of fresh meat; hwijIs the connection weight between the neuron of the input layer and the neuron of the hidden layer, owjThe connection weight between the neuron of the hidden layer and the neuron of the output layer;
the output layer being a node, Yl(l ═ 1) is the predicted value of BPANN, i.e. the predicted K value; w is aijAnd wjkThe neuron connection weight of the BPANN.
In order to reduce the complexity of the neural network, the PCA method is adopted to carry out principal component decomposition on the preprocessed spectrum matrix, a principal component set with the cumulative contribution rate of 95% is taken as the input of the neural network, and the node number of a single hidden layer is obtained by an empirical formula:
Figure BDA0001250160060000112
in the formula: number of nodes in hidden layer
Figure BDA0001250160060000113
In the formula: a is the number of nodes of the input layer, c is a constant in the range of 1-10, and 1 is the number of nodes of the output layer. The range of the number of nodes in the hidden layer is selected to be 4-17.
Before the BP neural network prediction, firstly, a network is trained, and the network has associative memory and prediction capabilities through training. The training process of the BPANN comprises the following steps:
1) network initialization: determining a and b according to the input dimension and the output dimension; initializing hwij、owjHidden layer threshold thajThe output layer threshold thb, wherein the weight and the threshold are random numbers in the range of-1, and the given learning rate η is 0.01;
2) hidden layer output calculation: output of hidden layer
Figure BDA0001250160060000121
In the formula, the function f (x) is a hidden layer excitation function, and the expression of the function is:
Figure BDA0001250160060000122
3) output layer output calculation: predicted output of BPANN
Figure BDA0001250160060000123
4) And (3) error calculation: prediction error of BPANN
Figure BDA0001250160060000124
In the formula, y is a K value actually measured by an HPLC method;
5) updating the weight value: hwij=hwij+ηHj(1-Hj)Iiowje,i=1,2,…,a;j=1,2,…,b;
owj=owj+ηHje,j=1,2,…,b;
6) Updating a threshold value: thaj=thaj+ηHj(1-Hj)owje,j=1,2,…,b;thb=thb+e;
7) And judging whether the iteration of the algorithm is finished according to the iteration times and the final iteration error, and returning to the step 2) if the iteration is not finished.
The main characteristics of BPANN are signal forward propagation and error backward propagation. In forward pass, the input signal is processed layer by layer from the input layer through the hidden layer to the output layer. The neuronal state of each layer only affects the neuronal state of the next layer. If the expected output cannot be obtained by the output layer, the backward propagation is carried out, and the network weight and the threshold are adjusted according to the prediction error, so that the predicted output of the BP neural network continuously approaches the expected output. After training is finished, the trained network can be used for predicting the K value.
As can be seen, the BPANN-based K-value prediction modeling comprises three steps of construction, training and prediction of the BPANN. The algorithm flow is shown in fig. 5-2.
Before training the neural network, the input spectrum matrix is normalized to be in a range of-1 to 1, and the actually measured K value is normalized to be in a range of 0 to 1. Because the node number of the hidden layer can influence the performance of neural network prediction, the node number of the hidden layer is selected by adopting a leave-one-out cross-validation method, namely the node number when the RMSECV value is minimum is selected. According to the principle, the numerical change situation of the RMSECV is considered when the prediction model is in the range of the preprocessing difference width 3-47, the node number of the hidden layer is in the range of 4-17, and the RMSECV is shown in the figure 5-3. It can be seen that when the preprocessing difference width is 13 and the number of hidden layer nodes is 9, the RMSECV value is 18.10% at the minimum, and the prediction performance of the BPANN model is the best. FIGS. 5-4 are K-value predicted scatter plots of the BPANN model under this parameter. 2.3.3BPANN improved model
The core idea of the improved algorithm is to train different predictors (weak predictors) for the same correction set, and then assemble the weak predictors to form a stronger final predictor (strong predictor). The algorithm is characterized in that the performance of the weak predictor is poor, but the performance of the strong predictor is excellent. The reason for this is that it is possible to reasonably partition the training set and to reasonably merge weak predictors to form strong predictors. The reasonable points are embodied in the following two aspects: using weighted selected training data instead of randomly selected training samples, thus focusing the training focus on training data samples which are difficult to predict; weak predictors are combined by using a weighted voting mechanism, so that weak predictors with good prediction effects have larger weights, and predictors with poor prediction effects have smaller weights. The improved model takes the BPANN as a weak predictor, repeatedly trains BPANN prediction samples to output, and combines a plurality of obtained BPANN weak predictors into a strong predictor through an improved algorithm.
The specific improved algorithm is as follows:
1. initialization: n is the number of samples in the correction set, and the distribution weight of training data in the initial correction set is Dt(j) Determining a threshold phi of a prediction error of a weak predictor, determining a BPANN weak predictor structure according to the input and output dimensions of a sample, and initializing network parameters;
2. training the weak predictor: training the weak predictor: when the t-th weak predictor is trained, the weak predictor is trained by using training data to obtain a predicted K value output by the weak predictor
Figure BDA0001250160060000131
Then calculating the prediction error
Figure BDA0001250160060000132
T is the number of weak predictors;
3. computing weights w for weak predictorst
Figure BDA0001250160060000133
4. The weight of the test data is adjusted,
Figure BDA0001250160060000134
Figure BDA0001250160060000135
in the formula, BtIs a normalization factor
5. Circulating the steps 2,3 and 4 for T times to obtain T weak predictionsAnd the prediction model is ft. Superposing according to the updated weight to obtain a prediction model of the strong predictor
Figure BDA0001250160060000136
The improved algorithm flow is shown in fig. 6-1, and it can be seen that it is an integrated algorithm, which integrates the prediction capability of each weak predictor in the set by the above-mentioned weighting strategy to obtain the improvement of the overall prediction performance. The invention adopts the BPANN improvement algorithm, uses the improvement algorithm to improve the prediction performance of the BPANN model on the basis of the BPANN model parameters established in the prior art, and constructs the BPANN improvement prediction model of the K value.
The threshold phi of the prediction error has obvious influence on the performance of the BPANN improved model, and the phi with the minimum RMSECV value is selected by adopting a leave-one-out cross verification method: first, phi is selected to be 0.05 to 0.25 in a wide range with a wide step size of 0.01, and the RMSECV value is shown in FIG. 6-2. It can be seen that the performance of the model is good in the range of 0.13-0.14. Then, in this range, 0.13-0.14 is examined with a smaller step size of 0.001, and it can be seen in fig. 6-3 that when Φ is 0.136, the RMSECV value is the smallest, and the BPANN improvement model can obtain the optimal prediction performance. The parameters were iterated through the BPANN improvement algorithm to obtain 10 weak BPANN predictors, whose performance is shown in table 2. It can be seen that: the performance of the weak predictor is different, wherein RMSEC of the weak predictors 2 and 8 is smaller, the prediction performance is better, and the prediction error sum epsilon is largertSmaller, iteratively generated weight wtLarger, meaning that the two weak predictors have a greater weight in the integrated strong predictor; in contrast, the weak predictor 6 has a larger RMSEC, the worst prediction performance, its prediction error sum εtMaximum, iteratively generated weight wtMinimum, meaning that the weak predictor contributes very little in the process of integrating the strong predictor; the prediction performance of the 10 weak predictors is not ideal, but when the weight w is usedtAfter the weak predictors are weighted and integrated, the obtained strong predictors have better performance, and the prediction performance of the strong predictors is better than that of each weak predictor. K-value predicted scatter plots for the BPANN improved model were obtained by this strong predictor, as shown in fig. 6-4.
TABLE 2 predictor Performance of BPANN improvement models
Figure BDA0001250160060000141
2.4 comparison of THz spectral prediction models
The three prediction models adopt K value prediction accuracy of a prediction set data verification model, and the method comprises the following steps: for the prediction set data, taking the prediction set sample spectrum data after first-order differential preprocessing:
Figure BDA0001250160060000142
in the formula, n' is the number of samples in the prediction set, m is the number of spectral wave points, and x is the ATR reflectivity of the samples in the prediction set;
and (3) compressing the spectrum information of the prediction set by using a linear transformation coefficient matrix P in a PCA algorithm of the correction set: is denoted as Un′×a=Xn′×mPm×aTaking each row data I of the matrix Ui(i-1, 2, …, a) as the input of the established prediction model, and the model output, namely the prediction K value, is obtained through the calculation of the prediction model. By the calculation of the formula (5) and the formula (6), the prediction ability of the model can be obtained as an index for evaluating the model performance.
Table 3 systematic comparison of BPANN improvement models with conventional BPANN and Principal Component Regression (PCR) models in terms of performance of THz spectra in predicting K values. It can be seen that the selection of the difference width for spectral data preprocessing 13(BPANN, BPANN-modified model) or 15(PCR model) is suitable, and too small or too large a difference width may degrade the prediction performance of the model. It can also be seen that the prediction performance of the non-linear model (BPANN, BPANN improved model) is significantly better than that of the linear model (PCR model). Furthermore, the BPANN improved model has improved K value prediction performance compared with BPANN. The above results can be explained by the following aspects:
TABLE 3 comparison of regression results for three K-value prediction models
Figure BDA0001250160060000151
1. From the nonlinear trend of the THz spectrum and the pork K value, the K value is determined by the ratio of the contents of 6 correlators of ATP, and the contents of the molecular substances are in a nonlinear relation with the spectral line data in the THz spectrum. The fitting effect of the non-linear model will be better than that of the linear model.
2. From the chemical mechanism of pork deterioration, pork deterioration is a complex chemical process, and under the action of a plurality of putrescence-causing bacteria, proteins in muscles are firstly hydrolyzed into polypeptides, then hydrolyzed into amino acids, and further decomposed into various organic substances. The research in the prior art shows that biological molecules such as protein, polypeptide and amino acid have respective characteristic absorption in the THz wave band, the characteristic absorption mainly comes from the collective vibration mode of the molecules, so that the THz spectrum can reflect the changes. However, the types of biomolecules in pork are various, and the characteristic spectra overlap at normal temperature, so that the THz spectrum does not only express the change information of the K value, but also comprehensively expresses the complex change of the chemical components of the pork. Therefore, the relationship between the K value and the THz spectrum should be a complex nonlinear relationship, and a linear model cannot explain such a complex regression relationship.
3. From the principle structure of a modeling algorithm, the nonlinear model is better in self-learning and self-adjustment than the linear model, and the network topological structure of the BPANN is more suitable for analyzing complex chemical components. The built prediction model is more excellent when complex regression relationships are encountered. The BPANN improvement algorithm enables a BPANN weak prediction model to be gradually improved into a strong prediction model by reasonably integrating the weak predictor in an iteration process. Therefore, the BPANN improvement model provides improved prediction performance over the BPANN model.
2.5 practical significance of the model
R in Table 3PThe RMSEP index represents the accuracy of the model to the prediction of the K value of the unknown sample. The unknown samples, i.e. the samples for which the K-value is unknown, are different from the samples of the calibration and prediction sets for which the K-value is obtained by the method of HPLC for physicochemical analysis. THz light of the processed unknown sample is compressed by a first-order differential preprocessing and PCA spectral information compression methodAnd spectrum data are substituted into the model, so that the K value of the unknown sample can be quickly calculated. According to the scheme, the samples of the correction set are used for modeling, the samples of the prediction set are used for judging the quality of the model, and the K value prediction of the unknown samples is the practicability of the model. The practical meaning of the model is to rapidly detect the K value of an unknown sample through the THz spectrum, rather than modeling and training only with known samples (including a correction set and a prediction set).
3. Conclusion
The research result shows that the THz spectral data of the surface of the lean meat of the cold fresh pig is obtained in the range of 0.1 THz-2 THz by adopting an ATR total reflection mode, and a mathematical model is constructed after first-order differentiation and SG smoothing filtering treatment, so that the K value of the pork can be rapidly detected in a nondestructive mode to evaluate the freshness of the pork. Comparative studies on three prediction models, namely PCR, BPANN and BPANN improved models show that the nonlinear BPANN and BPANN improved models can better fit a complex nonlinear relation between THz spectral data and a freshness K value than a linear PCR model. The BPANN improved model has certain advantages in processing a complex relation model, and the prediction performance of the BPANN is improved to RP0.78, RMSEP 13.45%.
Based on the same THz detection principle, the THz spectrum analysis method can also nondestructively detect the K value of other meats (such as chicken, beef, fish and the like). The invention provides a theoretical basis for further developing and designing portable rapid nondestructive testing equipment based on the method.

Claims (4)

1. A THz spectrum rapid nondestructive testing method for fresh meat freshness K value is characterized by comprising the following steps:
first) meat collection and processing: uniformly cutting fresh meat into meat slices, and refrigerating; when sampling and cutting, the fat and connective tissue of the meat are avoided;
II) measuring the THz spectral data of the meat: adopting THz Attenuated Total Reflection (ATR) mode to quickly detect the sample; when the spectrum is collected, the temperature and the humidity of the environment are kept stable;
thirdly), rapidly detecting the freshness K value of the fresh meat by adopting a relation prediction model between the THz spectrum data and the fresh meat K value: the relation model is a Principal component Regression, a PCR prediction model or a Back Propagation Neural Network (BPANN) prediction model;
substituting the THz spectrum data obtained in the step two) into a relation prediction model, and quickly calculating the K value of the unknown sample;
in the second step):
sampling the THz spectrum of the fresh meat sample in a THz Attenuated Total Reflection (ATR) detection mode; sampling the THz spectrum of the fresh meat sample to obtain sample spectrum data for processing;
all samples used for modeling are randomly divided into a correction set and a prediction set; the spectral data of the correction set is used for establishing a prediction model; the spectral data of the prediction set is used for checking the accuracy of the prediction model;
the sample spectral data matrix is
Figure FDA0002296950350000011
In the formula, n is the number of samples in the correction set, m is the number of spectral wave points, and x is the ATR reflectivity of the samples; xn×mRefers to the spectral data in the correction set;
firstly, compressing the spectrum information: the matrix Xn×mDecomposed into two matrices U and P, denoted as Un×a=Xn×mPm×a(ii) a In the formula, the matrix U represents data xijVector position in the new coordinate system, i ═ 1,2, …, n; j is 1,2, …, m; the column vector of the matrix P represents the linear transformation coefficient of the new coordinate system and the original coordinate system; the dimension a of U is less than the dimension m of X; taking data I of each row of the matrix UiAs an input to the prediction model of step three), i ═ 1,2, …, a;
in the third step):
A. the PCR prediction model comprises the following steps:
performing multiple linear regression on the matrix U to obtain a multiple linear equation: y is UB + E, namely a PCR prediction model, and is used for predicting the K value; in the formula, Y is the actually measured K value of the sample set, and E is a random error variable matrix obeying normal distribution;
the parameter B is obtained by the following method:in the modeling process of the PCR prediction model, randomly dividing a plurality of samples into a plurality of correction sets and a plurality of prediction sets; the spectral data of the correction set is used for establishing a prediction model; the spectral data of the prediction set is used for checking the accuracy of the prediction model; and (3) forming the actually measured K values of the correction set samples into a matrix Y to obtain the solution of the equation: b ═ U '(U' U)-1U 'Y, wherein U' is a transposed matrix of U; after the parameter matrix B is obtained, substituting THz spectrum data of the prediction set samples into the multivariate linear equation to obtain a prediction K value of the prediction set samples;
B. the modeling step of the BPANN prediction model comprises the following steps:
1) constructing a multilayer back propagation neural network:
the method comprises the steps that a multi-layer back propagation neural network BPANN is regarded as a nonlinear function, the input and the predicted value of the BPANN are respectively an independent variable and a dependent variable of the nonlinear function, and the BPANN expresses a function mapping relation from the independent variable to the dependent variable;
the topological structure of the BPANN consists of three layers of networks, namely an input layer, a hidden layer and an output layer, wherein each layer comprises a plurality of neurons, and the neurons of the adjacent layers are in one-way connection with weight values;
Iiis the input value of the input layer, i ═ 1,2, …, a; hjIs the output value of the hidden layer, j ═ 1,2, …, b; the output layer being a node
Figure FDA0002296950350000021
Figure FDA0002296950350000022
The value of (a) is a predicted value of BPANN, namely a predicted K value of fresh meat; hwijIs the connection weight between the neuron of the input layer and the neuron of the hidden layer, owjThe connection weight between the neuron of the hidden layer and the neuron of the output layer;
number of nodes in hidden layer
Figure FDA0002296950350000023
In the formula: a is the number of nodes of the input layer, c is a constant in the range of 1-10, and 1 is the output layerThe number of nodes;
2) training the model of step 1), comprising:
2.1) network initialization: determining a and b according to the input dimension and the output dimension; initializing hwij、owjHidden layer threshold thajThe output layer threshold thb, wherein the weight and the threshold are random numbers in the range of-1, and the given learning rate η is 0.01;
2.2) hidden layer output calculation: output of hidden layer
Figure FDA0002296950350000024
In the formula, the function f (x) is a hidden layer excitation function, and the expression of the function is:
Figure FDA0002296950350000025
2.3) output layer output calculation: predicted output of BPANN
Figure FDA0002296950350000026
2.4) error calculation: prediction error of BPANN
Figure FDA0002296950350000027
In the formula, y is a K value actually measured by an HPLC method;
2.5) weight updating: hwij=hwij+ηHj(1-Hj)Iiowje,i=1,2,…,a;j=1,2,…,b;owj=owj+ηHje,j=1,2,…,b;
2.6) threshold update: thaj=thaj+ηHj(1-Hj)owje,j=1,2,…,b;thb=thb+e;
2.7) judging whether the iteration of the algorithm is finished according to the iteration times and the final iteration error, and returning to the step 2.2) if the iteration is not finished.
2. The rapid THz spectrum nondestructive testing method for fresh meat freshness K value according to claim 1, wherein the modeling step of the BPANN model further comprises the step 3) of optimizing the BPANN model to obtain a BPANN improved model by optimizing the BPANN prediction model as a weak predictor to obtain a strong predictor, and the step comprises the steps of:
3.1) initialization: n is the number of samples in the correction set, and the distribution weight of training data in the initial correction set is Dt(j) Determining a threshold phi of a prediction error of a weak predictor, determining a BPANN weak predictor structure according to the input and output dimensions of the sample, and initializing network parameters, wherein j is 1/n, j is 1,2, …, and n;
3.2) training weak predictors: when the t-th weak predictor is trained, the weak predictor is trained by using training data to obtain a predicted K value output by the weak predictor
Figure FDA0002296950350000031
Then calculating the prediction error
Figure FDA0002296950350000032
T is the number of weak predictors, T is 1,2, …, T;
3.3) computing the weights of the Weak predictors
Figure FDA0002296950350000033
3.4) adjusting the distribution weight of the training data and normalizing:
Figure FDA0002296950350000034
in the formula (I), the compound is shown in the specification,
Figure FDA0002296950350000035
3.5) circulating the steps of 3.2-3.4) for T times to obtain T weak predictors, wherein the prediction model is ftSuperposing according to the updated weight to obtain a prediction model of the strong predictor
Figure FDA0002296950350000036
The phi value is the phi value when the performance of the model is best selected by adopting a leave-one cross verification method, namely, when the phi value is different in a certain range and is investigated by a fixed step length, the RMSECV value under the leave-one cross verification method in the modeling step is taken as the optimal threshold value when the RMSECV value is minimum;
the operation process of the leave-one-cross verification method is as follows: in n samples of the correction set, each sample is taken as a prediction set independently, and the rest n-1 samples are taken as the correction set; obtaining a predicted value of a sample of the prediction set through a modeling step by taking the correction set as a basis; obtaining the predicted values of all samples through n times of circulation;
obtaining the performance index RMSECV of a leave-one-cross verification method,
Figure FDA0002296950350000037
3. the rapid THz spectrum nondestructive testing method for fresh meat freshness K value according to claim 1, wherein in the second step), the spectrum frequency band is 0.1 THz-2 THz.
4. The rapid THz spectrum nondestructive testing method for fresh meat freshness K value according to claim 1, wherein in step two), sample spectrum data is preprocessed, and the preprocessed data is used as data of a matrix X; the method for preprocessing the sample spectrum data comprises the following steps:
first, preprocessing the spectral data by using a first derivative FD:
Figure FDA0002296950350000038
wherein x is the ATR reflectivity of the sample at the number of wave points i; when i is 1, its first order differential value is
Figure FDA0002296950350000039
… …, respectively; when i is m, its first order differential value is
Figure FDA00022969503500000310
And smoothing the spectrum by adopting an SG polynomial: selecting adjacent data points in the odd number of spectral data sequences, forming a window by the data points, and smoothing the central point of the window; in the smoothing process, fitting a polynomial to the data points in the window, and calculating the smoothed points according to the obtained polynomial; the expression of the data point after smoothing is obtained as follows:
Figure FDA0002296950350000041
wherein the width of the window is 2d +1, s is constant, and λjSmoothing the coefficient sequence by SG;
the fitting of the polynomial uses a least squares method.
CN201710167209.0A 2017-03-20 2017-03-20 THz spectrum rapid nondestructive testing method for fresh meat freshness K value Expired - Fee Related CN106960091B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710167209.0A CN106960091B (en) 2017-03-20 2017-03-20 THz spectrum rapid nondestructive testing method for fresh meat freshness K value

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710167209.0A CN106960091B (en) 2017-03-20 2017-03-20 THz spectrum rapid nondestructive testing method for fresh meat freshness K value

Publications (2)

Publication Number Publication Date
CN106960091A CN106960091A (en) 2017-07-18
CN106960091B true CN106960091B (en) 2020-05-05

Family

ID=59470384

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710167209.0A Expired - Fee Related CN106960091B (en) 2017-03-20 2017-03-20 THz spectrum rapid nondestructive testing method for fresh meat freshness K value

Country Status (1)

Country Link
CN (1) CN106960091B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109030405A (en) * 2018-05-25 2018-12-18 深圳市太赫兹科技创新研究院有限公司 The method of Electromagnetic Wave Detection cordyceps sinensis quality based on Terahertz frequency range
CN113218878A (en) * 2020-01-21 2021-08-06 青岛海尔电冰箱有限公司 Method for judging freshness of food materials of refrigerator, refrigerator and storage medium
CN113008172B (en) * 2021-03-03 2022-09-27 北京理工大学 Terahertz wave-based ice and snow track inspection device and method

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102507459A (en) * 2011-11-23 2012-06-20 中国农业大学 Method and system for quick lossless evaluation on freshness of fresh beef
CN103487382A (en) * 2013-09-27 2014-01-01 浙江工商大学 Method for judging tuna meat freshness by redness index
CN104655761A (en) * 2015-02-28 2015-05-27 华南理工大学 Method for measuring fish freshness index value K in online way based on multi-spectral imaging
CN105548029A (en) * 2015-12-14 2016-05-04 北京农业质量标准与检测技术研究中心 Meat product freshness detection method based on spectral imaging technology

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102507459A (en) * 2011-11-23 2012-06-20 中国农业大学 Method and system for quick lossless evaluation on freshness of fresh beef
CN103487382A (en) * 2013-09-27 2014-01-01 浙江工商大学 Method for judging tuna meat freshness by redness index
CN104655761A (en) * 2015-02-28 2015-05-27 华南理工大学 Method for measuring fish freshness index value K in online way based on multi-spectral imaging
CN105548029A (en) * 2015-12-14 2016-05-04 北京农业质量标准与检测技术研究中心 Meat product freshness detection method based on spectral imaging technology

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Non-destructively sensing pork’s freshness indicator using near infrared multispectral imaging technique;Qiping Huang 等;《Journal of Food Engineering》;20150117;第154卷;第69-75页 *
The meat freshness detection based on terahertz wave;Xin B. Huang 等;《Proceedings SPIE 9795,Selected Papers of the Photoelectronic Technology Committee Conferences held June–July 2015》;20151105;第97953D页 *
冷鲜猪肉的新鲜度无损检测技术现状及THz检测技术展望;齐亮 等;《食品与机械》;20160930;第32卷(第9期);第219-224页 *
基于光谱技术和支持向量机的生鲜猪肉水分含量快速无损检测;张海云 等;《光谱学与光谱分析》;20121015;第32卷(第10期);第2794-2798页 *

Also Published As

Publication number Publication date
CN106960091A (en) 2017-07-18

Similar Documents

Publication Publication Date Title
Zhang et al. Nondestructive measurement of soluble solids content in apple using near infrared hyperspectral imaging coupled with wavelength selection algorithm
Liu et al. Recent advances in wavelength selection techniques for hyperspectral image processing in the food industry
Kartakoullis et al. Feasibility study of smartphone-based Near Infrared Spectroscopy (NIRS) for salted minced meat composition diagnostics at different temperatures
Meza-Márquez et al. Application of mid-infrared spectroscopy with multivariate analysis and soft independent modeling of class analogies (SIMCA) for the detection of adulterants in minced beef
CN107014769B (en) A kind of fresh meat K value Fast nondestructive evaluation model based on THz spectrum analysis
Chen et al. A new approach to near-infrared spectral data analysis using independent component analysis
Zhang et al. Rapid evaluation of texture parameters of Tan mutton using hyperspectral imaging with optimization algorithms
Qu et al. Predicting pork freshness using multi-index statistical information fusion method based on near infrared spectroscopy
Bonah et al. Comparison of variable selection algorithms on vis-NIR hyperspectral imaging spectra for quantitative monitoring and visualization of bacterial foodborne pathogens in fresh pork muscles
Liu et al. A comparative study for least angle regression on NIR spectra analysis to determine internal qualities of navel oranges
Baek et al. Shortwave infrared hyperspectral imaging system coupled with multivariable method for TVB-N measurement in pork
CN106960091B (en) THz spectrum rapid nondestructive testing method for fresh meat freshness K value
Li et al. Combining Vis-NIR and NIR hyperspectral imaging techniques with a data fusion strategy for the rapid qualitative evaluation of multiple qualities in chicken
Che et al. Application of visible/near‐infrared spectroscopy in the prediction of azodicarbonamide in wheat flour
Fan et al. Rapid determination of TBARS content by hyperspectral imaging for evaluating lipid oxidation in mutton
Cao Calibration optimization and efficiency in near infrared spectroscopy
Fatemi et al. Identification of informative spectral ranges for predicting major chemical constituents in corn using NIR spectroscopy
Zhang et al. Rapid identification of lamb freshness grades using visible and near-infrared spectroscopy (Vis-NIR)
von Gersdorff et al. Method comparison between real-time spectral and laboratory based measurements of moisture content and CIELAB color pattern during dehydration of beef slices
Mishra et al. Improved prediction of minced pork meat chemical properties with near-infrared spectroscopy by a fusion of scatter-correction techniques
He et al. Hyperspectral imaging combined with chemometrics for rapid detection of talcum powder adulterated in wheat flour
Andersen et al. Prediction of water holding capacity and pH in porcine longissimus lumborum using Raman spectroscopy
Li et al. Raman spectroscopy combined with support vector regression and variable selection method for accurately predicting salmon fillets storage time
Jiang et al. Comparison of wavelength selected methods for improving of prediction performance of PLS model to determine aflatoxin B1 (AFB1) in wheat samples during storage
Shao et al. Multivariate calibration of near-infrared spectra by using influential variables

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200505

CF01 Termination of patent right due to non-payment of annual fee