1. Introduction
Food products are routinely analyzed for a variety of reasons, including to comply with regulatory and labeling requirements, to meet product specifications, to ensure food safety, and for research, development, and innovation purposes. Quality testing is therefore an integral part of food production.
Grab sampling means taking a small portion of a particular
lot (e.g., a batch or production cycle) to the laboratory for testing. The sample is assumed to be representative of the whole lot. Different sampling regimes such as random, stratified, systematic or composite sampling can be applied to minimize the sampling error [
1,
2]. Nevertheless, grab sampling of heterogeneous food products is challenging and often leads to substantial sampling biases. Other challenges with grab sampling are time and resources. The time delay between sampling and analysis results limits the possibility for using the measurements in real-time process monitoring and control systems.
In-line measurement systems are viable alternatives to grab sampling for many quality attributes. Sensors placed directly in the process stream, either over a conveyor belt or in/onto a pipe or processing tank can automatically measure a large proportion of the lot, in real time, and therefore overcome many of the challenges of grab sampling.
Near-infrared spectroscopy (NIRS) is a technology that is based on illuminating a sample and analyzing the amount of transmitted or reflected light. This information is used to measure the properties of a sample. It is well known that NIRS is suitable for the in-line measurement of a range of food quality attributes related to chemical composition. The technology is widely used in the food industry, typically for the estimation of fat, water, protein, and carbohydrates in products such as meat, fish, cereals, and fruits [
3]. Recent reported examples are industrial monitoring of dry matter in large batches of potatoes [
4] and in-line quality sorting of chicken fillets [
5]. NIR spectroscopy is a well-established method for the rapid determination of fat, protein, and moisture in cheese [
6], and NIRS is widely used in the dairy industry [
7]. Several studies have shown that dry matter and fat can be measured in cheese in the laboratory [
8,
9,
10], while in a modern dairy there is a need to measure these features in-line on large cheese blocks weighing e.g., 20 kg. Eskildsen et al. [
11] showed that it is possible to estimate the average dry matter and fat content of blocks of Swiss cheese with satisfactory accuracy by using in-line scanning NIRS in reflectance mode. Good results were obtained on cheese blocks both before and after pressing. This is possible when the surface chemistry of the cheese is representative of the interior.
Continuous measurements of quality attributes of food products define time series and different tools for analyzing the process variation in more detail becomes available. One such tool is the power spectral density (PSD), which provides a decomposition of signal power (variation) in the frequency domain and is particularly useful when periodic components exist. Sometimes frequency components may be related to distinct causes. In this article, we used the PSD on predictions from NIR spectra.
There are numerous examples of this. For instance, in [
12], the authors studied the flow of particles suspended in a gas through a pipe using electrical capacitance tomography at two imaging planes. They used PSD to analyze aspects of the flow and argued that the probable cause of one peak in the PSD was due to a slip/stick condition at the particle inlet into the pipe. They also linked a second peak to periodic variations in the radial distribution of the particle flow.
In [
13], experiments were carried out to relate the information in the PSD to process parameters in a laboratory-scale wet mill. They first divided the frequency domain into three bands and for each showed probable causes: the low-frequency band was linked to the mill, the mid-frequency linked to balls–mill shell lining interactions and the high-frequency range linked to ball–ball interactions. Further, via their experimental design, they related this information to parameters defining the operation of the mill which was subsequently used to formulate soft sensors. Similar uses of the power spectral density can be found in [
14], where it was used to distinguish between possible physical processes in liquid-solid pipe flows; in [
15] it was used as one part of a system for diagnosing stator faults in induction motors; and in [
16] where three signal sources for diagnostics were compared using the PSD.
To the best of the authors’ knowledge, no publication has so far used the power spectral density as a diagnostic tool in the food industry.
In [
17], a review of the role of NIRS for food quality analysis shows a list of achievements including but not limited to: milk and yoghurt attributes in the dairy industry; monitoring fermentation processes in wine and beer production; the drying of fruit and vegetables post-harvest; and dry ingredient mixing. However, the authors concluded that much had been performed to demonstrate the feasibility of applications in laboratory conditions but identify a gap regarding actual implementations at the industrial scale.
The aim of this paper is therefore to document some of the benefits of in-line measurements at the industrial scale, including the advantages over traditional grab sampling whilst also showcasing the possibilities of using PSD for continuous in-line measurements as a diagnostic tool. We achieve this by presenting an actual case study from cheese production, where in-line NIRS was implemented to replace traditional lab measurements of dry matter in Gouda-type cheese. In doing so, we will focus more on the temporal dimension of the approach rather than conducting an in-depth analysis of the spectral information.
2. Materials and Methods
This paper is based on data from a large cheese dairy that produces several types of semi-hard cheeses. Cheese making is a complex biochemical process, and its industrial production involves a series of highly controlled processing steps: a vat is filled with pasteurized milk with a standardized protein/fat ratio; starter culture and rennet are added thereby coagulating the milk and forming a curd; the curd is then cut and heated inside the vat, after which most of the whey is drained off. The remaining curd/whey mixture is transferred to a “casomatic” (Tetra Pak
® Tebel, Netherlands), which consists of vertical columns where the curd particles fuse together to form solid cheese at the bottom. Cheese blocks are then sliced off as they exit the column and are transported on a conveyor belt to a press. They are then placed in a brine bath and stored in a controlled environment until they reach the desired maturity.
Figure 1 shows a schematic overview of the process with emphasis on the steps that are most relevant in this work.
Among the most important quality attributes of semi-hard cheeses is the dry matter content as it affects ripening and eating quality as well as profitability. Even if the process is highly standardized and controlled, the dry matter may vary substantially from batch to batch due to uncontrollable variations in raw milk quality, starter culture and rennet.
2.1. Data Collection
Data were collected during two time periods: (i) during test installations of an in-line NIRS system in 2019 and (ii) after full implementation of the same system in 2022. The dairy produces several different cheese types, but this study is limited to the one with the largest production volume.
In May–July 2019, a test installation of the NIRS instrument was run for a period of 18 production days at control point 1, CP1, and 11 production days at CP2 (see
Figure 1). The same instrument was used at both time points. This test included data from the production of (i) 36,794 cheese blocks from 404 cheese vat batches at CP1, and (ii) 22,435 cheese blocks from 245 cheese vat batches at CP2.
In 2022, a permanent implementation of the same NIRS system was installed at CP1. About one week of continuous measurement from the permanent installation was included in this work, corresponding to 8327 cheese blocks from 94 batches.
2.2. Lab Measurements of Dry Matter
Grab samples were manually collected at a the control point directly after the brine bath, denoted as CP3 in
Figure 1. Samples were routinely taken from approximately every fourth batch, by cutting a cross-section of 2 cm diameter from the center of one cheese block. According to the dairy’s protocol, samples were taken from a cheese block in the middle of the batch to account for a within-batch gradient caused by the casomatic (see
Section 2.4). The dry matter was measured on a benchtop NIR system (FoodScan
TM, FOSS Analytics, Hillerød, Denmark) using commercial calibration models.
2.3. In-Line Measurements of Dry Matter
The NIR spectrometer (DA 7440, Perten Instruments, Hägersten, Sweden), designed for in-line measurements on conveyor belts, was placed approximately 25 cm above the moving cheeses at either control points CP1 or CP2 (see
Figure 1). These points were selected based on previously conducted trials in the dairy’s pilot plant [
11]. The cheese blocks were transported in plastic containers with each block weighing approximately 20 kg. Spectra in the wavelength range 950–1650 nm were acquired continuously with a time resolution of 4 spectra/second. Spectra were aggregated to a single spectrum per cheese and pre-processed using the standard normal variate (SNV) approach.
Partial Least Squares (PLS) calibration models for dry matter were developed individually for CP1 and CP2, with both models calibrated to predict dry matter at CP3. This meant estimates from all three control points were directly comparable, even if the actual dry matter varied between them. Average calibration spectra at CP1 and CP2 are shown in
Figure 2 with systematic differences between the two control points (which lead to slightly different models).
For the test installation in 2019, provisional calibration models based on approximately 50 and 30 samples were developed for CP1 and CP2, respectively. Samples were split in an 80–20% ratio between training and test sets. A 10-fold cross validation was used on the training set to select the number of components. The models obtained RMSEP values of 0.63% and 0.38% dry matter at CP1 and CP2, respectively, based on calibration samples spanning approximately 2% dry matter. These models are not very accurate as their RPD values were 1.3 and 2.3, respectively, but they still capture relative differences between cheeses that allow for greater understanding of some sources of variation in the process. In 2022, the calibration model for CP1 was significantly improved by the instrument vendor, by adding more samples and refining the model. The updated model spanned a range of 3.5% dry matter and obtained a RMSECV of 0.36%.
2.4. Data Handling and Integration
During the test installation in 2019, the spectrometer was not integrated in the dairy’s Manufacturing Execution System (MES) and a three-step procedure was therefore developed to assign single spectra to cheese blocks and further to the corresponding cheese vat batch:
Aggregating spectra for each cheese block: The instrument was set up to provide an indication of valid spectra (corresponding to “cheese” and not background) based on a prior model. The temporal separation between two cheeses was used to aggregate the sequence of spectra into a single spectrum per cheese. This association of spectra to cheeses is not perfect, and errors occur leading to a variation in the length of batches in terms of the number of cheeses. It will be shown later that this complicates the decomposition.
Assigning cheese blocks to batches: This was accomplished based on the shape of the time series of dry matter estimates; the casomatic (
Figure 1) induced a gradient in the dry matter content in the resulting cheeses. This effect of increasing dry matter within the batches was well known and is illustrated in
Figure 3. Transitions between batches were found by considering the expected negative jump in dry matter values by approximately 1%. This was discovered using an appropriate smooth derivative filter in the time domain and by looking for local maxima at temporal intervals approximately 90 cheeses.
Linking the identified batches to MES: The batches from step 1 were matched with the MES system based on time stamps and the approximate time delays between the cheese vat and CP1/CP2. Data on other production parameters such as temperature, pH and durations of different processing steps in the cheese vat could then be exported from the MES and laboratory databases and then combined with the in-line measurements.
All steps were performed in a semi-automated way, meaning that the automatically obtained results were manually checked and corrected if necessary. Nonetheless, the procedure is prone to error, particularly in determining the transition between batches in step 2. This is illustrated in
Figure 3 (observe e.g., the transition between light blue and dark red batches at approximately 10:15).
In the full implementation in 2022, the spectrometer was integrated with the MES to ensure full traceability between single cheese blocks and cheese vat batches. However, missing values were detected whenever the cheese ID was not a contiguous sequence, or a new batch did not start with cheese number 1, or the length of the batch was not an integer number of 4 cheeses (the number of columns in the casomatic). These missing values have been imputed in an attempt to make the series as consistent as possible with respect to the casomatic column. The imputation was based on the known triangular form of the casomatic contribution and performed separately for each column: , where were estimated from the cheeses in the same batch and from the same column.
2.5. Time Series Trend Analysis
The in-line NIRS measurements formed a time series, and it was already known that the dry matter content of a single cheese was given by the cheese vat batch mean as well as its order in the casomatic step (within a batch), which gave the cyclic trends shown in
Figure 3. It was therefore natural to decompose the dry matter variation into components representing the batch mean, a linear casomatic trend and residuals.
The production ran continuously for periods, with pauses due to cleaning or production of other cheese types in between. These pauses were used to split the entire data set into subsets of continuous time series, denoted “sequences”. The time series analyses were performed separately on each sequence, and were assumed to be continuous.
2.5.1. Decomposition of the Time Series
The dry matter time series signal (
) was decomposed by first considering the overall mean (
μ). The centered sequence (
) was then split into batches where least squares models including a constant and linear term were fit to each batch separately. The batch offsets (
) were equated to the constant terms while the casomatic contributions (
) were associated to the linear components. However, due to the uncertainty related to batch transitions, only the central 80% of the samples were used for fitting the model. The remaining signal was denoted the residual (
). The full decomposition was therefore expressed as:
Note, the batch mean was the sum of the first two terms: .
2.5.2. Power Spectral Density
Stationary processes such as this can be described using their power spectral density (PSD), which is the Fourier transform of a signal’s autocorrelation function:
where
is the autocorrelation function of the signal
,
the time lag,
the Fourier transform and
the PSD as a function of frequency
. Equation (3) leads to the relationship between the power of the signal and the PSD [
18]:
which implies that the power associated with any frequency interval can be found by integrating the PSD over that interval. The
is hence the expected power per unit frequency. For instance, the casomatic signal with a period of say 92 cheeses would have a peak in its PSD spectrum at a (fundamental) frequency of
, along with components at multiples (
) of the fundamental frequency
. This set (sometimes with the exception of the fundamental) is collectively referred to as the harmonics. This is a property shared by all periodic signals and can be used to identify other periodic sources of influence on the process. A perfectly periodic signal has non-zero PSD values only at its harmonic frequencies, while signals with small variations in period or shape will give non-zero PSD values over regions around these frequencies.
The implementation of the PSD is based on the discrete signal (sampled in time) and using Welch’s method by averaging the estimates of different windows of the signal within a sequence [
19]:
where
indexes the window used for the distinct estimates
of the PSD based on individual, overlapping windows of length
and
is the discrete time Fourier transform (DTFT) of the time series
. The restriction on
is due to sampling and has units of cycles/sample. Also, to increase contrast in the PSD the
Hann window (
) was used with appropriate scaling (
) of the PSD estimates.
Estimates of time series components corresponding to quasi-periodic signals can be obtained based on filtering in the frequency domain with narrow passbands around its harmonic frequencies. It is worth noting that this procedure reliably separates components only when the different components have a disjoint support in the frequency domain, i.e., the peaks in the PSD do not overlap.
2.6. Comparison of Sampling Regimes by Simulation
The main aim of the dry matter measurements was to estimate the batch mean, and in-line NIR was expected to give more precise estimates than the traditional grab sampling regime. The estimates may have been affected by cyclic trends in the time series, and a simulation study was set up to compare different sampling regimes for various magnitudes of cyclic trend components. Our results showed that there are two cyclic trends present at CP2, caused by the casomatic and press steps (described in
Section 3.1). Data were therefore simulated to resemble the observed dry matter time series at CP2, since this control point was assumed to be similar to CP3, where the grab sampling was performed. A continuous time series of N = 920,000 cheese blocks (10,000 batches with 92 units in each batch) was constructed by summing independent components representing batch mean (
), casomatic contribution (
), press contribution (
) and noise (
):
The batch mean was constant over all 92 units in a batch. The casomatic effect was simulated as a linear cyclic effect with positive slope and of period 92, and the press effect as a linear cyclic effect with negative slope and period 30. The noise was independent and normally distributed.
Four different time series were simulated, with varying size of the casomatic and press components:
- -
Simulation 1: the press component was set to zero, mimicking a case where there was only one cyclic trend which was aligned with the batch.
- -
Simulation 2: the size of the press and casomatic components were similar to the observed data, i.e., the variation of the casomatic component was approximately 3.5 times the press component (see
Table 2).
- -
Simulation 3: the magnitude of the press component was set to be 3.5 times the casomatic component, to mimic a situation where the dominating cyclic trend was not aligned with batches.
- -
Simulation 4: the casomatic component was set to zero, mimicking a case where there was only one cyclic trend which was not aligned with the batch.
For each of these time series, four different sampling regimes were performed:
Grab random: the batch mean was represented by a randomly selected unit.
Grab mid precise: the batch mean was represented by one of the 2 middle units in a batch.
Grab mid sloppy: the batch mean was represented by one of the 20 middle units.
In-line: the batch mean was estimated as the average over all units in a batch.
The precision of the batch mean estimates were evaluated by mean squared error (MSE) and the coefficient of determination (R
2). The variations in estimated batch means were also calculated, for comparison with observed data. The four simulated data sets are provided in the
Supplementary Materials.
4. Discussion
In-line measurements have many advantages over traditional grab sampling, also beyond the case illustrated here. With regard to precision, the added value of in-line measurements is higher when the lot is more heterogeneous. Appropriately designed sampling regimes, such as stratified or composite sampling, may also give representative grab samples with precise estimates of within- and between-lot variations. However, such regimes are often labor intensive, and in-line NIRS will in many cases be a more viable method. It is also clear from the results section that the average over a batch in this case improved the precision of batch estimates in addition to reducing the impact of asynchronous, periodic processes (the press signal) as well as high-frequency contributions (e.g., the columns).
A higher temporal resolution greatly improves process understanding and makes it possible to assess the variation within and between e.g., batches, production days, raw material providers and operator shifts. Also, portable and robust devices that can be moved between different measurement points gives the opportunity to distinguish between sources of variation in different steps of the process. Such data and insights are important for identifying and prioritizing improvement areas, e.g., as part of
lean manufacturing, six sigma or
sustainable production methodologies [
20]. The details that such continuous measurements provide may also improve the sensitivity of detecting faults or anomalous production conditions and the capacity to localize the source of such problems.
In the case study presented here, the use of in-line NIRS at two locations in the process has opened the door for a time-frequency and power spectral density analysis of the process. The differences observed along the production process have enabled an attribution of variation that would otherwise have been more difficult to document. This analysis has also identified components that were not previously known: the press’ and columns’ contributions to dry matter variation.
The success of a time-frequency analysis of the data in this case study has relied on the presence of periodic contributions as periodic signal components are very well distinguished: observe the spikes in
Figure 6 as opposed to the time-domain signal in
Figure 8 where the general slope is visible (casomatic) but other components are less so. This analysis allowed for the identification of a significant contribution from the casomatic based on the fundamental frequency and knowledge of the number of cheeses corresponding to a batch. For the same reasons, a second important contribution could be identified with a fundamental frequency of 30 cheeses, which in this case could be attributed to the press, both because of knowledge of the press (30 cheeses per arm), but also from comparing measurements at two locations. A third periodic component was also observed, and a probable association to the columns of the casomatic was suggested. Finally, it was possible to estimate a flat, background noise level. Such a flat spectrum is related to uncorrelated processes, meaning the part of the next dry matter measurement is uncorrelated to all other measurements. This noise level is best understood in terms of measurement uncertainty (or error) as it is unlikely that the process itself adds a unique, independent contribution to each cheese.
Having obtained the decomposition, partly based on prior knowledge and partly on observations in the time-frequency analysis, it was possible to quantify the sources of variation. This quantification showed the relative contributions of the components in the model, see
Table 1. This quantification, along with knowledge of its variability may allow for more precise detection of when an anomaly is occurring, e.g., small changes in the casomatic would be difficult to detect in the total dry matter signal given that the batch offset and residuals variations are large.
A second important point relates to the casomatic signal. Increasing dry matter through a batch is directly related to the position of the corresponding block in the columns of the casomatic and this relationship is basically due to a physical phenomenon. Therefore, while the signal induced variation in dry matter, it basically added the same variation every time as its cause was stable. This can be observed in how little this variance differed between “sequences” (the parentheses in
Table 1). At the other end, the estimate of variation due to batch offsets differed significantly from one “sequence” to the next: there were likely upstream causes that induced these differences. The fact that there were such differences suggested that its causes were unstable and may in part be open to optimization, after all there were sequences with low variance.
Having distinguished sources of variation, it may be desirable to reduce these. An obvious target for reduction would be the casomatic as it was found to be the processing step responsible for most of the within-batch variation. It might not be possible to remove this variation without changing the processing technology, but knowledge about it can give rise to mitigation strategies such as re-ordering of the cheese blocks in subsequent processing steps or sorting cheeses to different storage conditions. Another major contribution to variance is the batch offset. While its variance and instability has been quantified, further detective work would be necessary to relate this variation to possible causes, of which raw material variation and differences in process settings are usual suspects.
Among the advantages of time-frequency analysis mentioned above was the ability to distinguish signals that are periodic with different fundamental frequencies. Such signals are expressed as a series of spikes in the frequency domain, such as those observed in
Figure 6. However, time-frequency analysis is not limited to this category of signals: many processes occupy wide bands in the frequency domain. When describing signals as “smooth”, this normally translates into signals with a limited support in the frequency domain: there is little power above a certain frequency. In other words, the “speed” of the process is related to the highest significant frequency component, and variation e.g., due to different work shifts (approximately changing every 8 h) would occupy a larger band than seasonal variations (months).
This paper describes the retrospective analysis of in-line measurements. An even bigger potential of in-line NIRS lies in the real-time operationalization of the process measurements, usually coupled with results from retrospective analysis. Such applications may for instance be anomaly/fault detection, real-time process adjustments and control, and the sorting of products into different quality categories.