Multi-sensor data fusion and parallel factor analysis reveals kinetics of wood weathering

Understanding mechanisms of materials deterioration during service life is fundamental for their confident use in the building sector. This work presents analysis of time series of data related to wood weathering acquired at three scales (molecular, microscopic, macroscopic) with different sensors. By using several complementary techniques, the material description is precise and complete; however, the data provided by multiple equipment are often not directly comparable due to different resolution, sensitivity and/or data format. This paper presents an alternative approach for multi-sensor data fusion and modelling of the deterioration processes by means of PARAFAC model. Time series data generated within this research were arranged in a data cube of dimensions samples × sensors × measuring time. The original protocol for data fusion as well as novel meta parameters, such as cumulative nested biplot, was proposed and tested. It was possible to successfully differentiate weathering trends of diverse materials on the basis of the NIR spectra and selected surface appearance indicators. A unique advantage for such visualization of the PARAFAC model output is the possibility of straightforward comparison of the degradation kinetics and deterioration trends simultaneously for all tested materials.


Introduction
Weathering is a natural process occurring to all materials exposed to environmental conditions. Depending on the local climate, material type or construction details, the progress of deterioration, related to functional and aesthetical aspects, may differ significantly [1]. Weathering deterioration is especially relevant for wood and other materials with biological origin that are commonly perceived as not resistant to natural weathering. However, this is not exactly true, and for that reason, it is extremely important to properly understand their weathering mechanisms and predict the extent of material deterioration in order to assure satisfactory performance of bio-based products [2].
Recent advancements in diverse scientific fields has led to development of several instruments and methodologies for assessment and monitoring of material properties. The latest trend in monitoring material status and performance consists of using multiple sensors simultaneously, which proved to be more effective with respect to a single sensor approach, due to more accurate representation of reality [3]. Selection of optimal sensors, measurement strategy, signal processing and interpretation of results is challenging and demanding, especially when considering the complex and heterogenous surface of wood [4]. The fusion of complementary information available from different sensors will yield increased knowledge of the investigated phenomena, which is not attainable separately by any of the single sources. Moreover, it is expected that combining data from complementary analyses will improve the performance of statistical models. However, such multi-sensor monitoring introduces new issues and challenges in the data analysis phase, where the strategy for merging different sources of information is fundamental [5][6][7].
Data acquired by diverse sensors usually rely on the measurement of different physical phenomena that, in the majority of cases, are not directly correlated to each other or respond with different kinetics to material alterations. Thus, both data integration as well as interpretation of obtained results are challenging. Some sensors may be more suitable for particular applications than others, due to higher sensitivity in certain measurement conditions, better resolution or optimized scanning time. Finally, the collected data bring shared (correlated) as well as distinctive information, which makes it necessary to interpret them altogether [8,9].
In general, multi-sensors measurement generates a set of data collected by different platforms and available for the entire set of samples, with each set usually referred to as a data block. Multi-block analysis, and more generally data fusion techniques, allows joint extraction of information from different blocks, conveying it in a single model that is more efficient than creating individual models. Data fusion approaches can be distinguished according to the level at which data are fused in low-(fusion of raw data), mid-(fusion of features extracted by each single data block) and high-level (fusion of results/decisions derived by each data block) [10]. Multi-block methods usually operate at low-level [5,11]. This paper presents the results regarding in service performance of modified wood during natural weathering test. Investigated bio-based materials were characterized in a preceding research [12], during and after degradation by biotic and abiotic agents in order to provide experimental data to be used for better understanding of their performance/degradation as a function of time and/or weather dose. A multi-sensors measurement chain allowed the acquisition of properties at different scales (molecular, microscopic, macroscopic). Time series data generated within this research refers to multi-sensor measurements and are in the form of a three-way array. In that case, samples are characterized by several sets of variables acquired at different time, thus can be arranged in a data cube of dimensions samples × sensors × measuring time. A data analysis pipeline consisting of low-level data fusion and multiway data decomposition was applied in order to deal with a three-way multi-block. Moreover, an effective post-processing of the obtained fused model in order to better depict the weathering trends affecting wood samples that were subjected to diverse bulk modification procedures was proposed.

Experimental samples
Radiata pine (Pinus radiata D. Don) samples representing six different commercially available modification processes were selected for the demonstration. Weathering performance of these materials was compared with untreated wood of Scots pine (Pinus sylvestris L.), being usually considered in standards as a reference material (#1) [13]. Each material was represented by three replicated samples. The wood modification processes included: thermal treatment (#2), thermal treatment with penetrating oil (#3), thermal treatment with silicate treatment (#4), thermal treatment with coating (#5), furfurylation (#6) and acetylation (#7). All investigated modified materials are commercially available on the market and were investigated for their service life performance and aesthetical deterioration [12].

Weathering test
Natural weathering tests were performed in San Michele, Italy (46 • 11 ′ 15 ′′ N, 11 • 08 ′ 00 ′′ E), with the objective of providing reference data for simulation of the bio-based materials' performance as a function of exposure time. Samples were exposed on the vertical stand oriented South. The weathering experiment was carried out for 12 months and started in March 2017. Twelve samples for each type of modified wood were exposed at the beginning of the test. Three replicas wood samples from each set were dismounted from the stand after three, six, nine and 12 months. In addition, three intact samples were included in the collection that correspond to no-weathering, i.e. reference samples. All samples after collection from the stand were stored in a dark and climatized room (20 • C, 60% relative humidity) to terminate any progress of deterioration. The alteration of wood surface appearance is presented in Fig. 1, where all investigated materials are shown at the different time of exposure.

Multi-sensor characterization of weathered materials
All the weathered wood samples were assessed in order to highlight the appearance changes as well as to quantify the deterioration progress. Assessment was performed after conclusion of the natural weathering experiment and conditioning of the wood samples. The whole batch of all samples at diverse stages of deterioration was measured at once after conclusion of the weathering campaign in a random sequence to avoid any time-related or methodological bias variance in the data. Materials characterization included measurement of the NIR and colour spectra and derived parameters, wood surface imaging as well as gloss.

CIE L*a*b* colour
Changes in colour were assessed by a spectrometer following the CIE Fig. 1. Change of modified wood appearance along one-year exposure to natural weathering.
L*a*b* system, where colour is expressed with three parameters: L* (lightness), a* (red-green tone) and b* (yellow-blue tone). CIE L*a*b* colours were measured using a MicroFlash 200D spectrophotometer (DataColor Int, Lawrenceville, USA). The selected illuminant was D65 and viewer angle was 10 • . All specimens were measured on ten different spots randomly selected over the weathered surface.

Colour RGB and HSL
Experimental samples were scanned with an office scanner, i.e. HP Scanjet G2710, and colour texture was represented as RGB and HSI images. A set of descriptors linked to early-and latewood components were derived from each RGB/HSI image [14,15]. The frequency histograms of the intensity values corresponding to each colour channel (R, G and B) as well as to H, S and I were analysed to identify values corresponding to the tails of the colour distribution. These tails correspond to extreme wood tissues that form the ring structure of the standing tree, with the bright and low density wood originating at the beginning of the vegetation season (earlywood), which is contrasted to dark and high-density latewood originating in autumn. As a result, twelve descriptors were obtained corresponding to R early , G early , B early , H early , S early , I erly, R late , G late , B late , H late , S late , I late , respectively. A custom software for image analysis was implemented in LabView 2018 (National Instruments, Austin, USA) software package.

Gloss
The mode of light reflection from the surfaces was measured using a REFO60 (Dr. Lange, Düsseldorf, Germany) gloss meter with an incidence/reflectance angle of 60 • . Ten measurements were taken in five randomly selected locations over each specimen surface, following two light irradiance directions corresponding to along and across the fibres.

FT-NIR spectroscopy
Near infrared (NIR) spectra were collected with a Vector N-22 Fourier-transform NIR spectrometer produced by Bruker Optics GmbH (Ettlingen, Germany). The system was equipped with a fibre optic probe, and the measurement range was between 12,000 cm − 1 to 4000 cm − 1 (833-2500 nm). The spectral wavenumber interval was 3.85 cm − 1 with zero-filling equal to 2. The spectral resolution was 8 cm − 1 and 32 internal scans were averaged at each spectrum. The background was measured once per hour on reference Spectralon resin. Three FT-NIR spectra were collected at random locations on the weathered surface of each experimental sample. These three spectra were averaged, while the three samples of wood corresponding to replicated experimental conditions, i.e. same material and same treatment, were kept distinct in the data set assembly.
The interpretation of the spectra and derived loadings for the PAR-AFAC model, reported in Section 3 Results and Discussion, relies on literature references [12,16] and is summarized in Table 1.

Data fusion and pre-processing
The analysed data set includes all the registered data at each time point, namely: the twelve colour descriptors derived from digital images; the CIE Lab colour parameters; the NIR spectrum and the gloss parameters. Each data block was arranged in a distinct three-way array whose modes (dimensions) were (Fig. 2): • Mode 1: wood material type (7 wood types × 3 replica). This is common to all data blocks. • Mode 2: descriptors (X RGB/HSI : 12 pixels' quantifiers, X CIE Lab : 6 colour coordinates, X gloss : 2 gloss indices, X NIR : 2540 spectral points). • Mode 3: weathering time (reference + 4 exposure periods). This is common to all data blocks.
The three arrays, namely X RGB/HSI , X CIE Lab and X gloss , were concatenated at low-level data fusion and scaled to unit variance within the 2nd Mode to form X B2 [17].
The NIR spectra, X NIR , were pre-processed by second derivative, including smoothing (Savitzky Golay second order polynomial filter, 21 smoothing points) and formed array X B1 (Fig. 2).
Both data blocks (X B1 and X B2 ) were then fused at low-level. Scaling to equal block variance across the 2nd Mode was applied in order to have a balanced contribution in terms of variance. Finally, the array was centred across the 1st Mode prior to exploratory analysis.

PARAFAC decomposition
PARAllel FACtor (PARAFAC) analysis is a decomposition method for multiway arrays, based on the idea of parallel proportional profiles [18]. It assumes that a set of common factors can be used to describe simultaneously the variation occurring in several matrices with different weighting coefficients for each matrix. In this research, PARAFAC is used to highlight the weathering trends and kinetics (3rd Mode) that are characteristic of the various studied materials (1st Mode), as well as those which are the descriptors (2nd Mode) that are better capable of capturing the variation in the surface aspects and the chemical modifications that occurred.
In the case of a three-way array, the PARAFAC decomposition can be expressed by eq. (1): where: x ijk are elements of the array X (I × J × K) and e ijk are elements of the residuals array E containing the un-modelled part of the data array. The matrices A, B and C, whose elements are a if , b jf and c kf , respectively, hold the loadings for each of the three modes. The number of extracted components (factors) F is the same in each dimension. In analogy to Tucker3 model [18], g fff are the elements of the core array G, which is super-diagonal in PARAFAC. Core consistency and explained variance analysis were considered to determine the optimal number of factors F.

Post processing
It is possible to depict the weathering mechanism similarities and/or differences among the studied materials by inspecting the loadings scatter plot Mode 3, one factor vs. another, while exploration of Mode 2 loadings allows identifying the most salient descriptors. However, to have a clearer characterization of the weathering trend for each studied material the nested biplot representation has been used here [19]. In particular, the nested biplot representation allows recovering for each material its own weathering trend by combining in a single plot the loadings of Mode1 (samples) and Mode 3 (time) as detailed below.
The nested biplot consists of an alternative representation of the loadings for the different modes that is based on the particular rearrangement of the PARAFAC decomposition as expressed in eq. (2): where: s f , with elements s (ik)f , is a vector in which the rows consist of the fully crossed levels of two of the modes i, j or k for a given factor f. The new approach allowed examination of the behaviour for samples (Mode 1) at different weathering conditions (Mode 3) with respect to the measured variables (Mode 2). As Mode 3 corresponds to the time series, trajectories of the nested-mode biplot in the variable space are in fact a representation of how the aesthetical/functional aspects of assessed surfaces altered over time. The nested biplots were further elaborated to visualize the weathering trajectories and related deterioration kinetics. This was achieved by taking the sum of the loadings' values within the nested biplot over all the factors in the PARAFAC model. Subsequently, the value corresponding to no-weathered sample at the exposure time zero was removed from each trajectory. Finally, the plot of the absolute values versus time shows the weathering trajectories (as shown in the Results and Discussion section on Fig. 6).
The congruence loadings were computed to assess the relevance of the different surface state descriptors forming Mode 2 of the PARAFAC model. These congruence (or correlation) loadings reveal the relationship between the original variables and the components/factors resulting from the PARAFAC model [20], i.e. they represent the modelled variance by each variable. Congruence loadings were computed by CONLOAD function [21]. Congruence loadings are an easy and effective way for identification of the different tested sensors' role, and, thus, establish which are less relevant. This allows decision to exclude these sensors from the measurement portfolio to simplify monitoring in future routine analysis.

Results and Discussion
The fused data that describe natural weathering of modified woods were analysed by a three factor PARAFAC model. Fig. 3 presents the resulting third-mode loadings plot that expresses time-related kinetics of degradation. For clarity, the loadings values of the three replicas have been averaged in Fig. 3. The first factor shows a gradual decrease of its value from month three to month 12.
The analysis of Mode 1 loadings shown in Fig. 4a reveals that factor 1 distinguishes material #7 (high negative values) from all the others, laying relatively close to the origin of factor 1 on its positive side. In fact, weathering trajectory of acetylated wood (material #7) follows the trend captured by factor 1. This is due to the fact that the appearance of that wood surprisingly lightens in time, as visible in Fig. 1. That behaviour is related to chemical reactions of acetylated wood polymers induced by photodecomposition by sunlight and hydrolysis [22]. The trend of factors 2 and 3 with exposure time (Fig. 3, Mode 3 loadings) provides a highlight on the initial phase of the weathering changes; in particular, it indicates amplified effect of the changes in the six-and three-month periods, respectively. Such a trend can be noticed (Fig. 4b) for materials #4 and #7, which lightened sharply at month three (high negative loadings values on factor 2), and material #5, whose aspect changed till month six and then remained almost stable. On the contrary, material #2 seems to be the most stable as loadings values for all factors ( Fig. 4a and b) are close to zero. Material #3, #6 and #1 changed most at month three, showing significant loadings on factor 3.
These trends are even more evident on the nested biplot in Mode 1/3 (material type/exposure time) as presented in Fig. 5. The weathering trajectories corresponding to each modified wood material are presented regarding three scenarios: continuous deterioration (factor 1), moderate kinetic changes reaching stability at month six (factor 2) and rapid changes at the initial period of exposure (factor 3). The horizontal axis corresponds to exposure time, where vertical axis indicates relative change of the material state. Values recorded for three replica samples were averaged to better highlight the trend. In accordance with above observations, acetylated wood (#7) behaves differently than others, especially according to the trend revealed by factor 1. An additional kinetics component of the moderate velocity changes induced in the initial six months exposure of sample #7 is also recorded in factor 2, even if it is not much pronounced. Nearly no trace of the fast changes to acetylated wood (#7) are noticed in factor 3.
The extreme trends of the moderate kinetic changes (factor 2) are noticed for samples #4 and #5, even if both are in opposite directions. The fast changes due to weathering embodied in factor 3 are noticed for samples #1 and #3 and to minor extent for #6, again considering opposite signs for both cases as previously observed in Fig. 4b.
All the above observations of weathering trajectories can be    Fig. 6. A unique advantage for such visualization of the PARAFAC model output is the possibility of straightforward comparison of the degradation kinetics and deterioration trends simultaneously for all tested materials. The highest extent of changes is evident in sample #3 and #1 that corresponds to the reference wood (without any modification nor protection) and to the thermally modified wood impregnated with oil. Materials #4 and #5 present a similar trend of deterioration kinetics, even if both were finished with different techniques (silane and protective coating, respectively). The deterioration for a majority of samples seems to cease after three (#6) or six months (#4, #5). Some samples (#1 and #3) continued to deteriorate after initial intensive rise, even if the kinetics seem to be highly reduced after three months of weathering. The overall trend noticed for sample #7 (acetylated wood) is unique as its trend in the cumulative nested biplot is approaching the original stage after initial intensive rise. In contrast, thermally modified wood (#2) changed slowly but steadily along the whole 12 months of exposure.
It is important to notice that the multisensory approach used here disables direct link of the PARAFAC model loadings with the impressions of the surface aesthetical changes as presented in Fig. 1. It is due to the fact that not only a colour/gloss variation affects the model output, but also non-visible results of the chemical/physical changes as recorded in the near infrared spectrum. The role of different data blocks in the model and salience of each descriptor can be assessed by looking at the 2nd Mode loadings plots. In particular, for ease of interpretation, the loadings corresponding to the NIR spectrum (X B1 ) are presented in a separate plot (Fig. 7). The loadings corresponding to the other Fig. 6. Cumulative nested biplot for Mode 1/3 loadings of PARAFAC model for seven modified wood samples exposed to natural weathering. descriptors (X B2 ) are shown in Fig. 8. All the modelled variables that have congruence loadings higher than 0.6 are highlighted as a thick black line.
Detailed analysis of the Mode 2 loadings should be aligned with the PARAFAC Mode 1 loadings plot presented in Fig. 4. The process of wood acetylation results in decreased number of hydroxyls, followed by an increase of acetyl functional groups. This is recorded in Mode 2 loading (factor 1) as a series of peaks (8, 14, 15, 16, 25 and 26; see Table 1 for corresponding wavenumbers) assigned to CH 3 in acetyl ester groups of hemicelluloses [23]. Likewise, typical OH groups noticed for unmodified (chemically and/or thermally) woods are not present in either factor 1 or factor 2 loadings. The highest peak (6) corresponds to hydroxyl groups for all wood polymers that are alternated due to weathering [24]. This is especially related to the continuous decrease of lignin presence along the exposure duration that is depolymerized and leached out of the sample surface. According to Evans et al., losses of acetyl groups during weathering are a result of loss of acetylated lignin degradation products [25].
Factor 1 covering the NIR spectra seems to explain differences between chemical composition of samples, particularly of acetylated wood #7 as well as progress of chemical changes induced by natural weathering. In contrast, factor 2 loadings include several spectral features that are not directly interpretable by the band assignments table built on the basis of state-of-the-art know-how in wood NIR spectroscopy (Table 1).
It covers, therefore, an effect of diverse modifications and/or coatings implemented to study samples, particularly thermally modified and wood silicate impregnated (#4) or surface coated (#5) as well as acetylated woods (#7). Factor 3 discriminates samples #1, #2 from #3, #6 that correspond to different degree of wood modification extent. This is consistent with the high contribution of peak (4) assigned to holocellulose, which is known to be most affected by the thermal modification process. NIR spectra features of furfurylated wood are also most pronounced in factor 3 where, distinctive for this treatment, peak (19) shows high loadings.
Unfortunately, it is not a simple task to identify kinetics of the wood deterioration on the basis of Mode 2 loadings. However, several commonly recognized patterns of wood chemical changes due to natural weathering are present. These include spectral bands assigned to hygroscopic properties and hydroxyl groups in general (peak 10, 24), lignin (14,17), cellulose (5,14,20,24) and hemicellulose (14). All were identified as relevant in the preceding research of the authors [26].
Complementary information regarding the natural weathering of diverse woods is interpretable from Mode 2 loadings related to aesthetical properties of studied samples. Each aspect is recorded in the model to a certain degree, even if not particularly relevant properties are identified in factor 1 loadings as no variable possessed congruence exceeding 0.6. On the contrary, diverse colour or gloss aspects are important according to factor 2 and 3 loadings. Surface gloss, especially, contributes to factor 2, while diverse colour indicators are highlighted in both factor 2 and 3 loadings. Again, it is rather difficult to directly link these to the specific aspects of samples or to the kinetics of the deterioration. However, it is evident that additional (to NIR spectra) information greatly contributes to the model reliability and considerably increases the multi-sensor method's reliability.

Conclusions
Understanding the mechanism of changes to biomaterials due to natural weathering is an important aspect necessary for increasing the confidence of using such materials in construction. Multi-sensor evaluation is highly recommended to assure better reliability and generality of characterization. A challenge, however, is the proper integration of diverse data collected by different sensors. This paper presents an alternative approach for multi-sensor data fusion and modelling of the deterioration processes by means of PARAFAC model. The original protocol for data fusion is proposed as well as novel meta parameters, such as cumulative nested biplot. The pilot research reported here confirms suitability of this approach for modelling the natural weathering processes for seven types of wood samples that were treated and coated according to state-of-the-art industrial procedures. The data were presented in the three dimensioned block form and, consequently, analysed with an innovative chemometric approach.
It was possible to successfully differentiate diverse sample types on the basis of NIR spectra and selected surface appearance indictors. The characteristic trends for each material deterioration were identified. It may serve as an important design guidance for architects as well as a useful tool for engineers to predict aesthetical or functional changes of diverse materials during the service life of wooden objects exposed to natural weathering.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.