Can SAR Data Predict Harvesting Residues Mass?

1. Introduction

At the global level there has been an increasing demand for large-scale observations at high temporal frequency to monitor environmental changes and manage natural resources effectively, such as forests [1]. The increasing use of satellite-based technology offers a solution by enabling frequent, wide-area monitoring of forests, their status and biomass, along with other ecological parameters [2,3]. Traditional methods of forest monitoring, such as ground surveys and aerial photography, are labor-intensive and limited in spatial and temporal coverage [4]. Therefore, obtaining accurate and timely data is essential for informed decision-making, especially under the particularly crucial context of climate change and sustainable forest management.

Interferometric Synthetic Aperture Radar (InSAR) is a powerful remote sensing technique that generates high-resolution maps of surface deformation, topography, and other characteristics by analyzing the phase difference between two or more SAR images of the same area taken from slightly different positions or at different times [5]. Differential Interferometric SAR (DInSAR), an extension of InSAR, enhances this capability by comparing multiple SAR images over time to detect changes in the Earth's surface, providing detailed information on ground deformation and other dynamic processes [6].

DInSAR has become a valuable tool in environmental and forestry research, particularly since the BIOMASS mission, there is increasing interest in the use of this information for estimating forest biomass [7,8]. By capturing subtle changes in the forest canopy and ground structure, this technique can provide accurate measurements of forest height and volume [9], which are crucial for biomass estimation. Researchers leverage the sensitivity of SAR to the physical and dielectric properties of the observation objects to infer biomass from changes in signal phase and amplitude [10]. This method has been shown to be effective in various forest types and conditions, e.g. in the boreal forest of Canada [11] and Norway [12], in Brazilian Amazon [13] and tropical forests [14], making it an important technique for large-scale biomass monitoring and carbon stock assessment.

Harvesting residues, mainly consisting in branches and other non-commercial parts left on the ground after logging operations, are a significant component of forest ecosystems [15]. The role, presence and management of residues is a key factor for several reasons: they actively contribute to nutrient cycling within the forest [16], especially when it comes to influence Carbon and Nitrogen [17,18]; more to that, they also affect soil physical structure [19], biodiversity [20,21,22] and overall ecosystem health [23]. On the other hand, to reduce the impact and frequency of wildfires, the management of residues has helped to reduce the fire risk, as accumulated biomass can serve as fuel for wildfires [24]. At the same time, harvesting residues are a potential source of biomass for bioenergy production, an increasingly important element of renewable energy strategies [25,26].

Forest management procedures, such as clear-cutting, selective logging, and thinning, also depending on the machines and harvesting systems adopted, can generate varying amounts of harvesting residues [27]. The type and intensity of management practices influence the distribution and composition of residues, affecting subsequent decomposition rates and nutrient cycling.

Previous research on harvesting residues has primarily focused on field-based methods and statistical models to estimate the quantity and distribution of residues. Studies have explored various approaches at various scales and levels. Both Woo et al. [28] and Li et al. [29] have explored the possibility of integrating direct measurement from forest machinery during logging operations with allometric equations. Windrim et al. [30], on the other hand, developed an automated mapping of woody debris over recently harvested areas using UAV-borne data, high-resolution imagery and machine learning, and remote sensing techniques. However, these methods often face challenges in terms of scalability and accuracy, particularly over large and diverse forest landscapes.

The aim of this study is to understand the possibility to predict the quantity of harvesting residues over clear-felled areas for SAR data, in particular adopting InSAR and DInSAR techniques. The hypothesis is that along the observation period there will be little or no change in the residues spatial distribution (i.e., the geometry of the object), whereas most likely there will be changes in their properties, in particular the water content, that would be reflected in a variation of the dielectric constant of the woody material. Therefore only the change in moisture content will be investigated.

2. Materials and Methods

2.1. Study Area

The study area is part of an industrial forest plantation located in the KwaZulu-Natal province in South Africa (Figure 1). The complex topography of the province grants a varied and verdant climate: the area is characterized by a dry-winter humid subtropical climate, with warm summers and frequent rainfall. Winters are dry with high diurnal temperature variation, with possible light air frosts (South African Weather Service). The site, with a total surface of 17.54 ha, is divided into two parts which were both felled in the same period: the southern part has a surface of 7.34 ha, whereas the northern part is larger with a surface of 10.24 ha.

The species planted in the study area was Mexican weeping pine (Pinus patula Schiede ex Schltdl. & Cham.), hereafter identified as Pine. The harvesting operations were performed during February 2023, where the timber was mechanically felled by a harvester machine and then extracted and transported by forwarder to depot. The felling machine adopted was a Tigercat LH822D equipped with Log Max 7000XT head working with cut-to-length (CTL) technique, whereas the forwarder adopted was a Tigercat 1075C.

To better understand and inspect the area, the whole site was flown using a DJI Mavic Air 2S with a flying altitude of 65 m with 80% forward and 70% lateral overlap. The flight was performed in scattered cloud cover conditions, with Ground Control Points (GCP) positioned and recorded as well. The images were handled using a Structure-from-Motion (SfM) technique with Agisoft Metashape®, from which an orthomosaic was derived.

2.2. Field Data: Sampling Design and Residue Mass Estimation

The field measurements, used to estimate the harvesting residues mass, were first carried out in March 2023 in the southern part of the site, after the majority of the forwarding operations were completed and before the site was burnt and replanted at the end of October of the same year.

The sampling technique adopted was a line intercept sampling (LIS) method [31,32,33], used to estimate volumes and weights of downed woody material over completely clear-felled areas [34,35]. This technique considers residues on site to be randomly spread, with a random orientation under the sampling line and positioned horizontally, their circular shape, and the normal distribution of diameters.

Following the guidelines from Rizzolo [36], under each sampling line the diameter of every piece of woody debris was registered and classified in time lag classes (Table A 1). Each plot was composed of three sampling lines of 20 m each and organized as a perfect triangle with random orientation [37] (Figure 2); the number of triangles accounted for a total of 13 plots. For each triangle, the vertex positions were recorded using a GNSS system. The time lag division of woody material refers to the time required (in hours) for their moisture condition to fluctuate, Large woody debris, for example, will inherently take longer to dry out than finer material. According to this division, material finer than 76 mm in diameter belongs to fine woody debris (FWD) and larger material to coarse woody debris (CWD). The length has also been recorded for CWD larger than 203 mm.

Furthermore, 14 additional virtual plots were added across the area (using visual inspection in QGIS) to increase the information available for the training of the model, hereafter named “interpreted plots”, of which 8 were positioned in the northern part of the compartment. For each interpreted plot, the sampling design was similar to what was performed for the field plots: the plot corresponded to a perfect triangle with a side of 20 m and residues were sampled and digitized over the lines. The diameter for each piece was recorded and sorted in classes according to Table A 1. In particular, considering the resolution of the orthophoto (2.5 cm/pixel), it was possible to record residues with diameter with that dimension or higher.

The residue quantity estimation (Mg·ha^-1) was performed for all categories of residues: for 1h, 10h, 100h, and 1000h classes, the Brown’s formula was used [34,35], described in Equation 1. Whereas for bigger elements, the ones falling in the 1000h+ class, the estimator was computed using the Woodall formula [35], expressed by Equation 3.

$${\widehat{Y}}_{FWD, CWD}=\left({\displaystyle \frac{1.234 ·n·{\overline{d}}^2·SG ·c·a}{\sum L}}\right)·k_{decay}·10.000$$

(1)

$$c=\sqrt{1+{\left({\displaystyle \frac{Slop{e}_{\%}}{100}}\right)}^2}$$

(2)

$${\widehat{Y}}_{1000h+}={\displaystyle \frac{\pi }{2L}}·SG·{{\displaystyle \sum }}_{i=1}^n\left({\displaystyle \raisebox{1ex}{$y_i$}\!\left/ \!\raisebox{-1ex}{$l_i$}\right.}\right)·k_{decay}·10.000$$

(3)

Where: 1.234 is a computational constant to convert the volume to m³/ha; n is the number of elements for each class; d is the average squared diameter for the class; c is the corrected slope (Equation 2); a is the correction coefficient for the position of the elements, equal to 1.13 for FWD and equal to 1 for CWD; SG is the specific gravity of the wood, in this case an average value of 0.575 was adopted (575 kg m^-3); L is the length of the sampling line(s); k_decay is the decay coefficient as described in [35], in this case k_decay = 1 since the material was freshly cut; 10,000 are the square meters in 1 ha; y_i is the volume for the single CWD piece and l_i is the length of the piece.

2.3. Satellite Data and Interferometric Processing

The Sentinel-1 mission included a constellation of two synchronous satellites where currently only one is still operating properly, equipped with a C-band synthetic aperture radar (SAR), enabling them to get acquisitions regardless of the weather. For this study, 4 dual polarization Sentinel-1 SAR C-band acquisitions were downloaded from the Alaska Satellite Facility (Copernicus Sentinel data 2023, processed by ESA), date of download: 13/07/2024. The images were acquired in the interferometric wide swath mode, containing data with a 250 km swath at 5 m by 20 m spatial resolution (single look, SLC), and captured in three sub-swaths using Terrain Observation with Progressive Scans SAR (TOPSAR). Each SLC-IW image contains single-look complex data, from which amplitude and phase can be computed for each polarization: the amplitude describing the strength of the signal and the phase the number of cycles modulo 2Pi of oscillation that the wave executes between the radar and the surface and back again.

In this study, we used a technique known as SAR Interferometry (InSAR), in which two images acquired from slightly different positions (i.e., separated by a spatial baseline) are used to obtain a three-dimensional image of the Earth’s surface.

In order to measure the changes of the surface, images need to be acquired at different times (i.e., using a temporal baseline), a technique denoted Differential SAR Interferometry (DInSAR). In this case, a reference image (“master”) was selected (in our case, 12/04/2023) to compare all following changes, whereas the others were defined as “slave” images (06/05/2023, 17/07/2023 and 09/10/2023) (Figure 3).

2.4. Pre-Processing of Satellite Images

The SAR images were first processed using the Sentinel Application Platform (SNAP) from ESA (ESA Copernicus Hub). First of all, the images were inspected and the first sub-swath was selected, as it alone contained the study area. Moreover, among the two polarizations (VH and VV), the proposed methodology relied only on VV-polarization since it has been proved to be more sensitive to surface moisture [38,39], whereas the cross-polarization signal is more sensitive to changes in the geometry and physical properties, such as the roughness.

The main steps of the following methodology are summarized in Figure 4. To better fulfill the aims of this research, images were coupled in short (April-May), medium (April-July) and long-period (April-October). For each pair of SLC-IW single-pol images, three explanatory variables were derived: the phase, the amplitude and the coherence using the interferometric approach and the intensity of the signal expressing the backscatter. The coherence provides a measure of similarity for each pixel ranging from 0 to 1, where 0 indicates that change did happen and 1 where no change has happened. First, the images were coregistered, their orbit corrected and stacked; then, the interferogram was computed. After the Topographic Phase removal, the interferogram was multilooked (number of looks = 3) and Goldstein phase filtered (coherence threshold = 0.4). After this, the DInSAR phase (∆φ), amplitude and coherence images were geocoded and exported.

To compute the backscatter, the VV-pol images were handled individually. After orbit and radiometric correction, the images underwent multilooking (number of looks = 3) and speckle filtering (sigma = 0.9). The images were terrain corrected and exported as γ⁰ images describing the backscatter coefficient normalized by the illuminated area projected into the look direction.

2.5. Feature Extraction and Prediction Model for Residue Mass Estimation

The following steps in the elaboration were performed using R [40] and RStudio environment (R version 4.4.1; RStudio 2024.04.2+764 "Chocolate Cosmos"). After the pre-processing, the ∆φ, amplitude, coherence and γ⁰ images were stacked together and handled as raster, using the terra package [41]. The plots were paired with the mean quantity value obtained by averaging the estimations from their sample lines. The shapefile containing this information was used to extract the mean value from the stacked backscatter and then converted to data frame.

The data frame was used to feed two models: a generalized linear model (GLM) and a Random Forest model (RF). In the first case, the application of a linear model also had the intention to find possible correlations between the variables and to underline possibly significant variables available in the data frame. The RF model used is a machine learning algorithm built with the R package randomForest [42], with forest size of 500 trees (ntree = 500). All variables were used to initially feed the model, adequately setting the split decision parameter (mtry), and training and testing was performed using a leave-one-out cross-validation (LOOCV) procedure. The RF model was then iteratively adjusted by excluding the variables that reported a negative importance score.

2.6. Coherence Difference Analysis

To better gain insights on the results, a separate analysis was performed by considering the coherence values, comparing each image pair and by computing the differences between the pairings. Coherence (ϒ) values are calculated at pixel level and range from 0 to 1; a value of 1 indicates no alteration in the scattering properties between images [43]. However, when the observed surface changes, the complex backscatter is impacted causing a decrease in coherence, which is known as decorrelation. This approach moves from already established coherence difference analysis (CDA) methods [44] but still follows the same purpose, to define and to investigate the source of decorrelation. Consistent with [45] the interferometric coherence is affected by three main factors of decorrelation: radar thermal noise decorrelation (ϒ_N), geometric decorrelation (ϒ_G), and temporal decorrelation (ϒ_T). Considering that all images were acquired by the same antenna, it is reasonable to assume that the variation in thermal noise between images will not influence the coherence decorrelation [43] The geometric decorrelation, on the other hand, is highly influenced by the perpendicular baseline; if this parameter is major than the critical value, then ϒ_G becomes 0, causing maximum decorrelation [46]. Finally, the temporal decorrelation is the one that would contribute the most since the loss in coherence is mainly due to changes in the object's properties (i.e., geometry structure and dielectric properties).

2.7. Validation Assessment

Both models were validated by using the field plots as reference by comparing the predicted values against the estimates obtained through volume calculation. In this case, two indicators were considered: the squared Pearson index (R²) to evaluate the correlation between the field values and the predicted, and the Root Mean Squared Error (RMSE) to assess the model performance. At the end, the model was used to predict the residues mass distribution throughout the study site.

To further evaluate the goodness of the produced residues estimations and mapping, the field plots were also visually inspected, and the diameter of the material was registered, with the residues mass computed. After this, the bias and RMSE were computed.

3. Results

3.1. Field and Interpreted Plots Residues Quantity

The estimated quantities per hectare and per plot surface are reported in Table 1. Considering the average values, it can be shown that the Interpreted plots have higher values in both quantities per hectare and per plot surface (+113% compared to the field measurements). Nonetheless, if considering the total mass for each plot, the two series report comparable information (Table A 2 and Table A 3).

3.2. Coherence Difference Analysis

The coherence differences were computed for each image pair and between pairs (Figure 5). For the first pair (April-May), the coherence values were estimated between 0.036 and 0.793 (mean = 0.302; standard deviation = 0.13), whereas for the second pair (April-July), the values ranged from a minimum of 0.031 to 0.779 (mean = 0.258; standard deviation = 0.102). The last pair, April-October, showed coherence values between 0.037 and 0.679 (mean = 0.249; standard deviation = 0.097). The observed difference in the statistics results over time (short, medium and long period) denote an important decorrelation process between acquisitions. The coherence difference maps showed values between −0.531 and 0.555 (mean = −0.044; standard deviation = 0.152) for Figure 5C map, while the values reported in Figure 5E ranged from −0.651 to 0.482 (mean = −0.053; standard deviation = 0.156).

3.3. Prediction Models and Validation Assessment

The results for the GLM and RF models are shown in Figure 6, with the distinction between field plots and interpreted plots. In this case, the GLM model have used all the variables (Table A 5) with only one showing significant contributions at the end (p < 0.1), i.e. the phase from the April-May interferogram. The variables feeding the RF model, on the other hand, were filtered based on their importance score using only variables with a positive score (Figure A 1): in this case, the variables reaching the highest scores are the amplitude from the April-October interferogram and the phase from April-May, with also the contribution of the Coherence April-May and the γ⁰ from July; the rest of the variables show the same importance score. Overall, the GLM model performed better in terms of both R² and RMSE compared to RF. In the validation process, i.e. by considering only the field plots, the GLM model outperformed RF again (Figure 7).

For the bias assessment, the field plots have been inspected and interpreted, with the resulting volumes reported in Table A 5 and the error indices per single plot values (Mg) in Table 2. In this case, the errors computed underline how the residues estimated with the visual interpretation method are biased, on average, with a -0.05 Mg when considering just the CWD components (1000h and 1000h+ classes), and 0.14 Mg when considering all classes of material. The RMSE is lower than the one provided by both the GLM and RF models.

The representation of the residue distribution over the study area is then reported in Figure 8, with respect to the single plots quantity. In this case, the heat map built with the GLM model tends to represent all the spectrum of values, whereas the RF model tends to represent values converging towards the mean of the distribution. Regardless of the accuracy of the model used to depict such spatial representation, this tendency can be seen by using the plots values as references.

4. Discussion

This paper aimed to evaluate the possibility of using SAR data techniques (InSAR and DInSAR) to quantify the mass of harvesting residues left in clear-felled areas. The proposed approach was to compare two prediction models and evaluate their goodness-of-fit (R²), RMSE, and bias. This paper offers an innovative approach to estimate residues quantities more conveniently than physically weighing their mass.

The coherence difference analysis assumed that the presence of harvesting residues in clear-felled areas would be less stable over time than the surrounding solid soil (due to decomposition processes, residue pile volumes gradually collapsing to a solid ground layer etc.). Hence, residues would decorrelate more (measured as decreasing coherence) as time passes. The variation in dielectric constant was assumed to be lower for the surrounding clear-cut plantation, and hence present only smaller changes in the backscatter between SAR images. The perpendicular baselines ranged from -117 m to 80 m, which is a bit large for the sensitivity we aimed for. Nevertheless, the loss of coherence was useful in our predictions.

In the short term , the coherence exhibits higher values inside the harvesting area with more variations outside the perimeter, where the trees are still standing, as also reported by Akbari et al. [47], whereas the coherence gradually decreases with time, resulting in values close to 0. Considering the absence of inclement weather conditions throughout the time of observation, the changes in coherence in the short term are likely due to changes in volume scattering and changes of the dielectric constant throughout the drying process; changes in the surrounding trees are more likely related to fluctuations in the crown water content, due to the SAR operating at C-band [48]. These assumptions are also supported by the coherence differences , where lower changes are noticed outside the logged area (Figures 5C and 5E). Some of these findings were also noted when developing the linear prediction models, where the relation and contribution of various variables were studied. For example, in the GLM model the variables with the highest significance came from the short period interferogram (i.e., the April-May interferogram), that is also exhibiting higher variation in coherence values. On the other hand, the RF model prefers the phase difference between April and May, but also the amplitude from April-October (long term) and the γ⁰ from July.

The magnitude of residue masses estimated through the field survey are comparable with the values available in the literature for the species considered [37,49]. The higher residue mass values of the interpreted plots (+113% compared to the field plots) can be due to the influence of several factors. First of all, the interpretation and visual inspection of the logging area by the means of aerial photogrammetry: although available with fairly high resolution (2.5 cm/pixel), still the orthomosaic used for the visual assessment of the residues did not allow for the inspection and registration of the smaller woody material, with diameters lower than 3 cm (i.e., 1h and 10h classes). This could be reflected in the average quantities per plot that are higher in the interpreted plots compared to the field ones. Moreover, whereas the field survey also considers stumps, the visual inspections do not, therefore further reducing the CWD material estimations. A similar issue can be found also in [30], although using a different system for manual annotation. Second, the presence of residue piles in the area does not allow for an optimal visual annotation of material from QGIS, since it would be possible to register only the diameter of the elements on the very top of the pile, whereas with the field survey some material also from the deeper layer can be registered. Therefore, the volumes or masses considered results in different values, also for the same plots, as exposed by the bias evaluation (-0.05 Mg and 0.14 Mg, for CWD and FWD+CWD considered material).

Nevertheless, there are different factors and limitations to consider for this study and the proposed methodology. The main principle of the DInSAR technique adopted is that it allows for the retrieval of any changes on the surface of the area investigated, through the scatter of the beam signal; the rougher the object, the more scatter it should produce.

• Overall, this study investigated only the use of VV-polarization acquisitions, based on previous literature, considering any major detectable changes due to moisture variation in the residues and piles;
• A second general limitation is regarding the object of the investigation, the harvesting residues. Although not individually detectable using the SAR pixel size (10 x 10 m), they can be addressed in piles, providing a large (on average, 5 meters wide and 1.5 meters tall, at maximum), complex and rough object, when performing this kind of research. However, obtaining precise information of pile masses and volumes is still challenging, both with field measurements and remote sensing-based techniques [49,50,51];
• A topic more related to the prediction models is the number of available plots used to feed the model and its design. To increase the R² of the results and their robustness, it would have been necessary to have more field observations for both GLM and RF models. In fact, both models can be easily exposed to dataset with a limited number of observations, in particular RF [52], therefore more in-field observation are preferrable;
• The analysis of coherence provided insightful information about the signal scattering and its changes throughout the monitored period. Although minimized, different sources of decorrelation might have still played a role in the obtained results, especially the geometric decorrelation. Moreover, to properly address the volume scattering, the proposed methodology featured also the backscatter signal γ⁰, but saw limited contribution from the majority of them (Figure A 1). In this case, the possible addition of VH-polarization backscatter could add complexity to the response, delivering a more comprehensive output.

5. Conclusions

The paper evaluated the possibility for SAR data to predict the quantity of harvesting residues over clear-felled areas, performed over a study area of a pine plantation in South Africa. The models predictions provided results with an R² index of 0.47 and 0.13 for GLM and RF, respectively, when using all plots, and 0.23 and 0.006 when using only the field plots. The RMSE stabilized between 0.20-0.28 Mg. Overall, the bias regarding the interpretation of the woody material ranged between -0.05 and 0.14 Mg. Nevertheless, we can positively affirm that, even at this stage and with the limitations listed above, it is possible to derive harvesting residues mass from SAR data handled with InSAR and DInSAR techniques.

This information could be suitable to forest owners and managers for the possible retrieval of material to convoy to bio-energy producing plants, or to provide an indirect estimation of possible C-stock in woody material to help industrial companies to account for carbon or reduce their impact on biogeochemical cycles and nutrient fluxes in the re-establishment of timber plantations and the use of fertilizers.

Future research should address the necessity of defining a suitable number of plots or sampling design to derive residues information and train/test any prediction model. Moreover, the possibility of having different observations on different sites at the same time could be interesting to explore to address local variability between sites and/or species.