Using Multitemporal Sentinel-1 C-band Backscatter to Monitor Phenology and Classify Deciduous and Coniferous Forests in Northern Switzerland

Efficient methods to monitor forested areas help us to better understand their processes. To date, only a few studies have assessed the usability of multitemporal synthetic aperture radar (SAR) datasets in this context. Here we present an analysis of an unprecedented set of C-band observations of mixed temperate forests. We demonstrate the potential of using multitemporal C-band VV and VH polarisation data for monitoring phenology and classifying forests in northern Switzerland. Each SAR acquisition was first radiometrically terrain corrected using digital elevation model-based image simulations of the local illuminated area. The flattened backscatter values and the local area values were input to a temporal compositing process integrating backscatter values from ascending and descending tracks. The process used local resolution weighting of each input, producing composite backscatter values that strongly mitigated terrain-induced distortions. Several descriptors were calculated to show the seasonal variation of European beech (Fagus sylvatica), oak (Quercus robur, Quercus petraea) and Norway spruce (Picea abies) in C-band data. Using their distinct seasonal signatures, the timing of leaf emergence and leaf fall of the deciduous species were estimated and compared to available ground observations. Furthermore, classifications for the forest types ‘deciduous’ and ‘coniferous’ and the investigated species were implemented using random forest classifiers. The deciduous species backscatter was about 1 dB higher than spruce throughout the year in both polarisations. The forest types showed opposing seasonal backscatter behaviours. At VH, deciduous species showed higher backscatter in winter than in summer, whereas spruce showed higher backscatter in summer than in winter. In VV, this pattern was similar for spruce, while no distinct seasonal behaviour was apparent for the deciduous species. The time differences between the estimations and the ground observations of the phenological events were approximately within the error margin (±12 days) of the temporal resolution. The classification performances were promising, with higher accuracies achieved for the forest types (OA of 86% and κ = 0.73) than for individual species (OA of 72% and κ = 0.58). These results show that multitemporal C-band backscatter data have significant potential to supplement optical remote sensing data for ecological studies and mapping of mixed temperate forests.


Introduction
Forests provide valuable ecosystem goods and services for human well-being [1].Therefore, it is of great interest to gather information about the state of the forest ecosystems.Because forest inventories are costly, it is essential to develop cost-effective mapping methods to allow management of forests [2].In times of climate change, periodic assessments are in demand, as changes in phenology [3,4] as well as shifts in species composition [5,6] become apparent.It is important to discover efficient ways to continually monitor forested areas and subsequently gain knowledge about processes and changes.
Remote sensing offers valuable tools to implement mapping and monitoring approaches [7].Until now, mainly optical sensors have been used to monitor forests on a global scale, ranging from intra-annual studies [8,9] to others examining several years [4,10].Using optical sensors however, two main disadvantages are apparent.First, optical sensors are not sensitive to vegetation structure.Second, temporally consistent data acquisition is challenging, as the optical sensor's views are often impaired by cloud cover.To overcome these disadvantages, it is important to consider synthetic aperture radar (SAR).SAR is sensitive to vegetation structure [11] and is acquired independent of daylight and weather, so it can more easily produce temporally consistent data [12].
Until now, SAR data have often been used to estimate biomass [13,14] and map the extent of deforestation [15] in tropical forests.For boreal forests, several methods working with SAR data have been applied, e.g., for mapping forest types [16,17] or estimating growing stock volume [18].Because of frequent cloud coverage, SAR has to be considered in the tropics.In the northern boreal forests, however, SAR is also relevant, as long periods with sparse or no daylight pose a further barrier to using optical remote sensing data.It is useful to assess the potential of SAR data for mixed temperate forests, as only a few studies assessing this potential have been conducted until now [19,20].Given that significant structural changes due to phenology happen within a short time period in spring and autumn, reliable observations are critical.
The two recently launched Sentinel-1 (S-1) satellites, carrying C-band SAR sensors, offer new possibilities in the compilation of time series.Generally, they acquire backscatter in vertical-vertical (VV) and vertical-horizontal (VH) polarisation mode over landmasses [21].S-1 data have already been successfully applied to several land applications.Both VV and VH polarisations have been used to monitor crop growth of several cultures in Angola [22] and of rice in Myanmar [23].Balzter et al. [24] and Abdikan et al. [25] showed that using VV and VH polarisation data leads to satisfying classification of land cover classes, applying random forest (RF) and support vector machine (SVM) classifiers, respectively.Other studies exhibited the potential of S-1 data for burnt area mapping in forests [26], monitoring wetlands [27] or snow wetness estimation in the Alps [28,29].S-1 offers superior spatial and temporal coverage and much improved radiometric calibration stability in comparison to former spaceborne C-band sensors [30].
Previous work has shown that SAR forest backscatter is very complex, being dependent on several sensor and object properties.Generally, sensors with high frequencies such as X-or C-band are less able to penetrate the crown than are lower frequencies such as S-or L-band [31].Thus, the dominance of the crown return is stronger at high frequencies.The polarisation used by the SAR sensor also impacts the backscatter.Cross-polarised backscatter from forested areas leads to higher correlation with biomass than co-polarised backscatter [32].Object properties impacting forest backscatter can be divided into two groups.The first group consists of properties that influence the dielectric constant of the trees, such as the moisture content of branches and leaves [19], external moisture conditions [33], and the temperature of the wood [34,35].The second group consists of structural properties such as the size and orientation of branches and leaves [36], the spatial pattern of trees [37], and whether or not foliage is present [38,39].The last is especially interesting, as it could allow the monitoring of phenology at a land surface scale.
Until now, most studies with the aim of mixed temperate forest characterisation have focused on extracting additional information by using fully polarimetric SAR data to discriminate different backscatter behaviours of forest types [40,41] or estimate forest density using Wishart supervised and SVM classifiers [42].Collecting multiple SAR acquisitions into longer time series has the potential to reveal temporal signatures, as shown by Sharma et al. [33] with C-band horizontal-horizontal (HH) polarisation data for forest and Gamba et al. [43] with C-band VV polarisation data for several land cover types.However, only few studies have analysed longer C-band time series of forested areas.Ahern et al. [39] and Proisy et al. [19] found phenology effects on C-band SAR backscatter for coniferous and deciduous forest types using HH and VV polarisations, respectively.A recent study using S-1 data showed the impact of phenology of deciduous forests to be much higher in cross-polarised than in co-polarised data [20].
In this paper, we present an analysis of Sentinel-1 SAR radiometrically calibrated backscatter time series of an unprecedented set of VV and VH polarisation observations over mixed temperate forests in northern Switzerland.One coniferous and two broadleaved, deciduous species were investigated.We show how to build the C-band time series and to subsequently extract descriptors for the different forests investigated.A detailed analysis of these descriptors is provided, focusing on differences between tree species and between the investigated years.Finally, we show how our findings help to classify mixed temperate forests and how they could help in future monitoring of the seasonal vegetation cycles of forests.

Study Areas
The areas studied were the whole of the Canton of Zurich (ZH), and three study areas in the Canton of Bern (BE) in Switzerland.For the main analysis, two areas in the north of ZH and the three areas in BE were specifically selected based on their high species homogeneity.This resulted in five main study areas with forest stands of variable extents across northern Switzerland being evaluated (Table 1 and Figure 1).To be considered for our study, the five study areas had to consist of mixed temperate forest stands with a high homogeneity (≥80%) of single species.The three species examined were European beech (Fagus sylvatica), oak (Quercus robur, Quercus petraea) and Norway spruce (Picea abies).In addition, the extent of the study areas had to be large enough to allow an erosion of edge pixels (see Section 2.3.2 for details) to limit border effects.The elevations of the five study areas range between 390 and 670 m a.s.l.For some analyses, the entire forested area within the Canton of Zurich, depicted in the green rectangle, was used.The background aerial images were provided by the Swiss Federal Office of Topography swisstopo [44].Reproduced by permission of swisstopo (BA17116).Swiss map coordinates: CH1903+/LV95.

Data
We used Sentinel-1 C-band Interferometric Wide swath (IW) mode Ground Range Detected High Resolution (GRDH) images covering the study areas.The chosen time span was between 1 January 2015 and 24 May 2017.Because both available polarisations were considered, 433 VV-and 433 VH-images were used as input to the investigation [30].The IW mode specifications are listed in Table 2. To process the SAR data, a digital terrain model (DTM) was available with a spatial resolution of 2 m [45].For the purpose of this study, this DTM was resampled to a spatial resolution of 10 m to match the resolution of the S-1 images.
In-depth information about the forests in the Canton of Zurich were provided by the aerial image forest stand map released by the Office of Landscape, Agriculture and Environment (Amt für Landschaft und Natur) of the Canton of Zurich in 2010 [46].The vector dataset contains information derived from aerial imagery interpretation about several characteristics of the forest stands such as their outline, relative tree species abundance, tree cover density and stand height.For this study, the dataset was converted into a raster layer with a spatial resolution of 10 m using the Geospatial Data Abstraction Library (GDAL) tool rasterize [47].
Meteorological and phenological data of the study areas were also prepared.MeteoSwiss provided daily mean temperature and precipitation raster layers with a spatial resolution of 1.25', corresponding approximately to 2.3 km in northing and 1.6 km in easting [48].Both layers were produced by spatially interpolating data of the MeteoSwiss weather station network of Switzerland [49].Annual phenological observations were available, released by MeteoSwiss for selected weather station locations for the species beech and spruce [48].For beech, the timing of leaf emergence in spring, and leaf colouring and leaf fall in autumn is recorded by volunteers, for spruce, only needle shoot in spring.Observations of the two stations Rafz and Biel were selected for this study based on their proximity (<25 km) to the investigated forest stands of the study areas 1, 3 and 4 (see Table 1).For the study areas 2 and 5, no suitable reference emergence and fall observations were available.

Methods
In a first step, the SAR data were radiometrically calibrated and geocoded, and time series were built (Figure 2 and Table 3).Subsequently, masks of the outlines of the forest stands were generated to analyse the backscatter behaviour of the different species and forest types.Then, several descriptors were calculated to assess the temporal backscatter patterns within the specific forest stands (Figure 3).In a last step, classifications were implemented by making use of the significant differences in the temporal backscatter signatures of the different species (Figure 4).The acquired Sentinel-1 IW VV and VH acquisitions were processed in two different ways, as illustrated in Figure 2. First, the acquisitions were radiometrically calibrated and terrain-geocoded with the available DTM, which resulted in geometrically terrain corrected (GTC) γ 0 E images for both polarisations [50].These GTC images had a pixel spacing of 10 m.Second, to reduce the strong influence of topography on the backscatter [54], the VV and VH acquisitions were also radiometrically flattened and terrain-geocoded according to the methodology described in [52].This approach corrects not only the geometry of the scene but also its radiometry.This is achieved by normalising backscatter based on the local illuminated area as seen by the sensor rather than a rough ellipsoid-model-based incidence angle as in GTC processing.The resulting radiometrically terrain corrected (RTC) γ 0 T images exhibit an easier interpretability, which is especially valuable in areas with complex terrain.The pixel spacing of these RTC images was also 10 m.
In a second stage, a multitemporal compositing approach was applied to the RTC images of both polarisations to enhance spatial resolution, and reduce noise and the influence of speckle [53].The local resolution weighting (LRW) approach combines RTC images of ascending and descending acquisitions within a specific timespan into γ 0 LRW composite backscatter images.Because the local image resolution in an RTC image is inversely proportional to the local illuminated area, the approach weighs all the images' contribution to each pixel according to the local resolution within each contributing backscatter image.A timespan of 24 days and a temporal sampling interval of 12 days were found to adequately observe seasonal backscatter behaviour while also ensuring full areal coverage, even when fewer images than typical were available due to calibration activities in southern Germany.The compositing led to reliable time series of SAR observations of the study areas.Only for one single 24 day period could no composite of the study areas in ZH be produced.Between 6 June 2015 and 29 June 2015 no acquisitions were available in northern ZH, due to a calibration campaign of the German aerospace centre (DLR) in southern Germany [55].The evolution of the number of acquisitions available from 2015-2017 is shown in Figure 5.There are three reasons why the number of available acquisitions changed considerably over the investigated time period for both study areas.The low number of available acquisitions in summer 2015 (2 acquisitions) and 2016 (3 acquisitions) were caused by the DLR campaign.A gradual increase to a maximum of 18 acquisitions towards the end of the investigated time period was due to a (a) ramp-up of the production by the European space agency (ESA) payload data ground segment and (b) the addition of Sentinel-1B acquisitions in September 2016 [56].A backscatter calibration issue was mitigated by applying corrections.Both Miranda [57], and El Hajj et al. [58] reported positive Sentinel-1A antenna gain biases between 19 March 2015 and 25 November 2015.Hence, all backscatter values within that timespan were corrected.After [58], a correction of −0.9 dB was applied within that time period in 2015.

Species Mask Generation
To analyse the backscatter behaviour of γ 0 E images and γ 0 LRW composite time series, masks for the three different tree species and the two forest types were produced following a scheme shown in Figure 3.For the study areas in the Canton of Zurich, species and forest type masks were generated by segmenting the aerial image forest stand map according to its delineation of the forest stands.For the species masks, this was done by choosing all pixels with high species homogeneity for beech, oak and spruce.For the forest type masks, the considered forest stands consisted of either deciduous or coniferous species.Because there were not enough pixels of absolutely homogenous stands, pixels with species homogeneity greater than 80% were chosen.The threshold of 80% was chosen, as a trade-off between the available number of pixels and ensuring a considerable dominance of the species within the forest stand.Given that no aerial image-based forest stand map was available for the study areas in BE, the masks for these areas were produced based on information from local foresters [59][60][61].This spatial information about the forest stands was subsequently rasterised to produce the three masks.
All the masks were "eroded" to minimise edge effects and ensure the presence of the desired species or forest type.This was done using an erode function with a square mask of 3 × 3 times the sample interval.Similar to the homogeneity threshold, the 3 × 3 size was chosen as a trade-off between maximising the available number of pixels and mitigating edge effects.

Comparison of an End-of-Winter and Summer GTC Image (γ 0
E ) and LRW Composite (γ 0 LRW ) In a first analysis, the backscatter of the forested area in the whole of ZH of two γ 0 E images was compared.The chosen images were from 6 June 2016 for summer and leaf-on, and 16 March 2017 for end of winter and leaf-off.The images were masked with the generated forest type masks and the median, 25 and 75% percentiles were calculated within the masked areas ( x, Q1, Q3).Subsequently, the medians were compared for the two dates and the two forest types.RGB-overlays, histograms and correlation density plots were prepared to highlight differences in backscatter distribution between the two acquisition dates and between the two SAR processing methods.RGB-overlays and histograms were produced using the above-mentioned γ 0 E images; others were generated using γ 0 LRW composites of the timespans 24 May-16 June 2016 and 2 March-25 March 2017.In the RGB-overlays, the summer image was assigned to the red and blue channels, and the end-of-winter image was assigned to the green channel (Figure 6).
The advantages of the RTC processing (γ 0 T ) and the LRW compositing approach (γ 0 LRW ) over the GTC-based processing (γ 0 E ) are apparent.The colour pattern in the γ 0 LRW RGB-overlay (b) appears much smoother, leading to a better discrimination between the two forest types.In addition to that, the strong influence of topography on the γ 0 E RGB-overlay (a) is seen (e.g., steep slope right next to the river bed in the northwest).In areas with steep backslopes (c), backscatter was lower.In the γ 0 LRW RGB-overlay on the other hand, topographic effects were mitigated.

Time Series Analysis of LRW Composites (γ 0
LRW ) In the main analysis, the temporal behaviour of the forest backscatter was investigated.The descriptors listed on the first row of Table 4 were calculated for every γ 0 LRW composite within the eroded masks ( x, Q1, Q3).Subsequently, these descriptors were plotted in the temporal domain and their seasonal behaviour was inspected visually.To analyse this behaviour further, the temporal descriptors listed in Table 4 on the second row were calculated for the previously established descriptor time series.The medians within the time series were derived from all the calculated medians and 25 and 75% percentiles of the masked composites ( xTo, Q1To, Q3To).The medians of winter ( γ0 win ) and summer ( γ0 sum ) were calculated, as visual inspection of the plots indicated a change in backscatter between the two seasons.The timespan for the two seasons was chosen conservatively to be 1 June-15 September for summer and 1 December-15 March for winter.This was done to ensure that no 'winter' composites with leaf-off conditions appeared in the 'summer' subset and vice versa.Their difference (δ) was calculated according to Equation (1).
The timing of the change in backscatter was also of interest.Natural breakpoints in the time series were calculated using the algorithm described in [62].The breakpoints algorithm allows detection of structural changes in time series [63] and has successfully been used in previous ecological studies (e.g., [64,65]).The time series were split into annual sets to facilitate detection of the break dates.The algorithm assigns break dates where the OLS-based CUSUM test [66] registers the two highest empirical fluctuation process values [63].The first break date within each time series was assumed to be the leaf emergence (BD1), while the second break date was assumed to be the leaf fall of the respective year (BD2).The missing value due to the absent γ 0 LRW composite in summer 2015 for the study areas in the Canton of Zurich (1 & 2 in Table 1) was linearly interpolated [67] to facilitate use of the breakpoints algorithm.
The last step consisted of a comparison of the calculated descriptors for the different forests.General differences between deciduous and coniferous forests were compared as well as differences between the two deciduous species beech and oak.The phenological ground observations, when available, were juxtaposed with the breakpoints-based break dates.

Classification of Species and Forest Type
We conducted two classifications of the forested areas in the Canton of Zurich using a random forest classifier (RF) [68].In addition to a classification of the three species examined, another was performed on the forest types deciduous and coniferous.RF was chosen, as a model was needed that could handle the non-linear relations in our data.Moreover, a study using multitemporal C-band SAR data showed that RF outperformed other classifiers [69].In addition to the actual classification, it also ranks the importance of the used predictors.We used the MATLAB class TreeBagger to perform the classification [70].After an initial model validation phase of testing different parameter combinations, we chose the following parameters for both approaches.We generated a forest of 100 decision trees with a minimum terminal node size of 1.All the other parameters were set to the default values of TreeBagger.
The classification processing chain is illustrated in Figure 4. First, twelve independent predictors were calculated per pixel within the γ 0 LRW composite time series.The predictors used were (1-8) the respective median of the years 2015 and 2016 as well as the respective δ between winter (leaf-off) and summer (leaf-on) medians at both polarisations.The definition for the seasonal periods was again chosen conservatively, with 1 December-15 March for winter and 1 May-15 September for summer.Two break dates (9)(10)(11)(12) for the two years were derived per pixel VH time series, using the breakpoints algorithm.These break dates completed the set of predictors.
Second, the training and validation datasets were prepared.This was done in the same way as for the time series analysis (described in Section 2.3.2) by generating forest masks, applying the erode function and subsequently randomly selecting pixels for the two datasets.The mask generation and erosion for both classifications approaches were performed in the same way.Then, an RF model was trained with a set of randomly selected pixels.Another set of randomly selected pixels was then used to validate the model.The number of selected pixels for each class was balanced for both training and validation sets to avoid biasing in the model.This procedure was identically executed for the forest type and the species classification.
To ensure a representative estimation of the performances, multiple RF models with different training and validation datasets were trained, and mean accuracies were calculated after each run.Then, the number of runs necessary to stabilise the achieved mean accuracies was evaluated.It showed that going on to calculate more than 40 models did not significantly change the mean accuracies (±0.25%).The number of RF models was set to 40 for both approaches.The performances of the two classifications and their utility were then evaluated by generating the mean confusion matrix and calculating mean overall accuracy, producer's accuracy, user's accuracy and Cohen's kappa of the 40 different RF models.The "predictor importance" measure was also assessed.For every predictor, that measure specifies the loss in prediction accuracy when the values of that predictor are distributed randomly across the pixels [70].

Comparison of the Seasonal GTC Images (γ 0
E ) and LRW Composites (γ 0 LRW ) Table 5 shows the results of the calculation of the median, 25 and 75% percentiles backscatter within the forest type masks in the two γ 0 E images for both polarisations.Generally deciduous species showed higher backscatter than coniferous stands.With the exception of VH polarisation in summer, the difference between the two forest types was quite substantial.A seasonal change in VH backscatter within both forest types was observed.Coniferous forests had stronger backscatter in summer (6 June 2016) and lower at the end of winter (16 March 2017), whereas deciduous forests showed the opposite, with higher backscatter at the end of winter and lower in summer.This seasonal change is also apparent in a VH γ 0 E RGB-overlay (Figure 6a) and a γ 0 LRW counterpart (Figure 6b) of the study area in Teufen.Areas in magenta showed higher backscatter in summer, whereas areas in green showed higher backscatter at the end of winter.It can therefore be expected that areas in magenta consisted mainly of coniferous species, and areas in green mainly of deciduous species.Compared with the forest type map of the area (Figure 6d), a striking correspondence is visible between the colour patterns of the RGB-overlays and the forest type classes.
The distribution of the backscatter values of the two forest types in the two γ 0 LRW composites is depicted in Figure 7.The seasonal change between end-of-winter and summer backscatter can be observed in these correlation density plots.In coniferous forest (a), the highest densities were reached above the equivalence line with higher backscatter in summer.In deciduous forest (b) the opposite was observed, with highest densities reached below the equivalence line with higher backscatter at the end of winter.The sign of these trends is consistent with observations made in a previous work [20].The seasonal change is also depicted per pixel in the difference between end-of-winter and summer VH backscatter of the two forest types in both SAR products (Figure 8).Two features are apparent.First, values for deciduous and coniferous forests were mainly positive and negative, respectively.Second, the variance within the γ 0 LRW composites (b) was considerably lower.This is expected to allow a more reliable separability between the two forest types.

Time Series Analysis of LRW Composites (γ 0
LRW ) The temporal evolution of the spatial descriptors is shown in Figure 9. Inspecting the γ 0 LRW composite time series, the same features as in Section 3.1 were observed.In addition to the spatial descriptors, temporal descriptors were calculated for the approximately 70 γ 0 LRW composites.Examining these plots, the seasonal behaviour of the backscatter of the different species becomes even more apparent.An influence of temperature on the backscatter is seen at both polarisations and for all three species.This is manifested in minor changes in summer and major changes in winter, when temperatures were below the freezing point.The variance of the backscatter per composite was stable over time.The band between the 25 and 75% percentiles measured consistently approximately 1.5 dB for all species and study areas.Only for the study areas in ZH (Figure 9a,c-e) the band became broader in summer 2015 due to the fewer S-1 acquisitions available in that period.
The temporal descriptors for all the study areas are summarised in Table 6 for VH and Table 7 for VV polarisation.For the beech stands with available phenological ground observations and the two oak stands, Table 8 compares the extracted break dates from the VH polarisation time series with the reported leaf emergence and fall.Figure 9 shows that at VH polarisation, areas with deciduous species followed distinctive seasonal backscatter patterns.Backscatter was higher during leaf-off and lower during leaf-on conditions.For beech, the backscatter median ranged between −13.05 and −12.54 dB over the whole time series for the different study areas.The medians of winter and summer varied substantially between the two investigated years.Backscatter in winter was higher, with higher values observed in winter 2015/16 than winter 2016/17 (Table 6).Backscatter in summer was lower, with higher medians observed for the year 2016 (see 2016/17 in Table 6).Hence, the δs between winter and summer medians also differed between the years.For 2015/16, they were about 1.9 dB, whereas for 2016/17 the difference was about 0.2 dB.Most δs were positive, giving a clear indication of the higher backscatter in winter.Only study area 1 had a negative δ in 2016/17, caused by significantly lower backscatter due to temperatures below the freezing point in winter 2017 (see dotted red 0 • C line in Figure 9a).
Compared with beech, the two oak stands showed slightly lower median backscatter over the whole times series with values of −12.97 and −13.09 dB.The winter and summer medians of both years did not differ significantly from the values observed for beech.δ values of 1.96 and 1.45 dB for 2015/16 and 0.38 and 0.60 dB for 2016/17 were observed.As for beech, this δ was positive for both oak stands in both years, indicating a seasonal change in backscatter with higher values in winter.
Comparison of the break dates and their respective ground observations for the years 2015 and 2016 for the different beech and oak stands (Table 8 and Figure 9b) shows that they match approximately with the reported leaf emergence and fall.For most study areas, the differences between the break dates and the observations were within or close to the uncertainty (±12 days) caused by the temporal resolution of the γ 0 LRW composites.In 2015, the first break dates were predicted between 3 and 11 days later compared to the leaf emergence, while in 2016, they were predicted between 14 days in advance and 5 days later.
For two study areas, the extraction of the break dates did not work very well in one year.For beech in Rafz in 2016, the second break date seems to represent the leaf emergence, as the first break date was too early for leaf emergence.For oak in Büren a.d.A. in 2015, the first break date appeared far too late (approximately 80 days) for recording leaf emergence.
The second break date, however, captured the leaf fall even better than the first had the leaf emergence.In 2015, the break dates were predicted between 12 days in advance and 1 day later, while in 2016, they were predicted 3 days in advance.Visually these results are depicted with the vertical lines in the temporal plots in Figure 9.The predicted first break dates for the oak stands are also noteworthy.They had a tendency to be a few days later than the beech stands in the same study area (Rafz) or nearby (Büren a.d.A. to Galsberg and Feiberg).In VV polarisation, the backscatter of both deciduous species was higher than at VH polarisation with a median over the whole time series ranging between −7.75 and −7.07 dB.No distinctive seasonal pattern was observed.However, looking at the plotted time series of the beech stand in Teufen (Figure 9c), after a short increase in spring, a gradual decrease in backscatter from about −6.7 to −8 dB in summer was observed in both years.

Spruce
Spruce forests exhibited seasonal backscatter behaviour in both polarisations.However, compared to the deciduous species, the behaviour was inverted.Higher backscatter was recorded in summer, and lower backscatter in winter.Generally, the spruce stands showed lower backscatter of about 1 dB compared to the deciduous species.This was consistently shown in the median backscatter values over the whole time series.In VH, the two stands had a median of −14.08 and −14.02 dB, where as in VV, their median was −8.46 dB.
As for the deciduous species, the medians of winter and summer varied across the two years.In winter, VH backscatter was lower than in summer.Lower backscatter was observed in winter 2016/17 than in winter 2015/16 (Table 6).In summer, backscatter was generally higher.In summer 2016, higher values were observed than in summer 2015 (for 2016 see 2016/17 in Table 6).In VV, backscatter differences between winter and summer and the two years followed similar patterns (Table 7).This led to negative δ values at both polarisations.Compared with VH, the δ values were higher at VV. Hence, the seasonal difference was higher in VV.This can be seen in Figure 9d,e, where the seasonal difference in backscatter appears slightly higher in VV than in VH.

Classification Performances
The demonstrated differences in backscatter between the different species was next used to classify independent validation data into forest classes.Table 9 presents the performance of the forest type classification.Table 10 and Figure 10 show the performance of the species classification.The first three columns in Table 9 and the first four in Table 10 list the number of pixels used in the training and validation sets.The classification performance was assessed using the accuracy values in the following columns.For both forest classes, producer's and user's accuracies were fairly high, resulting in a high overall accuracy and a high Cohen's κ (Table 9).The predictors' strength listed in Figure 11 showed that all twelve predictors contributed to the achieved classification result.No single predictor appeared to dominate.
As expected, compared with the forest type classification, the achieved performance was lower for the species classification (Table 10).This was mainly caused by confusion between beech and oak (see Figure 10), resulting in lower PA and UA for both classes.The performance for spruce, however, was virtually the same as for coniferous forests.Due to the limited availability of oak forests in the Canton of ZH, the number of training and validation pixels was much lower for the species classification compared with the forest type classification.For the species classification, the strength of the predictors was not as consistent as in the forest type classification (Figure 11).The δ of the year 2016 of both polarisations, but especially VH, was very important.The break dates of 2015 played a subordinate role.

Discussion
Analysis of the SAR C-band backscatter time series of forests showed that the observed backscatter behaviour was robust for multiple sites over multiple years.Generally, annual median backscatter of forests consisting of deciduous species was higher than backscatter of forests consisting of spruce.Both types of forests showed distinctive seasonal backscatter patterns.The deciduous species had higher backscatter in winter and lower in summer, while spruce exhibited the opposite behaviour.
Temporal variation in SAR data is caused by signal noise, speckle and short-and long-term environmental changes.The former three were strongly mitigated because of the methodology used in the SAR data processing.The LRW approach combines the backscatter from multiple tracks (see Figure 5) into one composite.Hence, the approach suppresses the influences of signal noise, speckle and short-term variations, such as rainfall events before single acquisitions.As a result, the annual variation can mainly be attributed to long-term environmental changes occurring within the study areas.Moreover, the decision to use terrain-flattened γ 0 T images could contribute to better results as well.All study areas (bar Büren a.d.A.) investigated in this study are in hilly terrain.

Seasonal Signature of Deciduous Species
Both deciduous species showed seasonal changes in VH backscatter behaviour.Backscatter values dropped in the spring from a high winter level to a low summer level, rising again in the autumn, resulting in positive δ values.The predominant long-term change in deciduous forests is the seasonal vegetation cycle, manifested in the annual growth of trees and presence of foliage.The presence of foliage occurred at the time with the lower backscatter values observed in the summer months: a dependence of the backscatter on leaf presence was observed [19,20].
Potential causes for the changing cross-pol scattering behaviour remain uncertain.One reason might be the changing size of the scatterers present in the tree crown.In summer, leaves substitute branches as main scatterers.C-band is primarily sensitive to small branches and only secondarily to leaves [71].Leaves have been reported to act not only as scatterers but also as attenuators for C-band [32,72].Hence, leaf emergence in spring leads to a combination of effects, resulting in lower backscatter during the leaf-on period.(a) A scatterer substitution, where the sensitivity of C-band is reduced and (b) an attenuation of the microwave energy by the leaves.Leaf fall in autumn otherwise exposes the small branches again, leading to higher backscatter.
Another possible reason for higher backscatter values during the leaf-off period is the influence of ground scattering [73].With lower attenuation in winter than in summer, ground scattering contributes more strongly to the measured backscatter.Proisy et al. [19] showed that in winter, the contribution of soil and branches were comparable, whereas in summer, the contribution of the soil dropped because of the attenuation by the leaves.
Surprisingly, VV backscatter did not show the same seasonal behaviour as VH.Although slightly higher backscatter was observed in winter, the δs show that there was a less defined seasonal signature.An interesting feature was the constant decrease in backscatter during the leaf-on season.This could be attributed to the slow drying of leaves summer [74].
Leckie and Ranson [32] did not mention a difference in temporal signature between VH and VV backscatter in their comprehensive review of SAR backscatter of forest.Potential reasons for our observed differences between the two polarisations are discussed in the following.Volume scattering becomes more important when the scatterers present are more random and complex [75].Volume scattering is also known to have a larger influence on VH than VV polarisation [31,32].Leaf emergence probably leads to a more complex arrangement of the scatterers, leading to a stronger influence on the VH than the VV backscatter.Our empirical observation of the difference between leaf-off and leaf-on conditions in C-band backscatter of deciduous forests corresponded well with the simulation made by Chuah and Tan [38].The S-1 IW acquisitions processed in our study were measured at incident angles between 31 • and 46 • .Their modelled C-band backscatter values at 50 • incident angle were different at both polarisations, but they also showed a positive difference at VH and almost no difference at VV between leaf-off and leaf-on conditions.

Break Date Extraction to Monitor Phenology
The break date extraction yielded promising predictions, bearing in mind that the time series had a temporal resolution of 24 days.Even finer details such as differences between species were observable.The observed tendency of the oak stands to have their leaves emerge later than the beech stands agrees with the common theory concerning the timing of leaf emergence of the two species [76].However, it is important to consider that the independent ground validation data were produced by observations from volunteers.Hence, the reported values have to be treated with caution.
We showed that with our approach, the phenological cycles of deciduous forests can be monitored on a landscape scale.With longer time series available in the future, it will be possible to analyse phenology trends with S-1 data.The same type of phenology trend analyses that have usually used MODIS Normalised Difference Vegetation Index (NDVI) such as White et al. [77] in North America or Hamunyela et al. [78] in Europe could be conducted.Using S-1 instead of MODIS data, an enhanced spatial and temporal resolution would be available.
The results also showed that there are still some issues using the breakpoints algorithm for the extraction of the break dates.The problems encountered for the two study areas can be explained by the weak dependency of the backscatter on foliage within each respective year.Other influences (e.g., exceptional freezing temperature) on the backscatter seemed to disturb the seasonal signature to such an extent that a smooth break date extraction was not always possible.

Seasonal Signature of Spruce
Backscatter from spruce stands did not show similarly abrupt changes.At both polarisations, their seasonal signatures were less pronounced and more gradual.Lower values were observed in winter and higher in summer, resulting in negative δ values.Because of their perpetual foliage, no strong influence of changing foliage or ground scattering was expected [39].The backscatter of coniferous forests is mainly dependent on the number of needles and small branches [11].Because the leaf area index (LAI) of coniferous forests is usually highest in July [74], the highest backscatter can therefore be expected in that month.Dostálová et al. [20] argued that in addition to the larger number of needles, the development of the understory in summer might also influence backscatter in coniferous forests.
The observed lower annual VH median backscatter of spruce stands compared to beech and oak stands was mainly due to the higher backscatter observed in winter for the latter, caused by the factors mentioned in Section 4.1.At VV, this effect was not observed to the same extent.Nevertheless, higher backscatter in winter relative to the spruce stands was apparent.This observation confirms previous studies that investigated differences in backscatter between deciduous and coniferous species [19,39,79].

Differences between the Investigated Years
Clear differences were observed within the timespan investigated.Inter-annual differences were most probably due to different meteorological conditions.A clear influence of temperature on the backscatter can be seen in Figure 9.In summer, higher VH backscatter was observed when temperatures were also higher.In winter, there was an even stronger influence of the temperature.When the temperature fell below the freezing point, the backscatter dropped significantly.This was mainly observed in January 2017, when temperatures in Switzerland were very low for a long period, but also in January 2015.Temperature is known to influence the dielectric properties of the scatterers [35].Especially when temperatures are below 0 • C, dielectric properties of the frozen scatterers change fundamentally, resulting in lower backscatter [73].
Disparities were observed in the typical summer backscatter reductions of the deciduous forests at VH. Compared to 2016, the drop in 2015 was larger, which is well depicted in the plot of the study area Feiberg (Figure 9b).This can probably be attributed to different weather conditions in the springs and summers of 2015 and 2016.Compared to 2016, temperatures, sunshine duration and the amount of precipitation in the study areas in spring 2015 were higher than normal [80,81].These favourable growing conditions might have caused a better development of the tree crowns for the deciduous species [82].According to the points mentioned in Section 4.1, a well developed crown with a high LAI would lead to a higher attenuation of the microwave energy.Thus, a larger drop would be expected and was observed.
Although precipitation can have an influence on backscatter [33], no strong influence on short-term backscatter changes was detected in our data.As mentioned above, only a general influence on the growing season coupled with other meteorological factors was observed.This may be due to the applied compositing approach, which mitigates variations caused by short-term precipitation events.
Last, a possible influence of the backscatter calibration issue in summer 2015 cannot be disregarded.To correct the wrong backscatter values, a simple correction was executed by subtracting a fixed value in the affected images [57,58].Residual errors may have been larger in these earlier data, contributing to the differences between summer 2015 and summer 2016.However, the observed disparities were not the same for the different study areas and caused no "step" in the VV data, so the likelihood of radiometric miscalibration appears negligible.

Comparisons to Other Studies
The results of this study do not contradict previously published results but extend it with an unprecedented set of C-band cross-pol observations of mixed temperate forests.Similar backscatter values were retrieved in a simulation study of sensitivity of ERS-1 C-band VV backscatter over forest [73].The non linear dependence of backscatter on temperature below the freezing point had also been simulated in that study.Proisy et al. [19] analysed ERS-1/2 C-band VV backscatter of a mixed temperate forest for the timespan between April 1994 and February 1997.They presented similar VV backscatter ranges from deciduous and coniferous stands.Their time series consisted of one image per month and had no VH measurements.Although they used a lower temporal resolution, the shape of their temporal backscatter plot of a beech stand throughout the year was similar to the one shown in this study.A recent study, analysing intra-annual S-1 data of an Austrian mixed temperate forest between February 2015 and 2016, yielded results that were in line with ours [20].The ranges of backscatter observed from deciduous and coniferous forests at VV and VH as well as the seasonal signatures shown by the two forest types at both polarisations were consistent.

Classification of Forest Types and Species
The classification results show that the observed backscatter behaviour allows segmentation into two different forest types or even into three different species.Not only the VH data but also VV contributed to achieving the good classification performances.This was surprising, as the seasonal backscatter behaviour of deciduous forests was less distinct at VV than at VH polarisation.
As expected, accuracies were lower when classifying species vs. forest type.Substantial confusion was seen between the beech and oak class, whereas in the forest type classification, they are both part of the deciduous class.The accuracies for spruce, however, were almost the same as for the coniferous class.This was not surprising as the coniferous class is mainly composed of spruce.The high accuracies were attained due to the backscatter behaviour of spruce being substantially different than that of the two deciduous species.Spruce was also the only coniferous species that could be examined.Future studies might include more species and try to use their seasonal backscatter patterns for classification approaches.
In this study, only stands with species homogeneity higher than 80% were considered for the classification.Hence, it remains unclear how the classifiers perform on more heterogeneous mixed forests.Future studies could test how the classifiers perform in mixed forests.However, when we compared the classification results to aerial images, we observed that the classifier was able to detect small patches of clustered spruce trees within a deciduous stand.Therefore, the classifiers might even outperform the aggregated forest stand-based independent ground reference.With a sample interval of 10 m, the classifiers operate at a high spatial resolution, whereas the information of the ground reference is based on larger homogenised areal extents.Therefore, many misclassifications may have been caused by the nature of the independent ground reference information.This might be verified in the future using a ground reference with a higher spatial resolution.
There was a significant time difference between the ground reference and the SAR data acquisition.The ground reference was produced in 2010, so a considerable time gap of at least 5 years was a further possible source of errors.Logging activities during that period likely degraded the classification performance.Clear-cut logged spruce areas tended to be classified as deciduous (see Figure 12), most probably due to the high ground backscatter.
It is difficult to perform an in-depth comparison to other approaches using SAR but also different remote sensing technologies to classify mixed temperate forests.Study areas, the species investigated, and the number of training and validation samples differed substantially between the studies.Nevertheless, our achieved accuracies were comparable to other studies classifying deciduous and coniferous forests using even fully polarimetric C-band SAR data [16,83,84].Compared with other technologies such as airborne laser scanning (OA of 89-96% and κ = 0.61 − 0.92) [85][86][87] or imaging spectrometer data (OA of 83-99% and κ = 0.73 − 0.98) [88][89][90], our forest type classification performance (OA of 86% and κ = 0.73) was not as competitive.This could be due to the lower spatial resolution used in this instance.

Implications of the Findings
If one is aware of the mentioned different influences on C-band backscatter, C-band data can be used for several applications.Its potential for the classification of mixed temperate forests has been presented.Thus, a repetitive classification of an area allows change detection over time.For ecological applications such as seasonal vegetation cycle monitoring, the use of C-band data will become progressively more interesting, as longer time series of data become available.Nevertheless, the capability of the data to observe subtle differences between the years has been shown using only 2.5 years of data.Because the production of S-1 data is promised until at least 2030 [91], longer time series of data should be available in the future.Using them, changes over time, such as the shift of phenophases due to climate change, could be detected and quantified.
Compared with state-of-the-art methods using optical data, SAR has valuable advantages.The sensitivity to structural elements adds an additional independent dimension of information about forests.This structural information helps to improve the understanding of natural processes in forests.Furthermore, the high temporal and spatial resolutions offered by Sentinel-1 allow a much denser sampling of data.Since the end of the commissioning phase of Sentinel-1B in September 2016 [56], the number of products has been high enough to generate LRW composites with over twice the temporal resolution of those used in this study.An LRW composite every six days would offer even more insights into the temporal evolution of the backscatter of forests.To further enhance the temporal resolution, efforts should be made to combine S-1 and RADARSAT Constellation Mission data, due for launch in 2018 [92].
With high temporal resolution available, there lies also great potential in analysing larger spatial extents.The backscatter behaviour of forests described in this study was also observed over the whole Alpine region (Lat/Lon: 43.5-49 • N/5.5-17.5 • E).Applying the same methodology as described in Section 2.3 but with a temporal resolution of the LRW composites of 12 days, descriptors were calculated within the Coordination of Information on the Environment (CORINE) land cover classes for deciduous and coniferous forest [93].Plotted in the temporal domain, the same seasonal patterns are apparent at both polarisations and for both forest types, as shown in Figure 13.Future research could focus on more detailed analysis of the backscatter and the understanding of scattering processes.An in-depth study of the observed summer drop in VH backscatter might result in interesting new insights, especially with a higher temporal resolution.Differences, for example, between the deciduous species are hypothesized, but could not be tested in this study.A dependence of the magnitude of the observed drop to LAI could also be tested.If the drop in backscatter is caused by the attenuation of the leaves, a dependence is credible.This could be tested by correlating VH backscatter time series with periodic LAI estimations derived from hemispherical images of deciduous forest stands in spring or autumn [94].
This leads to the question of to what extent the findings of this study could be transferred to other ecosystems for time series analysis.In tropical or boreal ecosystems, where the use of SAR data is already well established, the use of multitemporal C-band data might contribute valuable new insights into the temporal evolution of these ecosystems.

Conclusions
This study showed that multitemporal C-band radiometrically calibrated and flattened backscatter data have the potential to supplement optical remote sensing data for use in forest mapping or ecological applications.Forests consisting of broadleaved, deciduous species, such as beech and oak, produce stronger backscatter in both polarisations than forests consisting of spruce.Typical seasonal phenological signatures were also observed within forested areas.In VH, areas consisting of deciduous forests showed a consistent and distinct intra-annual backscatter pattern with higher backscatter in winter and lower in summer.In coniferous forests, however, the opposite behaviour with lower backscatter in winter and higher in summer was observed.The same behaviour, but even more distinctive, was observed in coniferous forests at VV, whereas deciduous forests showed almost no seasonal cycle in VV.
These temporal signatures simplify the classification of mixed temperate forests because of the striking differences between the temporal signatures of deciduous vs. coniferous forests.Good classification performances were achieved for the forest types (OA of 86% and κ = 0.73) and for three different species (OA of 72% and κ = 0.58).The strong seasonal signal in VH allows the monitoring of phenology in deciduous forests.The dates independently reported for leaf emergence and leaf fall in three different study areas were satisfyingly matched.Almost all extracted break dates were within the error margin (±12 days) of the available temporal sampling.
These findings show that using multitemporal C-band data offers several new possibilities.Periodic and cost-effective mapping of mixed temperate forests will ease forest management to a great extent.In the near future, longer time series of S-1 and RADARSAT Constellation Mission data over several years will be available with high temporal resolution.Using them, monitoring at different sets of spatial and temporal scales of different ecosystems will be possible.These might lead to interesting ecological insights, such as an assessment of the shift of phenophases or a quantification of changes in species composition due to climate change.

Figure 1 .
Figure 1.The five investigated main study areas are located within the yellow rectangles in the Cantons of Zurich (ZH) and Bern (BE) in northern Switzerland.For some analyses, the entire forested area within the Canton of Zurich, depicted in the green rectangle, was used.The background aerial images were provided by the Swiss Federal Office of Topography swisstopo[44].Reproduced by permission of swisstopo (BA17116).Swiss map coordinates: CH1903+/LV95.

Figure 6 .
Figure 6.RGB-overlays of (a) γ 0 E images and (b) γ 0 LRW composites.The assigned bands are Red and Blue = summer and leaf-on (γ 0 E : 6 June 2016, γ 0 LRW : 24 May-16 June 2016) , Green = end of winter and leaf-off (γ 0 E : 16 March 2017, γ 0 LRW : 2 March-25 March 2017).The range of the backscatter values (−20 to −10 dB) was set identically for all channels in both RGB-overlays.Non-forested regions are masked out.In the lower row, maps of the main influencing factors on the backscatter are depicted.The slope map (c) of the area is shown as well as (d) the forest type map with the classes 'deciduous', 'coniferous' and 'not defined', consisting of all non-forested areas and areas with less than 50% homogenous species (Data source: GIS-ZH).The top row images contain modified Copernicus Sentinel data [2016/2017].Swiss map coordinates: CH1903+/LV95.

Figure 7 .
Figure 7. Correlation density plots comparing end-of-winter and summer VH backscatter of (a) coniferous and (b) deciduous forest in the whole of the Canton of Zurich for the γ 0LRW composites.The number of pixels (n) was considerably higher for coniferous forests.A higher Pearson correlation coefficient (ρ) was achieved for deciduous forest.

Figure 8 .
Figure 8. Histograms of the difference between end-of-winter and summer VH backscatter of the two forest types in the whole of the Canton of Zurich for (a) the γ 0 E images and (b) the γ 0 LRW composites.n = 36,621 for deciduous forest, 152,930 for coniferous forest.

Figure 9 .
Figure 9. Temporal plots of species backscatter from γ 0 LRW composites for selected study areas.(a-e) Blue shows the median, green corresponds to the 25 and 75% percentiles.The mean temperature for the 24 day period corresponding to each γ 0 LRW composite is shown in solid red.The horizontal dotted red lines mark the 0 • C line: forest backscatter was generally reduced when the temperature dipped below 0 • C. The solid vertical lines show phenological ground observations.Vertical green lines indicate leaf emergence for beech, and needle shoot for spruce.Red lines show the leaf fall date for beech.Dashed lines correspond to the backscatter derived break dates.The black bars on top of the figure display the periods used for the calculation of winter (1 December-15 March) and summer (1 June-15 September) medians.

Figure 10 .Figure 11 .
Figure 10.Map of the species classification for a subset of the study area in Rafz (ZH).It can be observed that the segmentation into the forest types 'deciduous' (Beech, Oak) and 'coniferous' (Spruce) is possible to a high accuracy.Confusion is seen between the two deciduous species shown in green tones.The background aerial image SWISSIMAGE was provided by the Swiss Federal Office of Topography swisstopo [44].Reproduced by permission of swisstopo (BA17116).Swiss map coordinates: CH1903+/LV95.

Figure 12 .
Figure 12.Map of the species classification for a subset of the study area in Rafz (ZH).The impact of logging activities on the classification result can be observed.A former logged spruce stand was mainly classified as beech or oak.The background aerial image SWISSIMAGE was provided by the Swiss Federal Office of Topography swisstopo [44].Reproduced by permission of swisstopo (BA17116).Swiss map coordinates: CH1903+/LV95. 06

Figure 13 .
Figure 13.Temporal plots of the per γ 0 LRW composite calculated spatial descriptors median (blue), 25 and 75% percentiles (green) for decidous and coniferous forests of the Alpine region (Lat/Lon: 43.5-49 • N/5.5-17.5 • E).The temporal resolution of the γ 0 LRW composites was 12 days and the time period depicted ranges from October 2014 to May 2017.The number of investigated pixels (n) was comparable for both forest types.

Table 1 .
The investigated study areas are distributed in northern Switzerland in the Cantons (Ca) of Zurich (ZH) and Bern (BE).The listed values in the columns elevation, slope and diameter at breast height (DBH) are all mean values from the respective study area.The core extent describes the investigated extent after the applied erosion.

Table 3 .
SAR processing schemes and backscatter products used in this study.

No. of IW acquisitions contributing to each 24 day composite
Classification scheme.This framework was used for both forest type and species classifications.

Table 4 .
Descriptors calculated within single composites and within the time series.

Table 5 .
Median ( x), 25 and 75% percentiles (Q1 & Q3) C-band VH and VV backscatter (γ 0 E images) for deciduous and coniferous forest stands of the whole of the Canton of Zurich.The investigated stands had a species homogeneity of at least 80%.

Table 6 .
VH polarisation time series descriptors for the five study areas (SA) with their respective occurring species (S).Column N shows the number of γ 0 LRW composites and Px the mean number of processed pixels.xTo,Q1To and Q3To are the median, 25 and 75% percentiles over the whole time series.γ0win, γ0 sum and δ show the medians for winter (1 December-15 March), summer (1 June-15 September) and the difference between the two.These descriptors were calculated for the years 2015/16 and 2016/17.The displayed backscatter values are in dB.

Table 7 .
VV polarisation time series descriptors for the five study areas (SA) with their respective occurring species (S).Column N shows the number of γ 0 LRW composites and Px the mean number of processed pixels.xTo, Q1To and Q3To are the median, 25 and 75% percentiles over the whole time series.γ0 June-15 September) and the difference between the two.These descriptors were calculated for the years 2015/16 and 2016/17.The displayed backscatter values are in dB.

Table 8 .
Comparison of the break dates (BD) extracted from the VH polarisation time series with the respective leaf emergence (LE) and leaf fall (LF) ground observation for the different study areas.For both years, the values are displayed in day of the year.The listed uncertainty is due to the compositing timespan of 24 days.

Table 9 .
Classification result for the two forest types.The mean number of pixels used for training in the 40 RF classifications for the two forest types is listed with the mean confusion matrix, producer's accuracy (PA), user's accuracy (UA), overall accuracy (OA) and Cohen's κ of the 40 different RF models.

Table 10 .
Classification result for the three different species.The mean number of pixels used for training in the 40 RF classifications for the three species is listed with the mean confusion matrix, producer's accuracy (PA), user's accuracy (UA), overall accuracy (OA) and Cohen's κ of the 40 different RF models.