DETECTING TROPICAL FOREST DEGRADATION USING OPTICAL SATELLITE DATA: AN EXPERIMENT IN PERU SHOW TEXTURE AT 3 M GIVES BEST RESULTS

Forest degradation is known to be widespread in the tropics, but is currently very poorly 1 mapped, in part because there is little quantitative data on which satellite sensor characteristics 2 and analysis methods are best at detecting it. To improve this, we used data from the Tropical 3 Forest Degradation Experiment (FODEX) plots in the southern Peruvian Amazon, where different 4 numbers of trees had been removed from four 1 ha forest plots, carefully inventoried by hand and 5 Terrestrial Laser Scanning before and after the logging to give a range of biomass change (∆AGB) 6 values. We conducted a comparative study of six multispectral optical satellite sensors (WorldView7 3, SkySat, SPOT-7, PlanetScope, Sentinel-2 and Landsat 8) at 0.3 – 30 m spatial resolution, to find the 8 best combination of sensor and remote sensing indicator for change detection. Spectral reflectance, 9 the Normalized Difference Vegetation Index (NDVI) and texture parameters were extracted after 10 radiometric calibration and image preprocessing. The strength of the relationships between the 11 change in these values and field-measured ∆AGB (computed in % ha−1) was analysed. The 12 results demonstrate that: (a) texture measures correlates more with ∆AGB than simple spectral 13 parameters; (b) the strongest correlations are achieved for those sensors with spatial resolutions 14 in the intermediate range (1.5 10 m), with finer or coarser resolutions producing worse results, 15 and (c) when texture is computed using a moving square window ranging between 9 14 m 16 in length. Maps predicting ∆AGB showed very promising results using a NIR-derived texture 17 parameter for 3 m resolution PlanetScope (R2 = 0.97 and RMSE = 1.80 % ha−1), followed by 1.5 18 m SPOT-7 (R2 = 0.74 and RMSE = 5.25 % ha−1) and 10 m Sentinel-2 (R2 = 0.71 and RMSE = 5.55 19 % ha−1). Texture models derived from 0.3 m WorldView-3 improved with increasing window 20 size, with highest R2 of 0.62 and RMSE = 6.35 % ha−1 for a window of 14 m in length. The 21 degradation in our field plots is invisible to the 30 m resolution Landsat data. Our findings imply 22 that, at least for lowland Peru, low-medium intensity disturbance can be detected best in optical 23 wavelengths using a texture measure derived from 3 m PlanetScope data. That such data are being 24 collected daily, and currently released free as monthly mosaics over tropical forests as part of the 25 Norway’s International Climate and Forest Initiative (NICFI), is excellent news for monitoring 26 such degradation. 27


Introduction
Tropical forests play a complex and essential role in modulating the global climate 31 system by removing carbon dioxide from the atmosphere through photosynthesis. They 32 are home to a rich biodiversity and provide an essential source of livelihood and socio- 33 cultural identity to local communities. However, pressures from extractive industries, 34 commodity-driven agricultural expansion and fires are increasingly threatening the 35 integrity of these forests [1,2]. As a result of factors such as deforestation, forest degrada- 36 tion and climate change, the ability of tropical forests to sequester atmospheric carbon 37 might be reduced or reversed in the near future [3,4] 38 Contrary to deforestation, which occurs when a forest has sufficient trees removed 39 that it no longer meets the local definition of forest [5], degradation is harder to define 40 [6]. In general, forest degradation is considered to occur when humans cause partial 41 deforestation, with more than 10-30% of forest cover remaining [7]. Half of tropical 42 forest degradation is thought to be driven by selective logging, followed by fuelwood 43 and charcoal production, fires and livestock grazing [8]. Studies aiming at quantifying 44 the carbon emissions and extent of degradation show that it is diffuse across the tropics 45 and affects land portions similar to, or even larger, than deforestation [9,10]. Since 46 degradation precedes deforestation, controlling the former may prevent the development 47 of the latter [6]. Monitoring forest degradation while protecting existing old-growth, 48 natural forests, has longer term climate and ecological benefits than planting new trees 49 [10][11][12][13][14], typically when forest lands are managed collectively [15][16][17][18][19] and the social 50 drivers of deforestation are addressed [20]. It is worth noting that the monitoring effort 51 should not imply the commodification of carbon stocks [21,22]. An abundant literature 52 has shown that market-based mechanisms reinforce exclusionary politics, such as the 53 centralization of forest tenure and the exclusion of certain actors, practices, knowledge 54 and claims to resources [23][24][25][26][27], which in turn exacerbates forest loss [25,28]. Addressing 55 deforestation and forest degradation requires trade-offs between economic, ecological 56 and social dimensions [29] and the recognition of the political and power dimensions 57 behind the interventions; an aspect that carbon finance, with its unidimensional approach 58 to forest management, tends to overlook, to the detriment of forests and forest-dependent 59 communities [30]. 60 While it has been possible to map deforestation even with moderate resolution (30 -61 250 m) satellite imagery [31,32], disturbances at finer scales, as those caused by selective 62 logging or other degradation processes, remain challenging to detect [2]. High resolution 63 optical data (e.g., WorldView-3) are cost prohibitive and are not always immediately 64 available in all regions. On the other hand, Synthetic Aperture Radar (SAR) systems 65 do not perform well in dense canopies due to biomass saturation effects [33]. Airborne 66 LiDAR has now gained widespread adoption in forestry applications because of its 67 ability to produce 3D point clouds with centimeter accuracy [34,35]. However, LiDAR 68 collection campaigns are very expensive and data availability is drastically limited in 69 space and time [36]. 70 In the last decade, the progressive release of remote sensing data that is both tem- 71 porally dense and high-resolution (≤10 m) has opened new frontiers for mapping and 72 monitoring changes in forest cover in usually remote and inaccessible areas. The use 73 of dense time series to detect forest disturbance events is particularly important in the 74 tropics, as rapid vegetation turnover can remove the transient signs of the disturbance 75 events rapidly [37], and cloud cover often masks images, meaning that many attempts 76 may be necessary before the trees are seen from space. Freely available, high-resolution 77 optical multispectral imagery from the European Space Agency (ESA)'s Sentinel-2 mis-78 sion and the monthly PlanetScope image mosaics (available across tropical forest areas 79 only [38]) allow for unprecedented monitoring of forest loss at frequent time intervals 80 and at minimum costs. While these satellite datasets are at a 3-4 m and 10 m resolution 81 respectively, commercial satellite data at higher resolutions (up to 0.3 m) are available 82 on request, and there remains a long archive of 30 m resolution Landsat imagery which, 83 if suitable for mapping degradation, would open up a long and consistent record of 84 forest disturbance. However, it is unclear which of these datasets, if any, is suitable for 85 mapping forest degradation. Hence, there is an increasing necessity of matching satellite 86 observations against ground truth data, showing the exact magnitude and location of 87 forest loss and thus reducing the critical uncertainties in biomass change estimates [39]. 88 The Normalized Difference Vegetation Index (NDVI) is one of the most widely 89 used indices for estimating tree aboveground biomass (AGB). NDVI is computed from 90 the ratio of the near-infrared (NIR) and red spectral bands. For healthy vegetation, the 91 presence of chlorophyll pigments causes a strong absorption in the red band, while 92 the internal leaf structure results in high reflectivity in the NIR. Therefore, healthy, 93 dense canopy cover tends to have a high NDVI value, while degradation of ecosystem 94 vegetation such that trunks, branches or the soil can be seen, or a decrease in the density 95 of green vegetation, would result in lower NDVI [40,41]. Unfortunately, NDVI saturates 96 in regions of high forest biomass, partially due to the low penetration capability of the 97 red band [42]. As multiple layers of leaves are not entirely visible to the sensor, forests 98 do not necessarily increase their canopy density as they become taller. In addition to this, 99 biomass estimation methods based only on spectral information are complicated by the 100 presence of canopy shadowing from large stands and the heterogeneity of forest stand 101 structures [43][44][45]. Vegetation indices have shown low to moderate agreements with 102 field data in temperate and boreal forests, where canopy cover is more sparse and there 103 is less species diversity [46,47]. However, in tropical and subtropical regions, with high 104 biomass levels and multi-storied forest canopies, and where there is a greater diversity 105 of tree species, vegetation indices have shown less sensitivity to canopy variations, with 106 low or insignificant results [33].

107
If changes in direct reflectance, or vegetation indices, are unlikely to show subtle 108 degradation, it is possible that changes in the texture of reflection or index patterns will 109 show after disturbance. The texture of images has shown potential for overcoming the 110 existing problems with biomass saturation [48], both in optical and SAR images [46,49-111 52], and thus it would follow that the same could be true for biomass change. Formally, 112 texture analysis is an image processing technique that measures the variability in pixel 113 intensity among neighboring pixels within a window of fixed size [48]. It contains 114 information on the structural and geometric properties of forest canopies [53]. However, 115 most of the previous texture-based biomass estimation methods have used medium 116 resolution images, such as the 30 m resolution Landsat TM data [54,55], with fewer 117 studies employing higher resolutions datasets [46,48]. Texture metrics derived from finer 118 spatial resolutions is expected to correlate better with field data, as the pixels will not 119 mix several large canopies (which likely have diameters in the 10-30 m range) [50,[56][57][58][59].

120
So far, texture-based biomass models have been developed to produce one-time maps 121 of forest biomass [48,54]. In order to determine forest gains and losses, two such maps 122 need to be subtracted from each other. The uncertainty on the final product is computed 123 from the propagation of the combined uncertainties of the two static maps, which means 124 that errors at the pixel level can be very large [60]. Considering that forest degradation 125 is a dynamic process, one way to reduce these uncertainties might involve deriving a 126 single regression model on the difference images. This approach requires experimental 127 data of biomass change (∆AGB) to derive the model, where the size and timeframe of 128 the disturbance is precisely known. quality of the best performing models for potential use in forest research and monitoring 138 applications. The South Eastern Peruvian Amazon, including the Madre de Dios region, has 142 been identified by the Peruvian government as a tropical 'Capital of Biodiversity' 1 .

143
The area is threatened by deforestation and degradation following the paving of the 144 tri-national interoceanic highway, opening up the area for development, and prolific 145 illegal goldmining operations [61]. Our study area is located in the North Eastern part  Scanning (TLS), providing digital three-dimensional models of every tree in the plot. 166 We also collected UAV laser scanning (LiDAR) and multispectral optical data (< 10 cm   stem. Wood density was derived from species information using a global wood density 187 database [66]. After the pre-logging measurements, Bélgica extracted a total of 24 trees, 188 in different proportions in each plot, to simulate a range of degradation events. For 189 the post-logging forest inventory, the direction of tree felling and the collateral damage 190 caused by the logging operations was also noted.

191
To convert forest inventory metrics into AGB we used the pantropical allometric 192 model from [67], using field-derived DBH (D) and wood density (ρ), and a regional Unfortunately, individual-tree biomass estimates from allometric equations are forms and biomass for large stems [69]. On the other hand, it was found that AGB 200 estimates from TLS measurements provide a better agreement with the reference data 201 from destructive harvest experiments, with a mean tree-scale relative error of 3% [70].

202
Considering that the logged trees in this study have an average DBH of 79 cm, and the number of extracted trees varies between 2 and 9 in each plot, we used TLS-derived AGB 204 for the logged stems, yielding likely a more accurate measurement of biomass change.

205
TLS data were collected using a Riegl VZ-400 scanner, following the protocols 206 outlined in [71]. Each plot was subdivided in squares of 10 x 10 m and two scans 207 were obtained at every grid intersection point. The point clouds of the logged trees 208 were extracted using the treeseg package [72]. From the quantitative structural models 209 (QSMs) [73,74], the AGB of each logged tree was derived by multiplying the model 210 volume by wood density. Total AGB loss in each plot was then obtained by summing 211 the TLS-derived AGB of the logged trees plus the change in AGB of the remaining trees 212 as measured from the field surveys (Table S1, Figure 2). Unfortunately, processed TLS 213 data volumes were not available for non-logged trees in the plots, so no comparison 214 between allometric and TLS biomass values could be made for these trees. We expect  (Table S2). 220 Figure 2. AGB loss in Mg ha −1 for the four 1-ha selectively logged plots. Forest loss amounts between 5% and 30% of initial biomass.

221
Commercial and open remote sensing data were collected over the field site before 222 and after the logging event (  calibrated products were then georeferenced using the processed WorldView-3 images.

270
As for the WorldView-3 data, the pre-logging and post-logging images were aligned to 271 match their extents.

272
The SPOT ("Satellite pour l'Observation de la Terre") series was developed within land surface class [86].

341
The image pre-processing pipeline is shown in Figure 3. In addition to the steps 342 described in the diagram, PlanetScope, Sentinel-2 and Landsat 8 were resampled to a 343 1 m resolution using a Nearest Neighbour convolution, to allow for easy and direct 344 comparison with the higher resolution data. 345 Figure 3. General methodology workflow used for producing biomass change estimates over the study area using fieldderived ∆AGB and remotely sensed data.

346
After preprocessing, the resulting red and NIR bands were used for the calculation 347 of the Normalized Difference Vegetation Index (NDVI): The choice of this particular index was motivated by its widespread use in land-349 cover change detection and monitoring of vegetation dynamics [87,88], and the fact 350 that it can be calculated from all the six satellite dataset types without the need of 351 sensor-specific calibration constants, making it suitable for multi-sensor comparison.

352
All the sensors had blue, green, red and NIR bands (BGRN), so these were also tested 353 individually, and stacked into a raster with the NDVI band for subsequent analysis. Texture is a function of local variance and is directly related to the resolution and 369 size of the dominant feature. In this study, this corresponds to the size and spacing of 370 the tree crowns. Low variance is expected in high resolution images where the pixel size 371 is smaller than the tree crown size, though it will rise around the edge of such crowns, 372 effectively acting as an edge detector. As pixel size increases to a scale that is comparable 373 to the tree crown, local variance increases too, and this should especially apply to 374 tropical forests, rich in tree species and with heterogeneous canopy structures [46]. A

407
A visual, multi-scale representation of some of the remotely sensed variables used 408 in this analysis is illustrated in Figure 5, with the RGB images showing plot C1 before 409 and after logging (Figure 5a). Figure 5b and Figure 5c show the difference in NDVI and 410 in the texture of the NIR band before and after logging for the same plot.  (Table S4).

412
Linear regression models were calculated using biomass change as the dependent 413 variable and all other parameters as the independent variables. For the spectral pa-414 rameters, the best fitting models are shown in Figure 6, with difference in simple NIR 415 reflectance illustrated in Figure 6a and difference in NDVI shown in Figure 6b. All and biomass change, is available in Table S3. 420 Figure 6. Multi-sensor comparison of regression models for (a) biomass loss and ∆NIR and (b) biomass loss and ∆NDVI. and for a window size of 3×3 pixels (R 2 = 0.97). The lowest performance for texture 425 measurement is observed in SkySat, with the best performing model obtained by the 426 texture of the NDVI and using a window size of 7×7 pixels (R 2 = 0.50). To illustrate how 427 spectral and textural-based models compare, Figure 7 shows the best-fitting models for 428 texture differences obtained from the NIR band ( Figure 7a) and NDVI (Figure 7b).

451
Biomass change maps were produced from PlanetScope, SPOT-7 and Sentinel-452 2 texture change data, as those were the optical sensors yielding the highest model 453 performance (Figure 8). Figure 9 compares the three biomass change maps of the 454 study area constructed from the best fitting models, where the textural parameter was 455 calculated from the NIR band on a window of 9×9 pixels for SPOT-7 (Fig. 9a), and 3×3 456 pixels for PlanetScope ( Figure 9b) and Sentinel-2 ( Figure 9c) and with a 1-ha resolution.

462
In this study we tested the ability of different high-medium resolution optical  as the texture measure to compute the variability in reflectance among neighbouring 495 pixels, to identify potential areas of abrupt change. We found that the best performing 496 texture-based models were derived from the NIR band (Figure 7a). This is due to the 497 higher penetration depth of this band compared to visible light [42]. value, and this is true for all the sensors analysed here ( Figure S1). This suggests that 519 the measurement of local variance in forest structures is enabled when crown sizes are 520 similar to the window resolution [46]. and hence the transferability of the model [55]. As a result, individual satellite-based 529 biomass maps over tropical areas normally have large pixel-level uncertainties [60,89].

530
Differencing such maps to produce a biomass change product means that substantial 531 spatial errors would be indistinguishable from the physical changes in forest cover. This 532 is the reason why differencing such maps is not recommendable, and there are currently 533 no trusted biomass change maps over the tropics.

534
In this study we used field measurements of biomass change where the entity and 535 the timescale of the disturbance events is well known, in order to reduce the uncertainty 536 caused by the temporal mismatch between satellite and field data collection. We used LiDAR data will improve the prediction and accuracy of our models. We hope that 556 the information presented here on the best combination of optical satellite and textural 557 parameter can be used as a template for future work on biomass change estimation.

559
The objective of this study was to link experimental data of biomass change to  Weaker relationships were observed for the very high resolution sensors (< 1.5 m).

581
Insignificant results were found for Landsat data at 30 m spatial resolution. This 582 suggests that local measurements of canopy variation are possible when the data is 583 neither too detailed (which increases the noise content) or coarse (which smooths 584 the information content);

585
• We found that there is an optimal window size for detecting forest degradation 586 when extracting textural information. This size slightly varies across the sensors; 587 in this study, the optimal range was found between 9 -14 m in length for a square 588 window, suggesting that local variation in forest canopies is best measured when 589 the window size is of comparable scale to the tree crown diameter.

590
Future research in the FODEX project will expand this analysis and attempt at re-591 ducing the noise content of the biomass maps by fusing multispectral optical information 592 with the remote sensing radar measurements acquired over the study area. Local LiDAR 593 data will also provide a larger sample size, improving the accuracy of the model. The

594
Acknowledgments: We would like to thank La Comunidad Nativa de Bélgica for allowing us 607 to conduct this research on their land, and for offering material and moral support during our 608 permanence in their community. We would like to specially thank the NGO AIDER -Asociación 609 para la Investigación y Desarrollo Integral, for offering their invaluable help with the campaign 610 preparation and for continuously supporting us with the logistics, ensuring the best outcome for 611 our fieldwork. This research would have not been possible without the work of the field assistants.

612
In particular, we would like to thank Roxana Sacatuma Cruz and José Sánchez Tintaya, and the 613 field assistants from Bélgica, Arturo Aspajo López, Leoncio Aspajo Lopez, Luis López Chapiama 614 and Kenny López Batista, for assisting us with the field inventory work, but also for their generous 615 guidance and advice. We also acknowledge the assistance provided with customs and logistics 616 by Prof. Eric Cosio and Dr. Norma Salinas from the Pontificia Universidad Católica del Perú. 617 We finally thank the Regional Government of Madre de Dios for issuing the Authorisation for