Cheese processing monitoring control method based on double spectrums and spectrum analysis condensation process

Abstract

The invention relates to the field of cheese processing control systems, in particular to a cheese processing monitoring control method based on a double-spectrum and spectrum analysis condensation process, which comprises the following steps: acquiring a terahertz spectrum signal sequence and a laser speckle image feature sequence based on a double-spectrum system, and constructing a molecular scale feature matrix and a micro scale dynamic feature matrix; and carrying out time point-by-time feature association on the micro-scale dynamic feature matrix and the molecular scale feature matrix, and establishing a cross-scale condensation feature matrix. According to the method, a topological data analysis mechanism is introduced to carry out topological modeling on the preliminary fusion feature matrix, so that the macroscopic stage conversion rule of the condensation process is clearly presented; and the constructed macroscopic trend characteristic matrix is complementary with the molecular and microscopic characteristic matrixes to finally form a cross-scale coagulation characteristic matrix which comprehensively covers key information of three dimensions of molecular structure, microscopic morphology and macroscopic trend, so that multi-dimensional characterization of the coagulation process is realized, and a comprehensive characteristic basis is provided for accurate judgment of the coagulation state.

Description

Technical Field

[0001] This invention relates to the field of cheese processing control systems, and more specifically to a cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process. Background Technology

[0002] Milk coagulation is the core process in cheese processing. It involves multi-scale dynamic reactions such as protein molecule polymerization, fat globule interface reconstruction, and microscopic viscosity changes, which directly determine the texture, flavor stability, and product qualification rate of cheese.

[0003] However, existing methods for monitoring coagulation in cheese processing have many technical limitations. For example, the accuracy of molecular feature analysis is insufficient; traditional spectral fitting methods struggle to accurately separate the characteristic parameters of core components such as protein amide bonds and phospholipid molecules at the fat globule interface; and redundant correlations between features are not effectively addressed, resulting in high redundancy in the feature matrix and interference with core information. Furthermore, current technologies focus on protein polymerization during coagulation, while microscopic signal processing exhibits significant nonlinearity and dynamic fluctuations. Traditional linear decomposition methods cannot effectively remove the essential feature components, accurately identify key nodes in the coagulation stage transition, or achieve refined closed-loop control of the coagulation process, leading to significant fluctuations in cheese product quality. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a method for monitoring and controlling cheese processing based on dual-spectrum and spectral analysis of the coagulation process.

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: a cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process, the method comprising: S101. Based on the terahertz dynamic evolution signal and speckle texture dynamic evolution signal obtained by the dual-spectrum system, time axis alignment and preprocessing are performed to obtain the terahertz spectral signal sequence and the laser speckle image feature sequence. S102. The spectral signal intensity of the terahertz spectral signal sequence is enhanced by coherent anti-Stokes Raman scattering spectroscopy. The enhanced terahertz spectral signal sequence is then fitted with Lorentz to construct a molecular-scale feature matrix. A nonlinear signal analysis strategy is embedded to decompose the laser speckle image feature sequence into effective intrinsic mode function components. After transforming the effective intrinsic mode function components, a micro-scale dynamic feature matrix is ​​constructed. The micro-scale dynamic feature matrix and the molecular-scale feature matrix are correlated at each time point to establish a preliminary fused feature matrix. Then, a topological data analysis mechanism is introduced to perform topological modeling. Topological invariants are extracted under a sliding time window. The extracted topological invariants constitute a macro-trend feature matrix. The micro-scale dynamic feature matrix, the macro-trend feature matrix, and the molecular-scale feature matrix are fused to obtain the final cross-scale condensation feature matrix. S103. Extract the covariance matrix from the cross-scale condensation feature matrix, perform eigenvalue decomposition on the extracted covariance matrix, extract the largest eigenvalue and then normalize it to obtain the comprehensive response value and the rate of change of the comprehensive response value, and determine the adjustment direction based on the rate of change of the comprehensive response value.

[0006] In a preferred embodiment, step S101 integrates a dual-spectrum system, including a terahertz time-domain spectral emission module, a terahertz receiving module, a laser speckle interferometric emission module, and an image acquisition module, onto an adjustable support above the cheese processing container. The detection directions of the terahertz time-domain spectral emission module and the laser speckle interferometric emission module are perpendicular to the milk surface, and the detection areas of the dual-spectrum system completely overlap. The terahertz time-domain spectral emission module and the terahertz receiving module are used to capture the terahertz wave amplitude attenuation, phase shift, and refractive index change signals caused by protein molecule polymerization and fat globule interface changes during milk coagulation, accurately capturing the molecular structure changes during milk coagulation. The laser speckle interferometric emission module and the image acquisition module are used to capture the dynamic evolution signals of speckle texture caused by micro-viscosity changes during milk coagulation, such as speckle particle size, grayscale distribution, and motion trajectory changes. By controlling the start-up and data acquisition timing of the dual-spectrum system, time axis alignment is completed. After baseline correction and noise suppression, a time axis-aligned terahertz spectral signal sequence and a laser speckle image feature sequence are obtained.

[0007] In a preferred embodiment, after acquiring the terahertz spectral signal sequence, step S102, based on coherent anti-Stokes Raman scattering spectral analysis, probes the milk surface with pump light and Stokes light. The pump light and Stokes light excite a coherent anti-Stokes Raman scattering effect. This effect amplifies the Raman scattering signal of the target component through a nonlinear optical process, thereby enhancing the weak signal. The coherent anti-Stokes Raman scattering effect is a third-order nonlinear Raman scattering process. Coherent anti-Stokes Raman scattering light is generated by the interaction of two laser beams with different frequencies with the sample, obtaining the enhanced terahertz spectral signal sequence, and then... The intensity value of the spectral signal after Lorentz fitting is obtained by dividing the amplitude of the characteristic peak by pi, multiplying by the square of half-width at half-maximum (FWHM) of the characteristic peak, and subtracting the sum of the squares of the center frequencies of the characteristic peaks. Based on the fitted curve, protein amide bond characteristic peak parameters including characteristic peak amplitude, FWHM, and peak position shift are extracted, as well as fat globule interface phospholipid molecule characteristic peak parameters including the rate of change of characteristic peak amplitude. Furthermore, the characteristic peak amplitude reflects the content and degree of polymerization of protein molecules, the FWHM reflects the uniformity of molecular vibration, the peak position shift reflects the change of the chemical environment of the molecule, and the rate of change of characteristic peak amplitude reflects the dynamic evolution of the fat globule interface structure. Using the time point aligned with the S101 time axis as the horizontal axis, and the first feature dimension, composed of peak amplitude, full width at half maximum (FWHM), peak position offset, and the rate of change of molecular peak amplitude at the fat globule interface, as the vertical axis, a preliminary molecular-scale feature matrix is ​​constructed. Each element in this matrix corresponds to the quantized value of a specific molecular feature dimension at a specific time point.

[0008] In a preferred embodiment, after obtaining the preliminary molecular-scale feature matrix, Pearson correlation analysis is used to pair each first feature dimension at the same time point to form a first correlation array. The covariance of the time-series signal data corresponding to the two first feature dimensions is divided by the product of the standard deviations of the two first feature dimensions, which is used as the first correlation coefficient of each first correlation array. Any first feature dimension at the same time point is removed from the first correlation array whose first correlation coefficient exceeds twice the mean of the first correlation coefficient. When the first correlation coefficient exceeds twice the mean of the correlation, it indicates that the two first feature dimensions are... The linear correlation is extremely strong, and retaining any one feature can cover the core information, which is convenient for subsequent processing. The results of pairing each first feature dimension at the same time point specifically include: protein amide bond peak amplitude and half-width at half-maximum (WHM), protein amide bond peak amplitude and peak position shift, protein amide bond peak amplitude and fat globule interface molecular peak amplitude change rate, protein amide bond WHM and peak position shift, protein amide bond WHM and fat globule interface molecular peak amplitude change rate, and protein amide bond peak position shift and fat globule interface molecular peak amplitude change rate. Finally, the molecular scale feature matrix is ​​constructed through the above steps.

[0009] In a preferred embodiment, after analysis by coherent anti-Stokes Raman scattering spectroscopy, features at the molecular scale are extracted and merged. However, features at the microscopic level require another approach. Therefore, after constructing the molecular-scale feature matrix, in step S102, the laser speckle image feature sequence, including the equivalent diameter of speckle particles, gray mean, gray variance, and texture entropy as the input of the nonlinear signal analysis strategy, is used as the input of the nonlinear signal analysis strategy. After feature extraction, the image feature sequence to be decomposed, including the equivalent diameter of speckle particles, gray mean, gray variance, and texture entropy, is obtained respectively. In the feature sequence of the image to be decomposed, if the feature value of the second feature dimension at the same time point is greater than the feature values ​​of the second feature dimension at each of the two adjacent time points, it is used as the determination of a local maximum point. If the feature value of the second feature dimension at the same time point is less than the feature values ​​of the second feature dimension at each of the two adjacent time points, it is used as the determination of a local minimum point. By traversing the time-series signal data of each feature sequence of the image to be decomposed, all local maximum points and local minimum points are identified, reflecting the fluctuation characteristics of the feature sequence signal of the image to be decomposed. The extracted sets of local maxima and local minima are fitted separately to obtain two continuous envelopes. The upper envelope covers all local maxima and the lower envelope covers all local minima. During the fitting process, it is ensured that the envelopes completely cover all current data points. The arithmetic mean of the upper and lower envelope values ​​at the same time point is calculated to obtain the mean of the envelope at the same time point. The mean of the envelope at all time points is then used to form a continuous mean curve. Subtract the mean value of the mean curve at each time point from the feature sequence of each image to be decomposed to obtain the corresponding candidate image feature signal sequence. Verify the candidate image feature signal sequence based on the intrinsic mode function to obtain the effective intrinsic mode function components. The verification logic includes that the absolute value of the mean of the upper and lower envelopes at all time points is less than or equal to 0.05 times the maximum signal amplitude, that is, the upper and lower envelopes of the candidate sequence are approximately symmetrical about the time axis, and the number of zero crossings in the candidate sequence is equal to or differs from the number of extreme points by no more than 1. The criterion for zero crossing is that the signal values ​​at adjacent time points have opposite signs, thus completing the decomposition of the laser speckle image feature sequence.

[0010] In a preferred embodiment, after obtaining the effective intrinsic mode function components, the transformed signal sequence, i.e., the orthogonal component, is obtained by integral operation based on the characteristic value of each time point of the effective intrinsic mode function components. The orthogonal component is perfectly aligned with the original effective intrinsic mode function components in time, and the phase difference between the two is always 90 degrees. The effective intrinsic mode function components are taken as the real part and the transformed signal sequence is taken as the imaginary part, and the combination is used to form an analytical signal corresponding to a set of complex data of real part values ​​and imaginary part values ​​at each time point. The real and imaginary parts of the analytic signal at each time point are squared, the two squared results are added together, and finally the square root of the sum is taken to obtain the instantaneous amplitude. The instantaneous amplitude is the amplitude intensity of the analytic signal at the current time point, which directly reflects the dynamic change intensity of the effective intrinsic mode function components. For the analytical signal at each time point, the phase angle is obtained by arctangent operation. That is, with the imaginary part as the numerator and the real part as the denominator, the phase angle difference between the current time point and the previous adjacent time point is calculated by using the adjacent time point difference method. Then, the difference is divided by the time interval between the two time points to obtain the instantaneous frequency of the current time point.

[0011] Furthermore, the instantaneous frequency value at the current time point is subtracted from the instantaneous frequency value at the previous adjacent time point to obtain the frequency difference. The frequency difference is then divided by the time interval to obtain the instantaneous frequency change rate at the current time point. A number of consecutive time points that can cover a sufficient time series range are taken as a sliding window. The instantaneous amplitude values ​​of the time points within the window are added together and then divided by the number of time points to obtain the mean value within the window. The difference between the instantaneous amplitude value and the mean value at each time point within the window is calculated sequentially. Each difference is squared and then summed. The sum is then divided by the number of time points minus 1. Finally, the square root of the result is taken to obtain the instantaneous amplitude fluctuation range at the current time point. For all time points of the effective intrinsic mode function components corresponding to the feature sequence of the image to be decomposed, calculate the phase angle of any two effective components at the same time point to obtain the instantaneous phase difference at that time point. Calculate the phase difference by pairing them up as described above to obtain an instantaneous phase difference time series curve and phase angle difference value corresponding to each pairing combination.

[0012] In a preferred embodiment, S102 calculates the instantaneous frequency change rate and instantaneous amplitude fluctuation amplitude using the sampling frequency of the dual-spectral system as the time interval, and obtains the instantaneous phase difference. With the time point as the horizontal axis and the instantaneous frequency change rate, instantaneous phase difference, and instantaneous amplitude fluctuation amplitude as the vertical axis, a microscale dynamic feature matrix with a third feature dimension is constructed.

[0013] The third feature dimension of this microscale dynamic feature matrix includes the instantaneous frequency change rate of the effective intrinsic mode function corresponding to the equivalent diameter of speckle particles, the instantaneous phase difference of the effective intrinsic mode function corresponding to the equivalent diameter of speckle particles, the instantaneous amplitude fluctuation of the effective intrinsic mode function corresponding to the gray mean, the instantaneous frequency change rate of the effective intrinsic mode function corresponding to the gray variance, the instantaneous phase difference of the effective intrinsic mode function corresponding to the gray variance, and the instantaneous amplitude fluctuation of the effective intrinsic mode function corresponding to the texture entropy.

[0014] After obtaining the microscale dynamic feature matrix and the molecular scale feature matrix, Pearson correlation analysis is performed on a time-by-time basis to pair all the first feature dimensions and the third feature dimension at the same time point. The first feature dimension represents the molecular scale feature, and the third feature dimension represents the microscale feature. Redundant features are removed by calculating the second correlation coefficient. Referring to the scheme disclosed in the first correlation coefficient in the above steps, after completing the time-by-time feature association, a preliminary fused feature matrix is ​​established.

[0015] In a preferred embodiment, step S102 calculates the pairwise Euclidean distances between feature vectors at all time points in the preliminary fused feature matrix, and constructs a minimum spanning tree based on the Euclidean distances. The minimum spanning tree is an algorithm for constructing a connected graph, ensuring that all vertices, i.e., the feature vectors at all time points, are connected and that the total edge length, i.e., the feature vector distance, is minimized. This is used to determine a reasonable distance threshold for topology modeling, and the distance threshold is set to 1.2 times the maximum edge length in the minimum spanning tree. This multiple ensures that the core connectivity of the feature vectors is preserved intact in the topology structure, while avoiding fragmentation of the topology structure due to an excessively small threshold. A sliding time window is used to divide the full-time-series feature vectors in the preliminary fused feature matrix to ensure sufficient coverage of the time series to characterize macroscopic changes. For each feature vector within the sliding time window, all other vertices with distances less than the distance threshold are connected and edges are generated to construct a topological complex. Topological invariants, including the Betti number, Euler feature number, and persistent homology barcode parameters, are extracted from the topological complex of each sliding time window.

[0016] The topological complex is used to establish the relationship between topological structure and time series. The Betti number, Euler feature number and persistent homology barcode parameters extracted from the topological complex represent the core invariant of the connectivity of the topological space, the comprehensive invariant of the overall structure of the topological space, and the invariant of the appearance and disappearance process of topological features, respectively. The key extraction methods are the 0th and 1st order Betti numbers. The 0th order Betti number corresponds to the number of connected components in the complex, and the 1st order Betti number corresponds to the number of pores in the complex. These are extracted by calculating the homology group of the complex. The Euler eigenvalue is calculated based on the simplex dimension, which is the number of vertices minus the number of edges plus the number of faces minus the number of volumes. It reflects the stage transition of the coagulation process, corresponding to the completion of the transformation of milk from liquid to gel state. The 0th and 1st order Betti numbers represent the degree of aggregation of microscopic gel particles during the coagulation process. The decrease in the number of connected components indicates an enhanced particle aggregation trend, reflecting the initial formation of the gel network. The number of pores first increases and then decreases, corresponding to the process of gel particles aggregating to form a network. The persistent coherence barcode first gradually increases the distance threshold from 0 to twice the set threshold, records the threshold when the feature appears and the threshold when the feature disappears for each Betti number, forming a barcode with the threshold as the horizontal axis and the duration as the vertical axis. The barcode length reflects the stability of the topological feature; the longer the length, the smoother the change in the coagulation process. The amplitude corresponds to the size of the Betti number, reflecting the strength of the topological feature.

[0017] In a preferred embodiment, after obtaining the topological invariants including the Betti number, Euler eigenvalue, and persistent cohomology barcode parameters, a macroscopic trend feature matrix is ​​constructed with the midpoint of the time corresponding to the sliding window as the horizontal axis and the topological invariants including the Betti number, Euler eigenvalue, and persistent cohomology barcode parameters as the vertical axis. An adaptive weight fusion strategy is used to fuse the microscale dynamic feature matrix, the macroscopic trend feature matrix, and the molecular scale feature matrix. Through the particle swarm optimization algorithm, with the objective function of maximizing the goodness of fit between the fused features and the degree of condensation, the optimized weight combination is finally output, and a cross-scale condensation feature matrix is ​​formed by a fusion logic of weighted summation at each time point.

[0018] In a preferred embodiment, S103 calculates the covariance matrix of the obtained cross-scale condensation feature matrix, performs eigenvalue decomposition on the covariance matrix, extracts the maximum eigenvalue, the maximum eigenvalue reflects the core information contribution of the covariance matrix, and can characterize the overall dominant change trend of cross-scale fusion features. After normalization, a comprehensive response value is obtained. The comprehensive response value obtained by the dual-spectral system in the previous and next cycles is used to calculate the rate of change of the comprehensive response value. If the rate of change of the comprehensive response value is greater than zero, the adjustment is reversed; if the comprehensive response value is less than zero, the adjustment is forward; if the comprehensive response value is equal to zero, the operating data of the condensation reaction component is kept unchanged, and the comprehensive response value is used as the adjustment amount of the operating data of the reaction component.

[0019] The beneficial effects of this invention are as follows: Through the feature association mechanism at each time point, the organic connection between molecular scale and microscale features is realized, avoiding the limitations of single-scale features. The introduction of a topological data analysis mechanism to perform topological modeling on the preliminary fused feature matrix can clearly present the macroscopic stage transformation law of the condensation process. The constructed macroscopic trend feature matrix complements the molecular and microscopic feature matrices. The resulting cross-scale condensation feature matrix comprehensively covers the key information of the three dimensions of molecular structure, micromorphology, and macroscopic trend, realizing the multi-dimensional characterization of the condensation process and providing comprehensive feature basis for accurate judgment of the condensation state. The intensity of the terahertz spectral signal was enhanced by coherent anti-Stokes Raman scattering spectroscopy. Combined with Lorentz fitting, key characteristic parameters of protein amide bonds and phospholipid molecules at the fat globule interface could be accurately separated. Then, highly redundant feature dimensions were eliminated by Pearson correlation analysis. The constructed molecular-scale feature matrix not only retains the core information of the dynamic evolution of molecular structure, but also avoids the analytical bias caused by redundant features, providing highly reliable feature support for the analysis of the molecular mechanism of the condensation process. Attached Figure Description

[0020] Figure 1 This is a flowchart of the present invention. Detailed Implementation

[0021] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0022] As attached Figure 1 As shown, this embodiment provides a method for monitoring and controlling cheese processing based on dual-spectrum and spectral analysis of the coagulation process. The method includes: S101. Based on the terahertz dynamic evolution signal and speckle texture dynamic evolution signal obtained by the dual-spectrum system, time axis alignment and preprocessing are performed to obtain the terahertz spectral signal sequence and the laser speckle image feature sequence. A dual-spectrum system, comprising a terahertz time-domain spectroscopy emission module, a terahertz receiving module, a laser speckle interferometry emission module, and an image acquisition module, is integrated into an adjustable support above a cheese processing container. The detection directions of the terahertz time-domain spectroscopy emission module and the laser speckle interferometry emission module are perpendicular to the milk surface, and the detection areas of the dual-spectrum system completely overlap. The terahertz time-domain spectroscopy emission module and the terahertz receiving module are used to capture the terahertz wave amplitude attenuation, phase shift, and refractive index change signals caused by protein molecule polymerization and fat globule interface changes during milk coagulation, accurately capturing the molecular structure changes during milk coagulation. The laser speckle interferometry emission module and the image acquisition module are used to capture the dynamic evolution signals of speckle texture caused by micro-viscosity changes during milk coagulation, such as speckle particle size, grayscale distribution, and motion trajectory changes. By controlling the start-up and data acquisition timing of the dual-spectrum system, time axis alignment is completed. After baseline correction and noise suppression, a time-axis aligned terahertz spectral signal sequence and a laser speckle image feature sequence are obtained.

[0023] S102. The spectral signal intensity of the terahertz spectral signal sequence is enhanced by coherent anti-Stokes Raman scattering spectroscopy. The enhanced terahertz spectral signal sequence is then fitted with Lorentz to construct a molecular-scale feature matrix. A nonlinear signal analysis strategy is embedded to decompose the laser speckle image feature sequence into effective intrinsic mode function components. After transforming the effective intrinsic mode function components, a micro-scale dynamic feature matrix is ​​constructed. The micro-scale dynamic feature matrix and the molecular-scale feature matrix are correlated at each time point to establish a preliminary fused feature matrix. Then, a topological data analysis mechanism is introduced to perform topological modeling. Topological invariants are extracted under a sliding time window. The extracted topological invariants constitute a macro-trend feature matrix. The micro-scale dynamic feature matrix, the macro-trend feature matrix, and the molecular-scale feature matrix are fused to obtain the final cross-scale condensation feature matrix. A1. After acquiring the terahertz spectral signal sequence, based on coherent anti-Stokes Raman scattering spectral analysis, pump light and Stokes light are probed onto the surface of milk. The pump light and Stokes light excite the coherent anti-Stokes Raman scattering effect. This effect amplifies the Raman scattering signal of the target component through a nonlinear optical process, thereby enhancing the weak signal. The coherent anti-Stokes Raman scattering effect is a third-order nonlinear Raman scattering process. By having two laser beams of different frequencies interact with the sample to generate coherent anti-Stokes Raman scattering light, the enhanced terahertz spectral signal sequence is obtained, and the characteristic peak amplitude of the terahertz spectral signal sequence is determined. The sum of the squares of the half-width at half-maximum (WHM) of the characteristic peaks, divided by pi and subtracted from the squares of the center frequencies of the characteristic peaks, is used to obtain the spectral signal intensity value after Lorentz fitting. Based on the fitted curve, protein amide bond characteristic peak parameters, including characteristic peak amplitude, WHM, and peak position shift, are extracted, as well as phospholipid molecule characteristic peak parameters at the fat globule interface, including the rate of change of characteristic peak amplitude. Furthermore, the characteristic peak amplitude reflects the content and degree of polymerization of protein molecules, the WHM reflects the uniformity of molecular vibration, the peak position shift reflects the change of the chemical environment in which the molecules are located, and the rate of change of characteristic peak amplitude reflects the dynamic evolution of the fat globule interface structure. A2. Using the time points aligned with the S101 time axis as the horizontal axis and the first feature dimension, composed of peak amplitude, full width at half maximum (FWHM), peak position offset, and the rate of change of molecular peak amplitude at the fat globule interface, as the vertical axis, a preliminary molecular-scale feature matrix is ​​constructed. Each element in this matrix corresponds to the quantized value of a specific molecular feature dimension at a specific time point.

[0024] A3. After obtaining the preliminary molecular-scale feature matrix, based on Pearson correlation analysis, each first feature dimension at the same time point is paired to form a first correlation array. The covariance of the time series signal data corresponding to the two first feature dimensions is divided by the product of the standard deviations of the two first feature dimensions, which is used as the first correlation coefficient for each first correlation array. Any first feature dimension at the same time point is removed from the first correlation array whose first correlation coefficient exceeds twice the mean of the first correlation coefficient. When the first correlation coefficient exceeds twice the mean of the correlation, it indicates that the two first feature dimensions are linearly correlated. The feature is highly accurate, and retaining any one feature can cover the core information, which is convenient for subsequent processing. The results of pairing each first feature dimension at the same time point specifically include: protein amide bond peak amplitude and half-width at half-maximum (WHM), protein amide bond peak amplitude and peak position shift, protein amide bond peak amplitude and fat globule interface molecular peak amplitude change rate, protein amide bond WHM and peak position shift, protein amide bond WHM and fat globule interface molecular peak amplitude change rate, and protein amide bond peak position shift and fat globule interface molecular peak amplitude change rate. Finally, the molecular scale feature matrix is ​​constructed through the above steps.

[0025] B1. After analyzing the coherent anti-Stokes Raman scattering spectrum, the features at the molecular scale are extracted and merged. However, the features at the microscopic level need to be processed through another method. Therefore, after constructing the molecular scale feature matrix, S102 uses the laser speckle image feature sequence, which includes the equivalent diameter of speckle particles, gray mean, gray variance, and texture entropy as the input of the nonlinear signal analysis strategy. After feature extraction, the image feature sequence to be decomposed, which includes the equivalent diameter of speckle particles, gray mean, gray variance, and texture entropy, is obtained respectively. B2. In the feature sequence of the image to be decomposed, if the feature value of the second feature dimension at the same time point is greater than the feature value of the second feature dimension at each of the two adjacent time points, it is used as the determination of a local maximum point. If the feature value of the second feature dimension at the same time point is less than the feature value of the second feature dimension at each of the two adjacent time points, it is used as the determination of a local minimum point. By traversing the time-series signal data of each feature sequence of the image to be decomposed, all local maximum points and local minimum points are identified, reflecting the fluctuation characteristics of the feature sequence signal of the image to be decomposed. B3. Fit the extracted sets of local maxima and local minima respectively to obtain two continuous envelopes. The upper envelope covers all local maxima and the lower envelope covers all local minima. During the fitting process, ensure that the envelopes completely cover all current data points. Calculate the arithmetic mean of the upper and lower envelope values ​​at the same time point to obtain the mean of the envelope at the same time point. Construct a continuous mean curve from the mean values ​​of the envelopes at all time points. B4. Subtract the mean value of the mean curve from each time point of the image feature sequence to be decomposed to obtain the corresponding candidate image feature signal sequence. Verify the candidate image feature signal sequence based on the intrinsic mode function to obtain the effective intrinsic mode function components. The verification logic includes that the absolute value of the mean of the upper and lower envelopes at all time points is less than or equal to 0.05 times the maximum signal amplitude, that is, the upper and lower envelopes of the candidate sequence are approximately symmetrical about the time axis, and the number of zero crossings in the candidate sequence is equal to or differs from the number of extreme points by no more than 1. The criterion for zero crossing is that the signal values ​​at adjacent time points have opposite signs, thus completing the decomposition of the laser speckle image feature sequence.

[0026] C1. After obtaining the effective intrinsic mode function components, the transformed signal sequence, i.e., the orthogonal component, is obtained by integral operation based on the characteristic value of each time point of the effective intrinsic mode function components. The orthogonal component is perfectly aligned with the original effective intrinsic mode function components in time, and the phase difference between the two is always 90 degrees. The effective intrinsic mode function components are taken as the real part and the transformed signal sequence is taken as the imaginary part, and the combination is used to form an analytical signal with a set of complex data of real part value and imaginary part value corresponding to each time point. C2. Squaring the real and imaginary parts of the analytic signal at each time point, then adding the two squaring results, and finally taking the square root of the sum to obtain the instantaneous amplitude. The instantaneous amplitude is the amplitude intensity of the analytic signal at the current time point, which directly reflects the dynamic change intensity of the effective intrinsic mode function components. C3. For the analytical signal at each time point, the phase angle is obtained by arctangent operation. That is, with the imaginary part as the numerator and the real part as the denominator, the phase angle difference between the current time point and the previous adjacent time point is calculated using the adjacent time point difference method. Then, the difference is divided by the time interval between the two time points to obtain the instantaneous frequency of the current time point. Furthermore, the instantaneous frequency value at the current time point is subtracted from the instantaneous frequency value at the previous adjacent time point to obtain the frequency difference. The frequency difference is then divided by the time interval to obtain the instantaneous frequency change rate at the current time point. A number of consecutive time points that can cover a sufficient temporal range are taken as a sliding window. The instantaneous amplitude values ​​of the time points within the window are added together and then divided by the number of time points to obtain the mean value within the window. The difference between the instantaneous amplitude value and the mean value at each time point within the window is calculated sequentially. Each difference is squared and summed. The sum is then divided by the number of time points minus 1. Finally, the square root of the result is taken to obtain the instantaneous amplitude fluctuation amplitude at the current time point. The phase angles of all time points of the effective intrinsic mode function components corresponding to the feature sequence of the image to be decomposed are calculated. At the same time point, the phase angle difference between any two effective components is calculated to obtain the instantaneous phase difference at that time point. The phase difference is calculated by pairwise pairing as described above to obtain an instantaneous phase difference time series curve and phase angle difference value corresponding to each pairing combination.

[0027] D1 and S102 use the sampling frequency of the dual-spectral system as the time interval to calculate the instantaneous frequency change rate and instantaneous amplitude fluctuation amplitude, and obtain the instantaneous phase difference. With the time point as the horizontal axis and the instantaneous frequency change rate, instantaneous phase difference, and instantaneous amplitude fluctuation amplitude as the vertical axis, a microscale dynamic feature matrix with a third feature dimension is constructed.

[0028] The third feature dimension of this microscale dynamic feature matrix includes the instantaneous frequency change rate of the effective intrinsic mode function corresponding to the equivalent diameter of speckle particles, the instantaneous phase difference of the effective intrinsic mode function corresponding to the equivalent diameter of speckle particles, the instantaneous amplitude fluctuation of the effective intrinsic mode function corresponding to the gray mean, the instantaneous frequency change rate of the effective intrinsic mode function corresponding to the gray variance, the instantaneous phase difference of the effective intrinsic mode function corresponding to the gray variance, and the instantaneous amplitude fluctuation of the effective intrinsic mode function corresponding to the texture entropy.

[0029] D2. After obtaining the microscale dynamic feature matrix and the molecular scale feature matrix, Pearson correlation analysis is performed on a time-by-time basis to pair all the first feature dimensions and the third feature dimensions at the same time point. The first feature dimension represents the molecular scale feature, and the third feature dimension represents the microscale feature. Redundant features are removed by calculating the second correlation coefficient. Referring to the scheme disclosed in the first correlation coefficient in the above steps, after completing the time-by-time feature association, a preliminary fused feature matrix is ​​established.

[0030] D3. By calculating the pairwise Euclidean distances of the feature vectors at all time points in the preliminary fused feature matrix, a minimum spanning tree is constructed based on the Euclidean distance. The minimum spanning tree is an algorithm for constructing a connected graph, ensuring that all vertices, i.e., the feature vectors at all time points, are connected and that the total edge length, i.e., the feature vector distance, is minimized. This is used to determine a reasonable distance threshold for topology modeling, and the distance threshold is set to 1.2 times the maximum edge length in the minimum spanning tree. This multiple ensures that the core connectivity of the feature vectors is preserved intact in the topology structure, while avoiding fragmentation of the topology structure due to an excessively small threshold. A sliding time window is used to divide the full-time-series feature vectors in the preliminary fused feature matrix to ensure sufficient coverage of the time series to characterize macroscopic changes. For each feature vector within the sliding time window, all other vertices with distances less than the distance threshold are connected and edges are generated to construct a topological complex. For each sliding time window, topological invariants, including the Betti number, Euler feature number, and persistent homology barcode parameters, are extracted from the topological complex.

[0031] The topological complex is used to establish the relationship between topological structure and time series. The Betti number, Euler feature number and persistent homology barcode parameters extracted from the topological complex represent the core invariant of the connectivity of the topological space, the comprehensive invariant of the overall structure of the topological space, and the invariant of the appearance and disappearance process of topological features, respectively. The key extraction methods are the 0th and 1st order Betti numbers. The 0th order Betti number corresponds to the number of connected components in the complex, and the 1st order Betti number corresponds to the number of pores in the complex. These are extracted by calculating the homology group of the complex. The Euler eigenvalue is calculated based on the simplex dimension, which is the number of vertices minus the number of edges plus the number of faces minus the number of volumes. It reflects the stage transition of the coagulation process, corresponding to the completion of the transformation of milk from liquid to gel state. The 0th and 1st order Betti numbers represent the degree of aggregation of microscopic gel particles during the coagulation process. The decrease in the number of connected components indicates an enhanced particle aggregation trend, reflecting the initial formation of the gel network. The number of pores first increases and then decreases, corresponding to the process of gel particles aggregating to form a network. The persistent coherence barcode first gradually increases the distance threshold from 0 to twice the set threshold, records the threshold when the feature appears and the threshold when the feature disappears for each Betti number, forming a barcode with the threshold as the horizontal axis and the duration as the vertical axis. The barcode length reflects the stability of the topological feature; the longer the length, the smoother the change in the coagulation process. The amplitude corresponds to the size of the Betti number, reflecting the strength of the topological feature.

[0032] D4. After obtaining the topological invariants including the Betti number, Euler eigenvalue, and persistent cohomology barcode parameters, a macroscopic trend feature matrix is ​​constructed with the midpoint of the time corresponding to the sliding window as the horizontal axis and the topological invariants including the Betti number, Euler eigenvalue, and persistent cohomology barcode parameters as the vertical axis. An adaptive weight fusion strategy is used to fuse the microscale dynamic feature matrix, the macroscopic trend feature matrix, and the molecular scale feature matrix. Through the particle swarm optimization algorithm, with the objective function of maximizing the goodness of fit between the fused features and the degree of condensation, the optimized weight combination is finally output, and a cross-scale condensation feature matrix is ​​formed by weighted summation at each time point.

[0033] S103. Extract the covariance matrix from the cross-scale condensation feature matrix, perform eigenvalue decomposition on the extracted covariance matrix, extract the largest eigenvalue and then normalize it to obtain the comprehensive response value and the rate of change of the comprehensive response value, and determine the adjustment direction based on the rate of change of the comprehensive response value.

[0034] The covariance matrix of the obtained cross-scale condensation feature matrix is ​​calculated, and the covariance matrix is ​​decomposed into eigenvalues ​​to extract the maximum eigenvalue. The maximum eigenvalue reflects the core information contribution of the covariance matrix and can characterize the overall dominant change trend of cross-scale fusion features. After normalization, the comprehensive response value is obtained. The comprehensive response value obtained by the dual-spectral system in the previous and next cycles is used to calculate the rate of change of the comprehensive response value. If the rate of change of the comprehensive response value is greater than zero, the adjustment is reversed; if the comprehensive response value is less than zero, the adjustment is forward; if the comprehensive response value is equal to zero, the operating data of the condensation reaction component is kept unchanged, and the comprehensive response value is used as the adjustment amount of the operating data of the reaction component.

Claims

1. A method for monitoring and controlling cheese processing based on dual-spectrum and spectral analysis coagulation process, characterized in that, The method includes: S101. Based on the terahertz dynamic evolution signal and speckle texture dynamic evolution signal obtained by the dual-spectrum system, time axis alignment and preprocessing are performed to obtain the terahertz spectral signal sequence and the laser speckle image feature sequence. S102. The spectral signal intensity of the terahertz spectral signal sequence is enhanced by coherent anti-Stokes Raman scattering spectroscopy. The enhanced terahertz spectral signal sequence is then fitted with Lorentz to construct a molecular-scale feature matrix. A nonlinear signal analysis strategy is embedded to decompose the laser speckle image feature sequence into effective intrinsic mode function components. After transforming the effective intrinsic mode function components, a microscale dynamic feature matrix is ​​constructed. The microscale dynamic feature matrix and the molecular-scale feature matrix are correlated at each time point to establish a preliminary fused feature matrix. Then, a topological data analysis mechanism is introduced to perform topological modeling. Topological invariants are extracted under a sliding time window. The extracted topological invariants constitute a macroscopic trend feature matrix. The microscale dynamic feature matrix, the macroscopic trend feature matrix, and the molecular-scale feature matrix are fused to finally obtain a cross-scale condensation feature matrix. S103. Extract the covariance matrix from the cross-scale condensation feature matrix, perform eigenvalue decomposition on the extracted covariance matrix, extract the largest eigenvalue and then normalize it to obtain the comprehensive response value and the rate of change of the comprehensive response value, and determine the adjustment direction based on the rate of change of the comprehensive response value.

2. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 1, characterized in that, S101 integrates a dual-spectrum system, including a terahertz time-domain spectral emission module, a terahertz receiving module, a laser speckle interferometric emission module, and an image acquisition module, onto an adjustable bracket above the cheese processing container. The detection directions of the terahertz time-domain spectral emission module and the laser speckle interferometric emission module are perpendicular to the milk surface, and the detection areas of the dual-spectrum system completely overlap. By controlling the start-up and data acquisition timing of the dual-spectrum system, time axis alignment is completed. After baseline correction and noise suppression, a time axis-aligned terahertz spectral signal sequence and a laser speckle image feature sequence are obtained.

3. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 1, characterized in that, After acquiring the terahertz spectral signal sequence, S102 uses coherent anti-Stokes Raman scattering spectral analysis to detect the surface of milk with pump light and Stokes light to obtain an enhanced terahertz spectral signal sequence. The characteristic peak amplitude of the terahertz spectral signal sequence is divided by pi and multiplied by the square of half the full width at half maximum (FWHM) of the characteristic peak, and the sum of the squares of the center frequencies of the characteristic peaks is subtracted to obtain the spectral signal intensity value after Lorentz fitting. Based on the fitted curve obtained by fitting, the characteristic peak parameters of protein amide bonds, including characteristic peak amplitude, FWHM, and peak position offset, and the characteristic peak parameters of phospholipid molecules at the fat globule interface, including the rate of change of characteristic peak amplitude, are extracted. Using the time point aligned with the S101 time axis as the horizontal axis and the first feature dimension composed of peak amplitude, full width at half maximum (FWHM), peak position offset, and rate of change of molecular peak amplitude at the fat globule interface as the vertical axis, a preliminary molecular-scale feature matrix is ​​constructed.

4. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 3, characterized in that, After obtaining the preliminary molecular-scale feature matrix, Pearson correlation analysis is used to pair each of the first feature dimensions at the same time point to form a first correlation array. The covariance of the time series signal data corresponding to the two first feature dimensions is divided by the product of the standard deviations of the two first feature dimensions to obtain the first correlation coefficient of each first correlation array. Any first feature dimension at the same time point is removed from the first correlation array whose first correlation coefficient exceeds twice the mean of the first correlation coefficient. Finally, the molecular-scale feature matrix is ​​constructed.

5. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 1, characterized in that, After constructing the molecular-scale feature matrix, S102 uses the laser speckle image feature sequence, which includes the equivalent diameter of speckle particles, gray mean, gray variance, and texture entropy as the input of the nonlinear signal analysis strategy. After feature extraction, the image feature sequence to be decomposed, which includes the equivalent diameter of speckle particles, gray mean, gray variance, and texture entropy, is obtained respectively. In the feature sequence of the image to be decomposed, if the feature value of the second feature dimension at the same time point is greater than the feature value of the second feature dimension at each of the two adjacent time points, it is used as the determination of a local maximum point. If the feature value of the second feature dimension at the same time point is less than the feature value of the second feature dimension at each of the two adjacent time points, it is used as the determination of a local minimum point. By traversing the time-series signal data of each feature sequence of the image to be decomposed, all local maximum points and local minimum points are identified. The extracted sets of local maxima and local minima are fitted to obtain two continuous envelopes. The values ​​of the two envelopes at the same time point are arithmetically averaged to obtain the mean value of the envelope at the same time point. The mean values ​​of the envelopes at all time points are used to form a continuous mean curve. Subtract the mean value of the mean curve from each time point of the image feature sequence to be decomposed to obtain the corresponding candidate image feature signal sequence. Verify the candidate image feature signal sequence based on the intrinsic mode function to obtain the effective intrinsic mode function components, thus completing the decomposition of the laser speckle image feature sequence.

6. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 5, characterized in that, After obtaining the effective intrinsic mode function components, the transformed signal sequence is obtained by integral operation based on the characteristic values ​​of each time point of the effective intrinsic mode function components. The effective intrinsic mode function components are taken as the real part and the transformed signal sequence is taken as the imaginary part, and the combination is used to form an analytical signal with a set of complex data of real part values ​​and imaginary part values ​​corresponding to each time point. The real and imaginary parts of the analytical signal at each time point are squared respectively, the two squared results are added together, and finally the square root of the sum is taken to obtain the instantaneous amplitude. For the analytical signal at each time point, the phase angle is obtained by arctangent operation. The difference between the phase angles of the current time point and the previous adjacent time point is calculated by using the adjacent time point difference method. The difference is then divided by the time interval between the two time points to obtain the instantaneous frequency of the current time point.

7. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 6, characterized in that, S102 uses the sampling frequency of the dual-spectral system as the time interval to calculate the instantaneous frequency change rate and instantaneous amplitude fluctuation, and obtains the instantaneous phase difference. With the time point as the horizontal axis and the instantaneous frequency change rate, instantaneous phase difference, and instantaneous amplitude fluctuation as the vertical axis, a microscale dynamic feature matrix with a third feature dimension is constructed. After obtaining the microscale dynamic feature matrix and the molecular scale feature matrix, Pearson correlation analysis is performed on a time-by-time basis to pair all the first and third feature dimensions at the same time point. Redundant features are removed by calculating the second correlation coefficient. After completing the time-by-time feature association, a preliminary fusion feature matrix is ​​established.

8. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 1, characterized in that, S102 calculates the pairwise Euclidean distances between the feature vectors at all time points in the preliminary fusion feature matrix, constructs a minimum spanning tree based on the Euclidean distances, and sets the distance threshold to 1.2 times the length of the maximum edge in the minimum spanning tree. A sliding time window is used to divide the full-time feature vectors in the preliminary fusion feature matrix. For each feature vector in the sliding time window, all other vertices with distances less than the distance threshold are connected and edges are generated to construct a topological complex. For each sliding time window, topological invariants including the Betti number, Euler feature number, and persistent homology barcode parameters are extracted from the topological complex.

9. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 8, characterized in that, After obtaining the topological invariants including the Betti number, Euler eigenvalue, and persistent cohomology barcode parameters, a macroscopic trend feature matrix is ​​constructed with the midpoint of the time corresponding to the sliding window as the horizontal axis and the topological invariants including the Betti number, Euler eigenvalue, and persistent cohomology barcode parameters as the vertical axis. An adaptive weight fusion strategy is used to fuse the microscale dynamic feature matrix, the macroscopic trend feature matrix, and the molecular scale feature matrix. Through the particle swarm optimization algorithm, with the objective function of maximizing the goodness of fit between the fused features and the degree of condensation, the optimized weight combination is finally output, and a cross-scale condensation feature matrix is ​​formed by a fusion logic of weighted summation at each time point.

10. The cheese processing monitoring and control method based on dual-spectrum and spectral analysis coagulation process according to claim 1, characterized in that, S103 calculates the covariance matrix of the cross-scale condensation characteristic matrix, performs eigenvalue decomposition on the covariance matrix, extracts the maximum eigenvalue, and obtains the comprehensive response value after normalization. The comprehensive response value obtained by the dual-spectral system in the previous and next cycles is used to calculate the rate of change of the comprehensive response value. If the rate of change of the comprehensive response value is greater than zero, the adjustment is reversed; if the comprehensive response value is less than zero, the adjustment is forward; if the comprehensive response value is equal to zero, the operating data of the condensation reaction component is kept unchanged, and the comprehensive response value is used as the adjustment amount of the operating data of the reaction component.

Images