DIGITAL ELEVATION MODELS FOR VEGETATION MONITORING: A CASE STUDY OF A FOREST DISTRICT IN POLAND

Monitoring vegetation cover is one of the prime aims of remote sensing. Spatio-temporal trajectories of the vegetation state are key findings used by various branches of sciences for research and planning of human activities purposes. Among available vegetation monitoring methods, including forests, are multispectral and microwave indices. However, there are other ways to monitor vegetation/forest cover using surveys performed from satellite orbit. These approaches, including LiDAR and InSAR-derived digital elevation data products (further DEMs), are relatively less discussed in the literature on forest monitoring. Here, we present the results of a study of forest monitoring using the difference between two DEMs captured approximately 15 years apart. DEMs involved are the Shuttle RadarTopography Mission (SRTM) and Copernicus. The study area is located in a forestry district in the southwestern part of Poland. Spatial constraints and forest stand age, and dominant species were available from a forest stand map of the district. The study concludes that DEMs difference can be used for preliminary forest change assessment.


INTRODUCTION
The last two decades witnessed the introduction of global/quasiglobal digital elevation data products, commonly referred to as the digital elevation model (DEM). Although photogrammetry can be considered a mature method for delivering global DEMs, the Synthetic Aperture Radar Interferometry (InSAR) technology should be considered a "gamechanger" for global topography measurements. This is because of the SAR instrument's allweather/time operation capabilities. Therefore, it can be anticipated that the InSAR-derived DEMs with ever-rising accuracy and spatial resolution will be prevalent in the coming years and decades.
For the bare earth's topography survey, the SAR data possess a drawback associated with the microwaves' impenetrability of vegetation (Dall, 2007;Becek, 2011, Schlund et al., 2019, resulting in pixels having overestimated elevation. As it was found, the magnitude of the overestimation depends on two vegetation parameters, i.e., vegetation height and vegetation density (Becek, 2011). It was demonstrated that the InSAR elevation is overestimated by 60% to 90% of forest height for Cband and X-band derived DEM, respectively. Because of this phenomenon, a note on terminology used in this paper is required. First, DEM is a generic term used to describe any elevation dataset. Second, Digital Terrain Model (DTM) is a term used to describe a data set representing the bare earth's topography without any objects above the surface. Third, Digital Surface Model (DSM) is a term used to describe an elevation dataset representing a vertical outline of topography, including all objects above the surface (Maune, 2007). Since an InSARderived DEMs over a forest do not represent its canopy, it cannot be described as a DSM. Hence, to avoid potential confusion in this paper, we use the generic term DEM, meaning the elevation of the phase centre of microwaves reflected from the vertical vegetation profile.
The mentioned drawback of InSAR-derived DEMs for modelling the bare earth topography may be used for forest aboveground biomass assessment (Becek, 2011). However, the InSAR-based forest aboveground biomass estimation method possesses a few features that make it superior to any multi-or hyperspectralbased method. The major problem with the spectral methods is that the spectral indices, such as Normalised Difference Vegetation Index (NDVI), saturate at low biomass levels. Because of this fact, in the biomass-reach tropical forest, the aboveground biomass cannot be effectively studied. An additional adverse factor is a persistent cloud cover over tropical regions, making cloud-free optical images acquisition nearly impossible.
The InSAR-based approach for forest biomass or forest monitoring, in the more general case, has largely been overshadowed by the Light Detection and Ranging (Lidar) method, including aerial laser scanning (ALS) (e.g., Liu et al., 2018), terrestrial laser scanning (TLS) (e.g., Sheridan et al., 2015;Chen et al., 2019), space-borne missions such as ICESat (e.g., Tulski & Becek, 2021), and follow-on ICESat-2 mission (e.g., Sun et al., 2020). However, the ALS survey is not always a viable option because of prohibitive costs for many forest inventory projects. In addition, the space-borne lidar does not provide sufficient spatial resolution of the samples for local-scale projects.
This study explores the difference between two InSAR-derived DEMs to study vegetation change that occurred during the time lapse between DEMs' acquisition.
Specific objectives of the study include calculating the difference between two DEMs, i.e., the Shuttle Radar Topography Mission (SRTM) and the Copernicus. These DEMs were captured some 15 years apart. Hence, the difference must contain a signal related to an increase/decrease in forest height or density. Subsequently, the calculated difference is analysed using basic statistical techniques to identify relationships between variations of the difference and forest age or fraction of the dominant species.
A forestry district located in southern Poland has been selected as a test site. Some 80% of forest stands in the district are Spruce as a dominant (coniferous) species. The study's results indicate the potential applicability of this approach for the preliminary assessment of vegetation cover change under various climatic conditions. It might be speculated that the method can be used as a stand-alone or in combination with other methods of vegetation monitoring, including SAR data (Bouvet et al., 2018;Collins & Mitchard, 2015;Dubayah et al., 2017;Ruiz-Ramos et al., 2020;Verhelst et al., 2021;Watanabe et al., 2018;Ygorra et al., 2021).

Area of Interest (AOI)
The study area is located in southwestern Poland, in the Sudety mountains. The elevation of AOI varies between 289 m and 982 m a.m.s.l. The topography is steep, occasionally mild. The area size of the district is 13,011 ha. Figure 1 shows the location of the forestry district Bystrzyca Klodzka. The Lat/Lon coordinates of the AOI are (Bottom Left; Top Right): 50°12'47", 16°26'04"; 50°25'20", 16°47'13".

Figure 1.
Location of the study area and forest stands (black outline). The area size of AOI is 335.1 km 2 , out of which approximately 76.15 km 2 is covered by forest. The forested area is divided into forest stands characterised by the type of dominant species, age and percentage of the dominant species. There are 1800 stands in the district. The forest stand is the minimum/maximum area size is 1 ha and 33.5 ha, respectively. Figure 1 also shows forest stands boundaries. The types of dominant species and their participation in the district are shown in Table 1. Figure 2 shows two Landsat images of a part of AOI taken in 2000 (left image) and 2020. In the left image, patches free of trees are visible. A careful photointerpretation of the images allows us to conclude that the patches resulted from timber harvesting, followed by reforestation.

Data
In this study, SRTM and Copernicus DEMs were used. Both DEMs were developed using the InSAR method. However, the SRTM instrument used C (5.6 cm) microwave band, while Copernicus was developed from data collected by the X-band (3.1 cm) instrument.  Table 4. Basic summary of the DEMs used in this project.

Method
The InSAR-derived DEM elevation over the forest is overestimated by 60% to 90% of forest height, depending on the microwave band used. The overestimation also is linearly dependent on the forest density. Selective logging of a forest stand decreases the overestimation to zero when all trees are removed from the stand. Figure 3 shows the relationship between different elevation readings over the area covered by forest. In the case of the InSAR-derived DEMs, available is the phase centre elevation only. The phase center elevation is located always above the terrain and below the forest canopy. Note also that the X-band elevation is always higher, but still below the lidar (not shown here) and DSM elevation. The elevation difference of DEMs captured at various times over a forest contains signals related to forest height and density changes, species type, and difference in penetration depth due to the wavelength used to create DEMs (e.g., Becek, 2011;Schlund et al., 2019). Hence, the difference between DEMs over the forested area can be expressed as per Equation (1): where DEM1/2 = older/never DEM, respectively e = forest height difference g = microwave's penetration depth difference ε = random error.
The random error -ε, is noise due to imperfection of instrument and measurement procedure and environmental factors. These errors are random, and their statistical properties, e.g., mean value and standard deviation, can only be estimated based on many measurements in similar conditions. This is the classic Gaussian error model. However, besides random errors, errors due to the slope of the terrain and pixel size impact the accuracy of the DEMs difference. According to (Becek, 2008), the variance of these error sources can be estimated using Equation (2) as follows: where p = pixel size s = slope at a given pixel.
Consequently, the variance of the error of the DEMs difference can be expressed as follows using Equation (3): where / 2 = instrument-induced error of SRTM and Copernicus modell, respectively p = pixel size s = slope at a given pixel.
Equation (3) proves that the accuracy of the DEMs difference, and any DEM, is not constant. Instead, it varies from pixel to pixel according to the slope at each pixel. It further means that DEMs difference is less accurate (or less sensitive), the larger pixel or steeper terrain (larger slope).
The DEM difference is a random process because of the ε term in Equation (1). As demonstrated (e.g., Becek et al., 2021), the DEMs' probability density function (pdf) follows the Laplace pdf, which for random variable d is expressed by Equation (4) (Norton, 1984): where m = location parameter l = scale parameter.
According to (Norton 1984), the maximum likelihood estimator of m is the sample median. Note that m corresponds to the mean value in the case of the normal distribution. The estimator of l is the mean absolute deviation from m as expressed by Equation (5) (Norton, 1984): where n = number of pixels in DEM m = location parameter.
The variance of the DEM difference can be calculated from Equation (6) where l = scale parameter.
The parameters m and σ were used instead of the mean and standard deviation.
The following data processing workflow has been adopted in this study:  For the Spruce dominant stands, mean DEMs difference and its standard deviation vs a fraction of dominant species have been calculated to assess the impact of the leaf-on/off status of deciduous species.  A histogram of DEMs difference has been prepared, and Laplace's probability density function (pdf) has been fitted.  The histogram has been analysed to identify timber extraction activities or some natural causes.  Identified stands have been assessed using photointerpretation of high-resolution aerial images (acquired in 2010) and recent Sentinel-2 imagery.

RESULTS
An initial step of evaluating the suitability of the DEMs difference for forest change assessment was to identify the impact that leaf-on/off status of trees and microwaves penetration depth (C-vs X band) has on the DEMs difference. Figure 4 shows the median DEMs difference for each dominant class vs the fraction of dominant species. As it can be seen, a downward trend is present, suggesting the data used for the Copernicus model were predominantly collected during the leaf-on period. This situation is due to the large number of scatterers (leaves) in the canopy cover that move the phase centre of SAR pixels upward during the leaf-on period. Annotations A through C pinpoints some discrepancies between the histogram and the Laplace pdf that might suggest changes in forest stands. In the next section, this issue will be further discussed. The location and scale parameters of the Laplace pdf are also marked in the figure.

DISCUSSION
The DEMs difference for forests stands vs the fraction of the dominant species shown in Figure 4 reveals that Copernicus exhibits higher impenetrability than SRTM. The mean value of the median difference for all fractions of the dominant species is 1.39 m. A downward trend suggests that the impenetrabilities of the C-and X-band are less pronounced for a higher fraction of coniferous species. This effect is because the impenetrability of the C-and X-band DEMs over coniferous forests remains unchanged throughout the year. Hence, one can estimate that the impenetrability over Spruce forest for Copernicus is approximately 1 m higher than that of SRTM. This conclusion is essential for deciding if the forest has changed. On the other hand, higher DEMs differences for lower fractions of the dominant species suggest that the Copernicus was developed using a mixture of data taken at various seasons (during leaf-on and off status).
A histogram of the DEMs differences shown in Figure 3 follows the Laplace probability density function (pdf). The pdf was drawn using the location and scale parameters estimated from the differences. This is proof that the DEMs difference follows the Laplace pdf instead of the Gauss one, which some authors attempt to use to model the histogram of DEMs differences.
The discrepancies between the Laplace pdf and the histogram marked as A in Figure 5 suggest some numbers of clear cuts among the investigated forest stands. However, visual inspection using recent satellite imagery suspected stands did not confirm this assertion. Therefore, the lack of visual confirmation might be caused by the reforestation of the harvested plots. On the other hand, the discrepancies between the histogram and pdf curve marked as B in Figure 5 are probably due to the forest height increase of stands dominated by younger trees.
Inconsistencies marked as C are probably due to the forest stands' uneven age structure: older trees gain height significantly slower.
Overall histogram of the DEMs difference does not provide clear evidence of apparent forest change between 2000 and 2015.
Knowing, however, that timber extraction in the forest district is routinely carried out, it can be concluded that tree harvesting between 2000 and 2015 was mainly absent.
A more suitable way to identify forest change using the DEMs difference offers a scatter plot in Figure 6. The trend line suggests that forest stands with younger trees increase their height (positive DEMs difference), while stands with older trees tend to exhibit negative difference, suggesting removal of trees.

CONCLUSIONS
This preliminary study uses the elevation difference between two DEMs captured some time apart, demonstrating the utility of the information content of freely available DEMs for monitoring vegetation cover at a local to global scale.
One of the critical factors controlling the sensitivity of the DEMs difference for monitoring vegetation and changes in topography is the slope, which is expressed by Equation (3). This study was carried out in relatively steep topography, implying a low sensitivity of the method (lower accuracy of the DEMs difference). However, it could still identify some natural and anthropogenic changes occurring in the forest district.
The presented method can provide the first evidence of change in vegetation cover. It is anticipated that further research in this direction might identify some quantitative indicators, both spatial and aspatial, of vegetation cover change.