3.2.4. Step 4—Radiative Heat Flux (Qrad)

The estimation of radiative heat flux (*Qrad*) from an area of IR frame mainly characterized by thermal anomaly (Region of Anomaly, *RoA*) is a newly proposed processing technique that can offer an interesting contribution to the investigation of possible variations of radiative thermal emissions.

In order to estimate *Qrad*, which is the thermal energy emitted per unity of area in a unity of time, the *RoA* has to include pixels whose temperatures are representative of the main thermal anomaly. Nevertheless, the *RoA* is usually not homogeneous and it is characterized by the presence of both high temperature sources (fumaroles) and low temperature sources (surrounding emission-free rocks). In addition, when sensor-target distance is more than approximately 10 m, the pixels of *RoA* can be several centimeters large, and therefore, some temperatures are underestimated if their pixels integrate both high and low temperatures [42,43]. Consequently, the variations trend of *Qrad* can be sensibly flattened. A solution to this problem is to calculate the Standard Deviation (SD) of pixels' temperatures of a specific *RoA* and then to use only temperature values (*TROAH*) greater than 2SD to estimate *Qrad*.

Finally, the *Qrad* of a specific *RoA* (W/m2) is calculated by using the Stefan-Boltzmann equation:

$$Q\_{rad}RoA = A \sum\_{i=1}^{n} \sigma \varepsilon (T\_{RoA}H\_i)^4 \tag{6}$$

where σ is the Stefan-Boltzmann constant, ε is the emissivity (for pyroclastic rocks is assumed to be 0.9) and *A* is the investigated area size (m2) obtained by multiplying pixel area and length n of *TROAH* time-series.

The detection of any possible change of *Qrad* trends, even though related to a specific *RoA*, allows to better characterize thermal behavior of the studied area if *RoA* is representative of the main thermal anomaly.

The use of de-seasoned time-series of temperature values (TdesTS) is essential in order to evaluate *Qrad* changes due to endogenous sources only. This means that only the dataset processed with STLd can be used.

The Matlab© code performing *Qrad* ('step04.m') is available as Supplementary Materials and its functionalities are illustrated in the Appendix A.

### 3.2.5. Step 5—Yearly Rate of Temperatures Change (YRTC)

The thermal variations, in a defined time interval, of every single pixels of IR frame, can be evidenced by evaluating the yearly rate of temperatures change (YRTC). This kind of elaboration produces a map of the IR frame, according to a color scale, of yearly rate of change of pixels' temperature. The yearly rate of temperatures change is represented by the values of slope coefficients of the linear fit of time-series temperatures of every pixel. On the other hand, the selected time interval has to be characterized by a progressive increase or decrease of maximum temperatures of IR frame, according to a correspondence as linear as possible. This needs a preliminary investigation of the temperatures trend over time.

The YRTC map is created by overlapping the values of slope coefficients on a picture (in the visible range) of the framed area. In order to show the yearly rate of change values of pixels whose temperature time-series best fit a linear model, a mask was applied. This mask allowed the display of values related to pixels whose linear regressions of temperature time-series had coefficients of determination (R2) higher than a user-defined threshold value.

The YRTC map gives the opportunity to evidence possible connections between temperature increase/decrease and geological features of the monitored site (Figure 7).

**Figure 7.** Yearly rate of temperature change maps of SF1 in the period 2016.02.01–2016.11.30. Map (**a**) has R<sup>2</sup> threshold value = 0.2. Map (**b**) has and R2 threshold value = 0.45. In this time-interval SF1 maximum temperatures decreased of about 10 ◦C as evidenced by the color map.

The Matlab© code performing YRTC data ('step06.m') is available as Supplementary Materials and its functionalities are illustrated in the Appendix A.

#### *3.3. System Automation and Graphic Interface*

The above-described methodologies were performed as steps by Matlab© functions, which can be executed with a command line or managed by a user-friendly graphic interface (GUI). Settings can be saved in user-defined configuration files. Due to the modular structure of the processing steps, they can be performed singularly or grouped in an automated sequence in order to execute the whole procedure at defined time by using the GUI that integrates the automation code. Automation is necessary if IR data processing is aimed to surveillance purposes.

The GUI Matlab code (asira\_gui.m) is available as Supplementary Materials and its functionalities are illustrated in the Appendix A.

#### **4. Results and Discussion**

In order to discuss the advantages and the limits of the above presented processing methodology, the results obtained by applying the five processing steps are reported. Two datasets were processed: (1) the first consisted of 2.901 IR JPEG frames acquired in the period 2016.01.27–2019.01.13 at Solfatara 1 (SF1) station; (2) the second consisted of 5.850 IR JPEG frames acquired in the period 2013.03.26–2019.01.13 at Pisciarelli (PS1) station.

#### *4.1. Data Quality Selection*

The relation (1), discussed in §3.2.1, was used to remove low-quality IR frames before starting the analysis of data. The efficiency of this procedure depended on the choice of the coefficient *c* which was influenced by the statistical distribution of data. Low values of Standard Deviation of IR frames temperatures were an indicator of low quality data and the lower the coefficient *c*, the higher the number of IR frames discarded as low quality ones. The analysis of data acquired by SF1 suggested c=1 as an appropriate value (Figure 8), as the visual inspection of discarded frames (about 11% of total frames) confirmed that they were mainly low-quality ones. This kind of preliminary analysis had to be made to every dataset from different stations as the coefficient *c* can be different depending on the physical and geometrical characteristics of framed area and IR sensor.

**Figure 8.** Frequency distribution of Standard Deviation values of IR frames temperatures (σF) acquired at SF1 station. The line 'Lower Threshold', which is defined by the relation (1) with *c* parameter equal to 1, splits good-quality frames (on the right of the line) and low-quality frames (on the left of the line).

#### *4.2. Seasonal Component Removal*

Two different methodologies of seasonal component removal are used in order to process IR datasets having different time-length. The background removal procedure (BKGr), previously proposed to seasonal correction [19,20], is suitable to very short datasets even though it has some limitations in the final output. The main limit was that the removal of seasonal component produces only residuals of maximum or median values of temperatures instead of absolute temperature values. Although this kind of analysis does not take full advantage of all the intrinsic information contained inside the IR frames, the BKGr method generated trends of temperature residuals which provide adequate information to characterize the thermal behavior of studied area.

The STL decomposition method (STLd), proposed for the first time in this work to remove seasonality to IR temperature time-series, needed a nearly two-year long datasets due to the statistical approach of the robust and widely applied algorithm. By using STLd method, it was possible to estimate the Seasonality as a separate component of temperature time-series. This feature allows the removal of seasonality to all pixels of IR frames, giving the opportunity to apply further analysis methods (e.g., radiative Heat Flux estimate), which needed the whole frame to be de-seasoned. In Figure 9, the background area boundaries (Figure 9a) inside the SF1 IR frame and the plot of Trend component of background area, obtained by applying STLd method are reported (Figure 9b). The constant and flat temperature values of background Trend confirmed the appropriate choice of this area.

**Figure 9.** Background area boundaries inside the SF1 IR frame (**a**) and plot of Trend component of background area which is obtained by applying STLd method (**b**). Trend values varie between 14.35 and 14.7 ◦C.

The Trend component evaluated by the STLd method was useful to estimate the long-term thermal behavior of the studied area even though it was not suitable for short-term observations. In order to describe short-term thermal behavior, aimed to surveillance purpose, it was necessary to merge both Trend and Reminder components to obtainT+R plots (Figure 10, blue line plot).

**Figure 10.** Plots of temperatures acquired at SF1 station removed of seasonal component. Red line = temperature residuals obtained by applying BKGr method; blue line = Trend+Reminder values obtained by applying STLd method.

Despite the reported limits of BKGr method, the comparison between temperature residuals plots with the BKGr method and T+R plots by STLd method of SF1 IR frames (Figure 10) showed a close similarity of data trends. This similarity confirms the effectiveness of the BKGr method to process datasets shorter than two years.

#### *4.3. Radiative Heat Flux Estimate*

The computation of radiative heat flux was available on IR frames where seasonality was removed by applying the STLd method. In order to obtain the correct trend of radiative heat flux of a definite area, a correct selection of area boundaries was necessary. The heat flux computation strongly depended on the Hpix BKGpix ratio of the selected area, where Hpix is the number of pixels related to thermal anomaly and BKGpix is the number of pixels related to emission-free rocks. The higher this ratio is, the more accurate the radiative heat flux estimation. This way, the choice of boundaries of processed areas had to be made in order to include as many Hpix as possible. The solution to attenuate the underestimate the heat flux due to the presence of BKGpix, proposed in the §3.2.4, was to select pixels whose temperatures were greater than 2σ of the frequency distribution of temperatures from selected area. The plots reported in Figure 11 show how efficient this kind of solution was. In this figure, plot a) reports heat flux time-series of Areas 1, 2 and 3 evaluated by selecting all the pixels inside each area; plot b) reports heat flux time series of the same areas, evaluated by applying the selection of pixels greater than 2σ of temperatures frequency distribution. The blue line plot is from Area 1, which only includes the major thermal anomaly of the SF1 frame, characterized by the higher Hpix BKGpix ratio. Red line plots of Area 2 and black line plots of Area 3 are representative of lower Hpix BKGpix ratios due to higher number of BKGpix included in the selected areas.

**Figure 11.** Heat flux plots of selected areas inside SF1 frames. Area 1 (c) includes the mayor thermal anomaly, Area 2 and 3 (c) include emission-free rocks. Plot a) = heat flux trends (smoothed with window = 29) of Areas 1, 2 and 3 evaluated by selecting all the pixels inside each area. Plot b) = heat flux trends (smoothed with window = 29) of Areas 1, 2 and 3 evaluated by applying the selection of pixels greater than 2 s of temperatures frequency distribution.

The comparison between Figure 11a,b evidences an underestimate of heat flux values when the computation includes all the pixels inside the selected areas (Figure 11a). Figure 11b was obtained selecting only the pixels whose temperatures values were greater than 2σ and showed a remarkable decrease of heat flux underestimate, better evidencing trend variations.

#### *4.4. Yearly Rate of Temperature Change Estimate*

As reported in §3.2.5, the final product of this processing step was a color scale map of the yearly rate of temperature change (YRTC) values overlapped to a picture (in the visible range) of the framed area. YRTC data were filtered according to a threshold value of the coefficient of determination (R2) of the linear regressions of pixels' temperature time-series. Two different examples of yearly temperature rate of change maps of PS1 area in the same time-interval (2016.03.10–2016.07.10) are reported in Figure 12. In this time-interval the PS1 temperatures were subjected to an increase of about 10 ◦C. The maps of Figure 12 only differ in the choice of R2 threshold value; hence, a correct choice of these parameter is critical to produce a map that is easy to comprehend. Map **b** (R2 = 0.7) shows better evidence of pixels whose temperatures rate of change time-series values best fit a linear model than map **a** (R<sup>2</sup> = 0.5). The ASIRA code allows the user to select both color scale limits and different values of R<sup>2</sup> threshold by using a user-friendly GUI in order to achieve the right balance between optimal visual result and reliability of data visualized.

**Figure 12.** Yearly rate of temperature change maps of PS1 in the period 2016.03.10–2016.07.10. Map (**a**) has R2 threshold value = 0.5. Map (**b**) has and R<sup>2</sup> threshold value = 0.7.

#### **5. Conclusions**

Relevant contribution to the surveillance of volcanic areas affected by thermal anomalies can be provided by monitoring the spatio-temporal evolution of surface temperatures field. The acquisition of IR image data by ground-based monitoring network is an effective tool to perform this task. However, the analysis of IR data time-series is not easy to accomplish due to the influence over IR temperatures of both exogenous and endogenous processes.

In this paper, we have presented a unique operational processing chain developed in Matlab© environment which allows the detection and quantification of possible changes in time and space of the ground-surface thermal features. This application (ASIRA, Automated System of InfraRed Analysis) performed a multi-step procedure that generated both trends of temperatures and heat fluxes as well as maps of yearly rate of temperatures change. The procedure implemented new algorithms based on improvements of previously proposed methods and also original techniques aimed to effectively remove seasonal component of IR temperature time-series and to evaluate radiative heat fluxes of thermal anomaly areas.

ASIRA can be performed as separate steps or executed in a fully-automated way by using a user-friendly graphic interface. The Matlab© code of ASIRA and the Operative Manual are included as Supplementary Materials.

The ASIRA code was applied to process IR data acquired by stations of TIRNet surveillance network operated by the Osservatorio Vesuviano, section of National Institute of Geophysics and Volcanology (INGV) at Campi Flegrei volcanic area (Italy). The results show the effectiveness of this method to provide a valuable contribution to the continuous monitoring of thermal anomalies related to studied areas.

This operative tool has been conceived for volcanic surveillance of diffuse degassing areas and low-temperature fumarole fields which variations may precede significant phases of volcanic unrest. Notwithstanding, the procedure can be applied to monitor different volcanic scenarios (i.e., lava-flows, active volcanic vents and eruptive fractures) but also different natural and environmental hazards (fires, waste-disposal sites, pollution discharges, landslides, etc.).

**Supplementary Materials:** the Matlab© code of A.S.I.R.A. (Automated System of InfraRed Analysis) which is described in Appendix A), and the Operative Manual (pdf file) are provided at the following link: http: //www.mdpi.com/2072-4292/11/5/553/s1.

**Author Contributions:** F.S. was responsible for the ASIRA design and implementation in the Matlab© software, data acquisition and analysis, and writing the manuscript. G.V. was responsible for the research design, data analysis and writing the manuscript. Both authors contributed to the software validation.

**Funding:** TIRNet monitoring network was partially funded by the 2000–2006 National Operating Program (NOP) and by SISTEMA project, which has been developed in the framework of the Campania Regional Operating Program (ROP) FESR 2007-2013.

**Conflicts of Interest:** The authors declare no conflict of interest.
