Analysis of Using Dense Image Matching Techniques to Study the Proce s of Secondary Su ce sion in Non-Forest Natura 2 0 Habitats

Secondary succession is considered a threat to non-forest Natura 2000 habitats. Currently available data and techniques such as airborne laser scanning (ALS) data processing can be used to study this process. Thanks to these techniques, information about the spatial extent and the height of research objects—trees and shrubs—can be obtained. However, only archival aerial photographs can be used to conduct analyses of the stage of succession process that took place in the 1960s or 1970s. On their basis, the extent of trees and shrubs can be determined using photointerpretation, but height information requires stereoscopic measurements. State-of-the-art dense image matching (DIM) algorithms provide the ability to automate this process and create digital surface models (DSMs) that are much more detailed than ones obtained using image matching techniques developed a dozen years ago. This research was part of the HabitARS project on the Ostoja Olsztyńsko-Mirowska Natura 2000 protected site (PLH240015). The source data included archival aerial photographs (analogue and digital) acquired from various phenological periods from 1971–2015, ALS data from 2016, and data from botanical campaigns. First, using the DIM algorithms, point clouds were generated and converted to DSMs. Heights interpolated from the DSMs were compared with stereoscopic measurements (1971–2012) and ALS data (2016). Then, the effectiveness of tree and shrub detection was analysed, considering the relationship between the date and the parameters of aerial images acquisition and DIM effects. The results showed that DIM can be used successfully in tree and shrub detection and monitoring, but the source images must meet certain conditions related to their quality. Based on the extensive material analysed, the detection of small trees and shrubs in aerial photographs must have a scale greater than 1:13,000 or a 25 cm GSD (Ground Sample Distance) at most, an image acquisition date from June–September (the period of full foliage in Poland), and good radiometric quality.


Introduction
Secondary succession is the process of re-establishing the original community after a disturbance [1].Currently, three-quarters of Earth's ice-free terrestrial biomes are being disturbed [2], but at the same time, many degraded areas are in a natural or/and anthropogenic recovery state [3].Secondary succession is consequently becoming the primary focus of terrestrial natural resource

Introduction
Secondary succession is the process of re-es disturbance [1].Currently, three-quarters of Earth's i [2], but at the same time, many degraded areas are in [3].Secondary succession is consequently becoming th

Introduction
Secondary succession is the process of re-establishing the original community after a disturbance [1].Currently, three-quarters of Earth's ice-free terrestrial biomes are being disturbed [2], but at the same time, many degraded areas are in a natural or/and anthropogenic recovery state [3].Secondary succession is consequently becoming the primary focus of terrestrial natural resource management for projected global land use changes [4].Generally, investigations of succession dynamics have been a central theme in plant community ecology [5][6][7][8].
Open areas together with their non-forest vegetation (especially less productive areas) that are currently abandoned result in secondary succession [9][10][11].Cessation of mowing or grazing causes species with clonal growth to complete their full development and induce changes in the quantitative and spatial structure of plant communities [12].The results of this process are the disappearance of some species groups (e.g., heliophilous) and the formation of shrub and forest communities created by species more adapted to poor light conditions.As a result, secondary succession also leads to changes on the landscape level.The abandonment of agricultural meadow use is one of the main reasons for the deterioration of biodiversity and conservation status [13,14].
In the past, a field inventory was the only reliable way to obtain information about tree parameters.However, in large or difficult-to-reach areas, the process is costly, and more importantly, of a discrete character.The development of remote sensing techniques allows for much faster measurements of some variables and is performed on the whole area of interest [15].Forest monitoring (quantitative and qualitative) is frequently based on multispectral [16] and hyperspectral imagery [17], as well as LIDAR (Light Detection and Ranging) [18] and SAR (Synthetic Aperture Radar) [19,20].
One of the most important and basic features observed is tree height.It is of great importance since it allows for tree segmentation [21], vegetation succession assessment [22], and biomass estimation [23].It can be obtained in the field by forest management planners; however, these measurements are very time consuming and are performed mostly on inventory plots [24].That is why an automatic method for extracting 3D information is needed.For the last twenty years, this problem has been solved mostly by using LIDAR data [25,26].This results in a point cloud that allows for a precise description of a forest stand.One of the most significant features of airborne laser scanning (ALS) is that the laser beam (LIDAR technology) can penetrate the canopy, enabling measurements under the treetops.This feature allows for simultaneously obtaining a digital surface model (DSM) and digital terrain model (DTM), thus allowing for a canopy height model (CHM) generation named, also as a normalized digital surface model (nDSM).However, LIDAR is a relatively young technology widely used for no more than twenty years in civil applications.Due to this fact, most archival remote sensing data are in the form of aerial images.Therefore, historical monitoring of past forest states using LIDAR is practically impossible [27].
Meanwhile, image processing has become more accessible.The decreasing cost of software, paired with the development in hardware and computing power, has caused a sudden spike in the popularity of dense image matching (DIM), as a cheaper source for independent 3D information.Several studies confirm that proper image matching techniques allow for a generation of DSM models with similar or better accuracy and much larger density than those extracted from a LIDAR point cloud [28][29][30].Semiglobal matching is especially regarded as a very accurate matching method and is performed in many vegetation studies [31][32][33].Creating a CHM from an image point cloud requires using a DTM from another source since stereophotogrammetry does not provide height information under the canopy.Usually, the DTM from a previously acquired ALS point cloud is used.It is widely assumed that the ground topography remains unchanged over the years [27,34].
What is more, in some cases DIM is the only technology that can be successfully used for automated detection of changes in the growth of forest over a long period as it requires the use of archival data collected mostly as analogue photographs.Such images, after proper scanning and processing, can be used successfully in retrospective vegetation analysis.However, much depends on the quality of the input data, mostly the scale of the image and its texture.The use of archival data is also limited by the lack of multispectral information and poor signal-to-noise ratio [35].Panchromatic images cannot provide multispectral information for vegetation indices calculations, but this disadvantage does not influence the possibility of using DIM techniques because they are based on a selected or simulated (average) image channel.The use of multitemporal images may cause some difficulties due to different image scales and qualities, different vegetation growth states, and different lighting conditions during the flight.Thus, it is recommended to perform flights on similar dates, with similar atmospheric conditions and similar image acquisition parameters [29].In some cases, the lack of proper ground control points for older images is also an important factor that constricts adequate aerotriangulation [34].The results may be improved by transforming scanned photos to a form similar to the images from a digital camera [36].
Image point clouds and derivative products are much better at the detection of major changes, such as clear cuts, than at smaller changes such as the removal of individual trees from the nondominant canopy layer [37].Similar assumptions can be made for the assessment of secondary succession (growth of trees and shrubs species).Smaller plants and stand-alone trees can be omitted in the CHM from an image point cloud, especially from lower quality data.
Secondary succession can be assessed in many ways using remote sensing techniques.The first attempts were based on photointerpretation with the added use of stereoscopic observations [15].In other research, Oikonomakis and Ganatsas [38] performed on-screen classification of succession areas based on the visual interpretation of land cover using orthophoto from 1945-2009.One class of openings and three classes of forest with varied density were created.Finally, a map of new forested areas was obtained.However, the labour intensity led to the development of more automated methods.
Maximum likelihood classification was used in Hernández et al. [39].Native forest, arborescent shrubland, dense espinal, good espinal, very degraded espinal, agriculture, water, urban and barren land, and plantation land cover classes were identified.Landsat satellite data allowed for analysis over a long period of time .Transition matrices were used to examine changes in different periods.Batistella and Lu [40] classified 31 sample plots into three secondary succession stage classes and forest, based on field collected vegetation structure variables and Landsat TM imagery.The authors stated that part of the field work could be replaced with LIDAR statistics.However, using only satellite imagery required merging two classes due to the similarities in spectral response between initial and intermediate succession.Differentiating the succession stages of mangroves was examined by Aslan et al. [22].The canopy height was modelled based on ALOS PRISM data.A histogram analysis (focusing on skewness) of the resulting CHM was used to assess whether such data could be used for the described task.Szostak et al. [41] proposed using geographic object-based image analysis (GEOBIA) in monitoring secondary forest succession.Their study shows that applying GEOBIA to high-resolution multispectral imagery can produce similar results to traditional on-screen photointerpretation techniques.
There are also papers in which the use of LIDAR-derived metrics and statistical modelling to predict forest succession stages are proposed [42].Regular ALS campaigns or the use of alternative technologies, such as stereo-matching of aerial photographs or radar technologies, give a good chance to manage and monitor the changes in rural areas [43] and generate 3D spatial indices describing the vegetation structure [44].LIDAR-derived metrics provide stable and repeatable measures, confirming the suitability of this data source for vegetation monitoring [45].
The aim of this study was to use DIM techniques to examine and assess the dynamics of the secondary succession process on non-forest Natura 2000 habitats (using archival aerial photographs) and to analyse whether their effectiveness depends on the species of trees and shrubs.These techniques enable automatically obtaining information on the range and height of trees and shrubs, which allows for faster and more comprehensive assessments of the dynamics of this process.In addition, since archival aerial photographs are characterized by different geometric and radiometric quality, the second objective of the study was to determine the minimum parameters that archival photos must meet in order to effectively determine the range of trees and shrubs in the early stages of the succession process.

Study Area
The research was carried out on the Ostoja Olszty ńsko-Mirowska Natura 2000 protected site (PLH240015).The study area is in the southern part of Poland, in the Silesian Voivodeship, near Czestochowa city (50 • 45" N; 19 • 17" E).It encompasses an area of about 25 km 2 .The area includes a complex of limestone hills (absolute heights range between 278 and 374 m).The hills are covered by natural forest communities (mostly beech and hornbeam forests) or semi-natural species-rich grasslands [46].Areas adjacent to the hills are occupied by cultivated fields and pine forests.The area is characterized by a high biodiversity of habitats with a considerable wealth of plant and animal species [47], so it is protected as a Natura 2000 site (Figure 1).The most important threatened non-forest habitats in this area are two types of grassland (codes 6120 and 6210) [48].The first type, habitat 6120, consists of dry, frequently open grasslands on more or less calcareous sands with a subcontinental centre of distribution, comprising semi-natural, moderately open to closed, relatively low-grown meso-xeric grasslands on slightly calcareous sands.This habitat is dominated by tussock-forming, narrow-leaved grasses of the Festuca ovina aggregate, often accompanied by Agrostis capillaris, Poa angustifolia, or Carex praecox.The floristic diversity of the habitat 6210 grasslands is created by many rare species of vascular plants [49].
Remote Sens. 2019, 11, x FOR PEER REVIEW 4 of 27 complex of limestone hills (absolute heights range between 278 and 374 m).The hills are covered by natural forest communities (mostly beech and hornbeam forests) or semi-natural species-rich grasslands [46].Areas adjacent to the hills are occupied by cultivated fields and pine forests.The area is characterized by a high biodiversity of habitats with a considerable wealth of plant and animal species [47], so it is protected as a Natura 2000 site (Figure 1).The most important threatened non-forest habitats in this area are two types of grassland (codes 6120 and 6210) [48].The first type, habitat 6120, consists of dry, frequently open grasslands on more or less calcareous sands with a subcontinental centre of distribution, comprising semi-natural, moderately open to closed, relatively low-grown meso-xeric grasslands on slightly calcareous sands.This habitat is dominated by tussock-forming, narrow-leaved grasses of the Festuca ovina aggregate, often accompanied by Agrostis capillaris, Poa angustifolia, or Carex praecox.The floristic diversity of the habitat 6210 grasslands is created by many rare species of vascular plants [49].One of the main threats to maintaining the proper level of species and protecting the semi-natural habitats, which are the subjects of conservation in the above-mentioned Natura 2000 site, is the cessation of pastoralism, and lack of grazing and mowing, that is the direct cause of secondary succession on grasslands [48].This process has been observed since the 1990s when there was a political change in Poland and agricultural activity in this area gradually lost its importance.Therefore, there is a need to monitor the succession process to plan effective activities within the framework of active protection of Natura 2000 habitats, such as grubbing shrubs and trees, and controlling sheep grazing, which were implemented in the research area, among others, from 2012-2016, as part of the LIFE11 NAT/PL/432 project "Protection of valuable natural non-forest habitats of the Eagle's Nest Landscape Park" [50].The main taxa appearing as a result of secondary succession in this area are Pinus sylvestris, Betula pendula, Prunus spinosa, Juniperus communis, Rhamnus cathartica, Crataegus spp., Cornus sanguinea and Rosa spp.(Figure 2).One of the main threats to maintaining the proper level of species and protecting the semi-natural habitats, which are the subjects of conservation in the above-mentioned Natura 2000 site, is the cessation of pastoralism, and lack of grazing and mowing, that is the direct cause of secondary succession on grasslands [48].This process has been observed since the 1990s when there was a political change in Poland and agricultural activity in this area gradually lost its importance.Therefore, there is a need to monitor the succession process to plan effective activities within the framework of active protection of Natura 2000 habitats, such as grubbing shrubs and trees, and controlling sheep grazing, which were implemented in the research area, among others, from 2012-2016, as part of the LIFE11 NAT/PL/432 project "Protection of valuable natural non-forest habitats of the Eagle's Nest Landscape Park" [50].The main taxa appearing as a result of secondary succession in this area are Pinus sylvestris, Betula pendula, Prunus spinosa, Juniperus communis, Rhamnus cathartica, Crataegus spp., Cornus sanguinea and Rosa spp.(Figure 2).

Data
Aerial photographs from 1971, 1982, 1996, 2003, 2009, 2012, and 2015 were obtained from the State Geodetic and Cartographic Resource for the study area.They were both analogue and digital, black and white (panchromatic), and colour (RGB).They differed in parameters, with various scales or pixel sizes (scale: from 1:13,000 to 1:32,000 scanned with 1814 dpi (14 µm) or 2400 dpi, pixel size: 24 and 25 cm), and acquisition dates (leafless period-spring/autumn/winter or the period of vegetation development-summer).The basic parameters characterizing these photos are presented in Table 1, and Figure 3 presents fragments of photos for the same area.In the case of aerial photos from 1971 and 1982 scratches are visible on the film for the photos from 1971.The photos from 1971 and 1982 also have visibly worse radiometric quality, which may negatively affect the efficiency of DIM algorithms which was also tested in the experiment.

Data
Aerial photographs from 1971, 1982, 1996, 2003, 2009, 2012, and 2015 were obtained from the State Geodetic and Cartographic Resource for the study area.They were both analogue and digital, black and white (panchromatic), and colour (RGB).They differed in parameters, with various scales or pixel sizes (scale: from 1:13,000 to 1:32,000 scanned with 1814 dpi (14 µm) or 2400 dpi, pixel size: 24 and 25 cm), and acquisition dates (leafless period-spring/autumn/winter or the period of vegetation development-summer).The basic parameters characterizing these photos are presented in Table 1, and Figure 3 presents fragments of photos for the same area.In the case of aerial photos from 1971 and 1982 scratches are visible on the film for the photos from 1971.The photos from 1971 and 1982 also have visibly worse radiometric quality, which may negatively affect the efficiency of DIM algorithms which was also tested in the experiment.In addition, ALS data from June 2016 were used in the study.The data were acquired using the Airborne Laser Scanner (FWF) Riegl LMS-Q680i (wavelength 1.55 µm; FOV max: 60 degrees; overlap: 62.7%; overlap width: 855 m).The point cloud was characterized by a density of 7 points per square metre.
The second source of data for the study was direct data from botanical field campaigns carried out in 2016-2017.During the botanical campaigns, information about 1000 trees and shrubs was obtained (Figure 4).The location of individual trees and shrubs was determined using a GNSS Mobile Mapper 120 with the real-time differentially corrected and a measurement accuracy of up to 0.2 m.For each tree or shrub, botanists determined the following parameters: species, height, diameter, and density of the crown.The selection of measuring points was made in accordance with established criteria.Only isolated individuals were measured, representing various species of trees and shrubs.Moreover, the full spectrum of height, size, and density of tree and shrub crowns was covered.In addition, since some individuals have grown since the 1960s, they could be used to assess the accuracy of the interpolated heights obtained from the DIM techniques in images from all dates.In addition, ALS data from June 2016 were used in the study.The data were acquired using the Airborne Laser Scanner (FWF) Riegl LMS-Q680i (wavelength 1.55 µm; FOV max: 60 degrees; overlap: 62.7%; overlap width: 855 m).The point cloud was characterized by a density of 7 points per square metre.
The second source of data for the study was direct data from botanical field campaigns carried out in 2016-2017.During the botanical campaigns, information about 1000 trees and shrubs was obtained (Figure 4).The location of individual trees and shrubs was determined using a GNSS Mobile Mapper 120 with the real-time differentially corrected and a measurement accuracy of up to 0.2 m.For each tree or shrub, botanists determined the following parameters: species, height, diameter, and density of the crown.The selection of measuring points was made in accordance with established criteria.Only isolated individuals were measured, representing various species of trees and shrubs.Moreover, the full spectrum of height, size, and density of tree and shrub crowns was covered.In addition, since some individuals have grown since the 1960s, they could be used to assess the accuracy of the interpolated heights obtained from the DIM techniques in images from all dates.

Methods
The workflow of the study is presented in Figure 5.The methodology consists of the following steps: 1. Interior and exterior orientation of archival aerial photos within a bundle adjustment (aerotriangulation) based on control points, 2. testing the impact of parameter selection on the quality of generated point clouds, 3. point cloud generation using DIM techniques, 4. DSM creation, 5. accuracy assessment of DSMs obtained using DIM techniques compared to stereo and field measurements, 6. comparison of CHM heights interpolated from the DSM (2015) and the LIDAR data (2016), 7. effectiveness analysis of tree and shrub detection considering the relationship between the date and parameters of aerial photos acquisition and the species and parameters of trees and shrubs.

Methods
The workflow of the study is presented in Figure 5.The methodology consists of the following steps: 1.
Interior and exterior orientation of archival aerial photos within a bundle adjustment (aerotriangulation) based on control points, 2.
testing the impact of parameter selection on the quality of generated point clouds, 3.
point cloud generation using DIM techniques, 4.
accuracy assessment of DSMs obtained using DIM techniques compared to stereo and field measurements, 6.
comparison of CHM heights interpolated from the DSM (2015) and the LIDAR data (2016), 7.
effectiveness analysis of tree and shrub detection considering the relationship between the date and parameters of aerial photos acquisition and the species and parameters of trees and shrubs.Professional photogrammetric systems, such as Trimble Inpho, require using approximate values for the exterior orientation parameters (EO) of photos in the process of aerotriangulation.Such data are rarely accessible for archival datasets, so photos were first processed in Agisoft software, which does not require as much input.The process of aerotriangulation in Agisoft is based on structure-from-motion (SfM) algorithms that are not dedicated for scanned analogue photograph processing and cause some errors in tie point extraction.To avoid these negative effects, the photos were pre-processed in SAPC software according to the methodology developed by Salach [36].This methodology transforms scans to a form similar to the images captured by a digital camera.Output images possess the same principal point position and the same resolution through cutting out the black photo frames.EO parameters obtained using Agisoft software were then used as approximate EO values in Trimble Inpho, which allowed the process of image block bundle adjustment to occur.Ground control points (GCP), necessary for the process of aerotriangulation, were calculated based on the oriented block of photos from 2009 (processed by data provider) as manual pass points.Their coordinates were calculated in bundle adjustment.Total accuracy of those points was assessed at around 50 cm.Roughly 15 points were used for each dataset.The number varied slightly depending on the visibility of some points on the photos.GCPs were divided into two groups: control points and checkpoints in a ratio of 2 to 1.
Cloud generation was performed based on images with known, calculated, and checked exterior orientations.For this purpose, DIM in Trimble Inpho was performed using feature based matching and least square matching algorithms in the Match-T DSM module.The overall assumptions for generating the point cloud required the creation of a relatively dense cloud of points with a low level of smoothing and a high permissible threshold of parallax in order to precisely capture individual objects, such as trees.Thus, many variants with different parameters were tested to achieve the best result (Table 2).Professional photogrammetric systems, such as Trimble Inpho, require using approximate values for the exterior orientation parameters (EO) of photos in the process of aerotriangulation.Such data are rarely accessible for archival datasets, so photos were first processed in Agisoft software, which does not require as much input.The process of aerotriangulation in Agisoft is based on structure-from-motion (SfM) algorithms that are not dedicated for scanned analogue photograph processing and cause some errors in tie point extraction.To avoid these negative effects, the photos were pre-processed in SAPC software according to the methodology developed by Salach [36].This methodology transforms scans to a form similar to the images captured by a digital camera.Output images possess the same principal point position and the same resolution through cutting out the black photo frames.EO parameters obtained using Agisoft software were then used as approximate EO values in Trimble Inpho, which allowed the process of image block bundle adjustment to occur.Ground control points (GCP), necessary for the process of aerotriangulation, were calculated based on the oriented block of photos from 2009 (processed by data provider) as manual pass points.Their coordinates were calculated in bundle adjustment.Total accuracy of those points was assessed at around 50 cm.Roughly 15 points were used for each dataset.The number varied slightly depending on the visibility of some points on the photos.GCPs were divided into two groups: control points and checkpoints in a ratio of 2 to 1.
Cloud generation was performed based on images with known, calculated, and checked exterior orientations.For this purpose, DIM in Trimble Inpho was performed using feature based matching and least square matching algorithms in the Match-T DSM module.The overall assumptions for generating the point cloud required the creation of a relatively dense cloud of points with a low level of smoothing and a high permissible threshold of parallax in order to precisely capture individual objects, such as trees.Thus, many variants with different parameters were tested to achieve the best result (Table 2).The resulting point cloud from each variant for each dataset was used for generating the DSM.This operation was performed in ArcGIS software with the parameters presented in Table 3.The interpolation of height in an elevation model from a few points to a model cell is particularly noteworthy, when the DSM grid is calculated.There are many publications on this topic related to ALS technology.Regarding equally dense image matching data, it should be noted that using the interpolation with maximum height selection is appropriate.
Finally, to assess the usefulness of the DSMs created for secondary succession, a comparison of heights obtained from DSM and stereoscopic observations of several dozen trees were made.Such measurements were collected for trees which have grown in the area of study since the 1960s and for which detailed information was obtained during the field campaign.Then, the heights measured as a result of stereoscopic observations were compared with the heights obtained automatically to find the best variant of the algorithm for the DIM.For each term analysed, the height difference, ∆, was calculated between the height from stereoscopic measurements and that obtained from the DIM algorithm.
Further analysis included the analysis of effectiveness for tree and shrub detection with respect to the date and the parameters of aerial image acquisition used in DIM to provide the conclusion about their use and an accuracy assessment of the DSM.Moreover, the comparison of heights interpolated from the DSM (2015) and LIDAR data (2016) on 1000 trees and bushes was analysed with respect to species division to evaluate the DIM technique in automatic detection of trees and shrubs in a DSM.

Results
In this section, the results and influence of dense image matching on point cloud generation, the DSM, and detection of vegetation were discussed.

Analysis of the Influence of Parameters during Image Matching on the Quality of Point Clouds
The impact of selected control parameters was analysed by comparing the following parameters: the number of filtered 3D points, time per generated DSM post, RMS (Root Mean Square) difference (Z) for the control used and adjusted points, estimated internal height accuracy (DSM), and elapsed time for DSM generation.The vertical accuracy did not change significantly regardless of the selected variant.A very large difference was visible in the density of the cloud, where variant w3 had, on average, 9 times more points than w1, w2, or w5, providing a point cloud density theoretically equal to one point per GSD value compared to other variants with image matching every 2, 3, or 5 pixels.This, of course, had an impact on the time needed to generate the cloud, but this was not a large enough effect to consider important.In this research, it was more important to be able to obtain a higher DSM accuracy.Additionally, in order to compare the impact of selecting the cloud parameters on the resulting DSM simply, images generated from different point clouds were subtracted from each other.This way, images of the differences between generated products were obtained.
The analysis of these images (Figure 6) revealed that small differences in the selection of parameters result in significant changes in the resulting cloud obtained and thus in the final rasters.
of selecting the cloud parameters on the resulting DSM simply, images generated from different point clouds were subtracted from each other.This way, images of the differences between generated products were obtained.
The analysis of these images (Figure 6) revealed that small differences in the selection of parameters result in significant changes in the resulting cloud obtained and thus in the final rasters.
Changing only the value of parallax (pair w2-w1) resulted mainly in point differences at the boundaries of high objects due to the mapping of an object or lack thereof.Changing the level of smoothing (w5-w1) and cloud density (w3-w1) affects the contours of objects, which differ as a result of smoothing.This effect is caused by stronger interpolation, while in case of cloud condensing, objects grow their perimeters.This increases the probability of tree detection as the cloud is denser and the crown has a larger area.This visual analysis confirms the choice of w3 as the best option when generating a cloud from images in Trimble Inpho software with low smoothing and the densest point cloud generated for points found for each corresponding pixel.

The Accuracy Assessment of the DSMs Obtained using Dense Image Matching Technique
Table 4 presents the accuracy assessment of the heights obtained based on archival photographs using dense image matching techniques versus stereo measurements.The visual presentation of these results is shown in Figure 7.The highest accuracy was obtained with the DSM from 2015, which was created from digital colour images of the best radiometric and geometric quality (GSD = 25 cm).These aerial photos were collected during the fully developed growing season; therefore, all species of trees and shrubs were visualised correctly.The standard deviation of the height difference Changing only the value of parallax (pair w2-w1) resulted mainly in point differences at the boundaries of high objects due to the mapping of an object or lack thereof.Changing the level of smoothing (w5-w1) and cloud density (w3-w1) affects the contours of objects, which differ as a result of smoothing.This effect is caused by stronger interpolation, while in case of cloud condensing, objects grow their perimeters.This increases the probability of tree detection as the cloud is denser and the crown has a larger area.This visual analysis confirms the choice of w3 as the best option when generating a cloud from images in Trimble Inpho software with low smoothing and the densest point cloud generated for points found for each corresponding pixel.

The Accuracy Assessment of the DSMs Obtained using Dense Image Matching Technique
Table 4 presents the accuracy assessment of the heights obtained based on archival photographs using dense image matching techniques versus stereo measurements.The visual presentation of these results is shown in Figure 7.The highest accuracy was obtained with the DSM from 2015, which was created from digital colour images of the best radiometric and geometric quality (GSD = 25 cm).These aerial photos were collected during the fully developed growing season; therefore, all species of trees and shrubs were visualised correctly.The standard deviation of the height difference determined by both methods equals 1.27 m (Table 4).Analysing the differences obtained for particular species, it was noticed that conifer species are relatively small (generally 0.2-0.4m, standard deviation of the height difference determined by both methods is 0.19 m), while for deciduous species the variation is much higher (standard deviation of the difference in height determined by both methods is 1.45 m, 0.54 after eliminating Populus balsamifera with 10% crown density).The largest differences in height were observed for individuals with open canopy closure (10%-20% of height).Analysing the data from 2015, it can be noticed that the size and the closure of its crown have a notable impact on the accuracy of the tree heights obtained using the DIM algorithm.In the case of less dense crowns, the contrast between the crown and its surroundings significantly affects the operation of DIM algorithms.
The least accurate results were obtained based on aerial photographs from 2012, 2009, 1996, and 1982.In 2012 and 2009, this was caused by the image acquisition date (March and April), which was at the beginning of the growing season in Poland when most of the trees and shrubs were still leafless or in the early stages of development.This issue limited the effective detection of trees and shrubs.Most values of ∆ are positive, which means that these trees in DSM were not contained with their height.Only in the case of coniferous trees was high accuracy obtained for most individuals (problems only occur in the case of pines with irregular crowns).The R 2 coefficient for coniferous trees was for 1971: 0.51, 1996: 0.36, 2003: 0.89, 2009: 0.64, 2012: 0.93, and 2015: 0.99.
Aerial photographs from 1996 and 1982 were made in smaller scales, 1:26,000 and 1:32,000, respectively, which lowered the accuracy of the DSMs generated, for both deciduous and coniferous species.The standard deviation of the height differences is 3.57 m for 1982 and 3.88 m for 1996 (Table 4).In both these terms, significant outliers (up to 13 m) can be observed for a few individuals.This applies particularly to deciduous trees (the biggest errors were obtained for Acer spp.).For photographs from 1971 and 1982, the poor radiometric quality and visible film damage (Figure 8) caused significant errors in the elevation model.
Summing up the performed accuracy analysis, the DSMs with the highest accuracy were obtained using images acquired in the fully developed growing season, with a scale of at least 1:13,000 (analogue cameras) or with a GSD less than or equal to 25 cm (digital cameras).

Analysis of the Effectiveness of the Tree and Shrub Detection using DIM
In order to evaluate the effectiveness of the DIM technique with respect to its use in the analysis of the succession process dynamics, an experiment with multitemporal data, various acquisition terms, and varied image qualities was carried out to determine what features of the examined objects affect the results of their detection.These analyses were conducted for nDSMs obtained using aerial photos from 2015 because the date of their acquisition was the nearest to the date of the field botanical campaigns (2016-2017) and the registration date of ALS data (2016).In addition, aerial photos from 2015 were acquired during the growing season, with GSD = 25 cm, therefore producing the best parameters from the analysed set of archival data.
The effectiveness analysis included an impact analysis of the parameters of trees and shrubs (height, crown diameter, crown density) on the tree and shrub detection using the DIM technique.The comparison of the height of trees and shrubs obtained based on ALS data and as a result of DIM techniques depending on the species, height, size, and crown density of the examined trees and shrubs are presented in Tables 5, 6, 7, 8, and in the sub-figures in Figure 9.The influence of crown diameter and density (by species) on the accuracy of the height of trees and shrubs using DIM techniques are illustrated in Figure 10.
The height differences indicate that a higher accuracy for determining the height of trees and shrubs is acquired for conifers (Table 5, µ|ΔALS-DIM| = 0.43 m).In the case of deciduous trees and shrubs, the mean value of height differences is 0.73 m; it is strongly differentiated (from 0.37 m to 1.64 m, and µ|ΔALS-DIM| = 0.73 m) and depends on the species and size of the measured object.In general, it was found that tree species, especially those with loose crowns (Betula pendula, Populus tremula), were characterized by higher values of height differences (Figure 10).In the case of Fraxinus sp., Robinia pseudoacacia, and Quercus rubra, larger differences arise from the specifics of objects present in the study area-some of them were characterized by a small crown density (30%-50%) (Robinia pseudoacacia) or smaller crown diameter (Figure 10), moreover Fraxinus sp. and Quercus rubra did not occur in large numbers.
Shrubs are characterized by a much smaller error in determining the height than is the case for woody species (the exception here is Corylus avellana).This applies to both deciduous and coniferous shrubs (Juniperus communis).The average height differences of shrubs obtained based on LIDAR data and using DIM techniques are 0.49 m for deciduous shrubs and 0.34 m for coniferous, respectively.However, considering their lower heights compared to woody species, the relative

Analysis of the Effectiveness of the Tree and Shrub Detection using DIM
In order to evaluate the effectiveness of the DIM technique with respect to its use in the analysis of the succession process dynamics, an experiment with multitemporal data, various acquisition terms, and varied image qualities was carried out to determine what features of the examined objects affect the results of their detection.These analyses were conducted for nDSMs obtained using aerial photos from 2015 because the date of their acquisition was the nearest to the date of the field botanical campaigns (2016-2017) and the registration date of ALS data (2016).In addition, aerial photos from 2015 were acquired during the growing season, with GSD = 25 cm, therefore producing the best parameters from the analysed set of archival data.
The effectiveness analysis included an impact analysis of the parameters of trees and shrubs (height, crown diameter, crown density) on the tree and shrub detection using the DIM technique.The comparison of the height of trees and shrubs obtained based on ALS data and as a result of DIM techniques depending on the species, height, size, and crown density of the examined trees and shrubs are presented in Tables 5-8, and in the sub-figures in Figure 9.The influence of crown diameter and density (by species) on the accuracy of the height of trees and shrubs using DIM techniques are illustrated in Figure 10.
The height differences indicate that a higher accuracy for determining the height of trees and shrubs is acquired for conifers (Table 5, µ |∆ALS-DIM| = 0.43 m).In the case of deciduous trees and shrubs, the mean value of height differences is 0.73 m; it is strongly differentiated (from 0.37 m to 1.64 m, and µ |∆ALS-DIM| = 0.73 m) and depends on the species and size of the measured object.In general, it was found that tree species, especially those with loose crowns (Betula pendula, Populus tremula), were characterized by higher values of height differences (Figure 10).In the case of Fraxinus sp., Robinia pseudoacacia, and Quercus rubra, larger differences arise from the specifics of objects present in the study area-some of them were characterized by a small crown density (30%-50%) (Robinia pseudoacacia) or smaller crown diameter (Figure 10), moreover Fraxinus sp. and Quercus rubra did not occur in large numbers.Fraxinus sp. and Quercus rubra are rarely found in the analysed area.Fraxinus sp. is usually planted as a roadside tree, Quercus rubra grows spontaneously in abandoned fields and grasslands, and sometimes is planted in forests.Most individuals of both Fraxinus sp. and Quercus rubra growing in abandoned fields and grasslands have well-developed, spherical crowns.But individuals of Fraxinus sp. are usually adult, and individuals of Quercus rubra are mostly young.
Analysis of the similarity degree of heights obtained based on LIDAR data from 2016 and using DIM algorithms (Figure 9c), depending on the size of crowns, showed that trees and shrubs with a crown size above 4 m are characterized by the highest compatibility of both heights (R 2 above 0.9, with the exception of 5 m crown diameter, which is due to the small sample number of such objects).
Analysing the influence of species on the results, for smaller shrubs and brushwood (with a crown diameter less than 2.4 m (Table 7) and height less than 2 m (Table 6)), there was no influence of type, species, or other traits on the value of µ|ALS-DIM|.Both coniferous and deciduous species are at a similar level.For larger trees, µ|ALS-DIM| values are reduced for coniferous species.No influence was observed for the density of crowns of trees and shrubs on the obtained mean differences in height, µ|ALS-DIM (Table 8).Because the analysed area has physiognomically different grades with different heights and crowns, an analysis was carried out with differences in height depending on: heights Table 5. Mean differences between the height of the trees and shrubs obtained using airborne laser scanning (ALS) (2016-reference data) and DIM (2015) taking into account the tree and shrub species.Shrubs are characterized by a much smaller error in determining the height than is the case for woody species (the exception here is Corylus avellana).This applies to both deciduous and coniferous shrubs (Juniperus communis).The average height differences of shrubs obtained based on LIDAR data and using DIM techniques are 0.49 m for deciduous shrubs and 0.34 m for coniferous, respectively.However, considering their lower heights compared to woody species, the relative error of determining the height is quite high, because it reaches 20%-28% of the height of the shrub (Table 6).In the case of trees, this error is generally less than 16%.

Tree and Shrubs
Analysing the effect of object features on the accuracy of the tree and shrub detection, greater variation in the height difference, µ |ALS-DIM| , occurs in deciduous trees and increases with decreasing tree height and reducing crown density (Figure 9a-d).Coniferous trees are generally characterized by a better height mapping than is the case for deciduous trees (Table 5, Figure 9a).It can also be noticed that the larger the tree crown (diameter over 3-4 m), the better the height mapping using DIM techniques (Figure 9c).The graphs in Figure 9 show outliers that represent deciduous trees with a height of approx.5-6 m.These are Betula pendula, Fraxinus sp., and Quercus rubra with tree crowns with 3-4 meter diameters and 50%-60% density (Figure 9c,d).Betula pendula is common in the analysed area.During the secondary succession process, Betula pendula enters abandoned fields and unused grasslands as one of the first tree species.Most individuals have well-developed, slender crowns.Over a dozen-year-olds, medium height individuals of Betula pendula dominate, which is associated with the massive abandonment of cultivation and grazing over the last 20 years.
Fraxinus sp. and Quercus rubra are rarely found in the analysed area.Fraxinus sp. is usually planted as a roadside tree, Quercus rubra grows spontaneously in abandoned fields and grasslands, and sometimes is planted in forests.Most individuals of both Fraxinus sp. and Quercus rubra growing in abandoned fields and grasslands have well-developed, spherical crowns.But individuals of Fraxinus sp. are usually adult, and individuals of Quercus rubra are mostly young.
Analysis of the similarity degree of heights obtained based on LIDAR data from 2016 and using DIM algorithms (Figure 9c), depending on the size of crowns, showed that trees and shrubs with a crown size above 4 m are characterized by the highest compatibility of both heights (R 2 above 0.9, with the exception of 5 m crown diameter, which is due to the small sample number of such objects).
Analysing the influence of species on the results, for smaller shrubs and brushwood (with a crown diameter less than 2.4 m (Table 7) and height less than 2 m (Table 6)), there was no influence of type, species, or other traits on the value of µ |ALS-DIM| .Both coniferous and deciduous species are at a similar level.For larger trees, µ |ALS-DIM| values are reduced for coniferous species.No influence was observed for the density of crowns of trees and shrubs on the obtained mean differences in height, µ |ALS-DIM (Table 8).Because the analysed area has physiognomically different grades with different heights and crowns, an analysis was carried out with differences in height depending on: heights and crown diameters (Table 9), heights and crown density (Table 10), crown diameter and density (Table 11).The results summarized in Table 9 indicate that the mean height differences, µ | ALS-DIM | , for trees and shrubs from a certain height range slightly increase when the crown diameter of a tree or bush decreases.Only in the case of crowns with a larger diameter is there no such a clear tendency, which may result from other features, such as species or shortening of the crown.Figure 10 clearly shows that shrub species (especially Rhamnus cathartica, Juniperus communis, Crataegus spp.) are characterized by large differences in height depending on the diameter and density of the crown.For Crataegus spp., the smallest height deviations were noted for higher crown densities.The largest deviations in height for Juniperus communis were from shrubs with a smaller crown diameter (<1.6 m).The differences observed in height obtained using DIM vs. ALS results mainly from the specificity of these species and their habit.
For smaller shrubs and trees (up to a height of 2 m), there is no significant relationship between crown density and the value of µ |ALS-DIM| (Table 10).Clear trends are visible only for trees and shrubs higher than 2 m.In most cases, for lower crown density (below 60%), differences in height with respect to ALS data are higher, but two exceptions were found here.For example, high values for trees with a height over 10 m and 100% crowns are for individuals with narrow crowns, which affects the results of the DIM algorithm.
A similar situation is seen for the analysis of the relationship between the crown diameter and crown density (Table 11).In the case of small trees and shrubs with crown diameters below 1.6 m, the influence of the crown density on the correctness of vegetation height cannot be noticed.For larger individuals with less than 60% crown density, generally higher differences with the ALS data in the value of µ |ALS-DIM| were found.These results were confirmed with R 2 determination coefficients obtained to determine the degree of similarity of the heights based on the ALS data from 2016 and as a result of the DIM algorithms (Figure 9d).Lower coefficient values were obtained for trees and shrubs with crown densification up to 60%.The Pearson correlation coefficient of trees heights obtained using the DIM and ALS (for all analysed trees and shrubs) is 0.96 (R 2 = 0.92).

Discussion and Conclusions
The experiment presented proved the possibility of using DIM techniques for modelling and detecting individual trees and shrubs, and consequently used this information in assessing the dynamics of the secondary succession process.This process can have different characteristics and dynamics in different areas [51].The course of succession is influenced by, among others, the state of abandoned land (e.g., black fallow, mowed meadow, stubble), soil fertility in nutrients, recently grown plants, the location of the ground (including proximity to the forest), slope, and insolation [52][53][54].This knowledge is important in the context of properly planning active protection of the Natura 2000 habitats.Using dense image matching techniques, it is possible to obtain information on the extent and height of succession trees and shrubs.However, the quality of the results obtained depends on the quality of archival materials and on the parameters for generating the point clouds.When using image matching algorithms during point cloud generation, it is important to select parameters with low smoothing and the densest point cloud generated containing points found for each pixel.There is then a chance to detect individual trees and shrubs.
The most accurate results from DIM data compared to ALS were achieved for data from a modern photogrammetric sensor (digital images from 2015) collected one year before the LIDAR data.The determination coefficient for vegetation height was R 2 = 0.921 for 1000 trees (Pearson coefficient r = 0.96).An equally strong relationship (r = 0.91 for trees younger than 40 years) between the field-based top height and the height calculated from the CHMs (prepared using semi-global matching) was obtained by Stepper et al. [24], Maltamo et al. [55], and Hitara [56].White et al. [57], in their research on Vancouver Island (British Columbia, Canada), also saw a strong relationship between the canopy height obtained from ALS and DIM (r = 0.89).They tested the semi-global matching algorithm for aerial images with GSD = 0.3 m.
In our study, the average height differences between DIM (2015) and ALS (2016) results for coniferous and deciduous trees and shrubs were 0.19 m and 1.45 m, respectively.Worse results for deciduous trees are due to the significant variation in crown density for various species.After eliminating one, Populus balsamifera, with an exceptionally low crown density, the average height difference for the deciduous trees was 0.54 m.Heurich et al. [28] compared the heights measured from LIDAR data and stereo digitalization (DMC camera, GSD = 11 cm) and obtained an average difference for forest stands (subalpine spruce forests with Norway spruce and partly mountain ash; mixed mountain forests with Norway spruce, white fir, European beech and sycamore maple; spruce forests with Norway spruce, mountain ash and birches) equal to 0.99 m (R 2 = 0.94).In the same study, the mean difference of height interpolated from the DSM derived by image correlation vs. stereo digitalization was from 0.51 (R 2 = 0.90, for fagus and acer) to 3.37 m (R 2 = 0.37, for picea, fagus, acer, tilia).The mean deviation between the image-based heights and the ALS-based heights for control polygons obtained by Stepper et al. [31] was 0.17 m.However, in this case, the study area included a broadleaf-dominated (Fagus sylvatica L. 42%, Quercus petraea L. 22%, Pinus sylvestris L. 10%) forest with a great variety of stand development stages.In our research, we analysed the possibility of determining the height of individual (separated) trees and shrubs, which is particularly interesting from the point of view of studying the early stages of the succession process.
Analysing trends in the experiments, better results and relationships between direct measurements from a botanical campaign were for individuals with a large crown and great crown density, depending on the tree and shrub species.Granholm et al. [58] also indicated that the accuracy of tree heights estimated using semiglobal matching (SGM) is low when canopy cover values are low.Regarding LIDAR technology, which also provides a CHM, Leckie et al. [59] achieved better results, similar to terrain measurements, when the tree density is higher.In this research, the average differences in the height of trees and shrubs from a certain height range slightly increase with decreasing crown diameter.A distinct effect of the crown density is visible for trees and shrubs higher than 2 m.There was no clear trend for smaller ones.Higher accuracy height determination (regardless of the date of image acquisition) was obtained for coniferous trees and shrubs (Pinus sylvestris and Juniperus communis), which do not undergo significant changes during the growing season.For deciduous trees, the tree species, especially those with loose crowns (Betula pendula, Populus tremula), were characterized by higher values of height differences.Sexston et al. [60] had the same observations regarding higher accuracy height determination for coniferous trees using remote sensing methods (LIDAR, GeoSAR).
Referring to archival photogrammetric data, the quality of which is related to technology development and time influence, DIM algorithms were also proven to be useful in this application.Based on the extensive analysed material, the effective detection of small trees and shrubs in aerial photographs is possible when individuals taller than 1 m are the subject of analysis.Lower vegetation (small trees and shrubs) solitary located are difficult to find in point clouds, even though they are visible to the photo-interpreter in aerial photographs.DIM also requires certain conditions associated with image parameters.The most important is the scale of aerial photos, which should be greater than 1:13,000 for analogue cameras and have a GSD lower than 25 cm for digital.The best results from a DIM can be found with image acquisition during the period when full foliage is observed (June-September for the tested area).Good radiometric quality of imagery is also required.A small dynamic range, low contrast, or damage to the photos causes errors in the operation of DIM algorithms.

Figure 1 .
Figure 1.Location of the study area (prepared based on UNEP/GRID-Warszawa, © EuroGeographics for the administrative boundaries, orthophotomap prepared by MGGP Aero).

Figure 1 .
Figure 1.Location of the study area (prepared based on UNEP/GRID-Warszawa, © EuroGeographics for the administrative boundaries, orthophotomap prepared by MGGP Aero).

Figure 4 .
Figure 4. Location of measured trees and shrubs as referenced data (on the right, close-up of the study area fragment).

Figure 4 .
Figure 4. Location of measured trees and shrubs as referenced data (on the right, close-up of the study area fragment).

Figure 6 .
Figure 6.The impact of selected cloud parameters on the resulting DSM.Upper left image (a) presents a normalized digital surface model (nDSM) from 2015, the other images show the results of subtracting a DSM from 2015 in approach w1 from: (b) DSM w2, (c) DSM w3, (d) DSM w5.

Figure 6 .
Figure 6.The impact of selected cloud parameters on the resulting DSM.Upper left image (a) presents a normalized digital surface model (nDSM) from 2015, the other images show the results of subtracting a DSM from 2015 in approach w1 from: (b) DSM w2, (c) DSM w3, (d) DSM w5.

Figure 7 . 27 Figure 7 .
Figure 7. Normalized digital surface models (right) and orthophotos (left) for a fragment of the study area-visible errors on DSMs from 1971 and 1982 resulting from the low radiometric quality of aerial photos.

Figure 8 .
Figure 8. Fragment of the orthoimage from 1971 with visible scratches on the film, which influenced the nDSM with outliers.

Figure 8 .
Figure 8. Fragment of the orthoimage from 1971 with visible scratches on the film, which influenced the nDSM with outliers.

Figure 9 .
Figure 9.Comparison of height [m] of trees and shrubs obtained using DIM (2015) and ALS (2016) taking into account: (a) tree and shrub type (deciduous/coniferous), (b) height of trees and shrubs from botanical campaigns, (c) crown diameter, (d) crown density.

Figure 9 .
Figure 9.Comparison of height [m] of trees and shrubs obtained using DIM (2015) and ALS (2016) taking into account: (a) tree and shrub type (deciduous/coniferous), (b) height of trees and shrubs from botanical campaigns, (c) crown diameter, (d) crown density.

Figure 10 .Figure 10 .
Figure 10.The influence of crown diameter and density (by species) on the accuracy of the height of trees and shrubs using DIM techniques.

Author Contributions:
Conceptualization, K.O.-S., K.B., D.M.-H.and D.K.; formal analysis, K.O.-S., .J. and K.B.; investigation, K.O.-S., D.M.-H.and J.W.; methodology, K.O.-S.and K.B.; project administration, K.O.-S.; Supervision, K.O.-S.; visualization, K.O.-S., .J. and K.B.; writing-original draft, K.O.-S., .J. and D.M.-H.; writing-review and editing, K.O.-S., K.B., J.W. and D.K. Funding: The study was co-financed by the Polish National Centre for Research and Development (NCBiR), project No. DZP/BIOSTRATEG-II/390/2015: The innovative approach supporting monitoring of non-forest Natura 2000 habitats, using remote sensing methods (HabitARS).The Consortium Leader is MGGP Aero.The project partners include the University of Lodz, the University of Warsaw, Warsaw University of Life Sciences, the Institute of Technology and Life Sciences, the University of Silesia in Katowice, and Warsaw University of Technology.The leader of WP5 (methodology for the identification of the trees and shrubs succession process in the area of non-forest natural habitats Natura 2000) is the Warsaw University of Technology.

Table 1 .
The characteristics of aerial photos used in the study.

Table 1 .
The characteristics of aerial photos used in the study.

Table 2 .
Settings of variants in the generation of point clouds in Trimble Inpho software (variants w1-w6).

Table 3 .
Parameters for generating digital surface models (DSMs) from point clouds.

Table 4 .
The accuracy assessment of the height obtained based on archival photos using the dense image matching (DIM) technique (versus stereo measurements).Shades of red indicate positive values, shades of blue indicate negative values.

Table 4 .
The accuracy assessment of the height obtained based on archival photos using the dense image matching (DIM) technique (versus stereo measurements).Shades 1 of red indicate positive values, shades of blue indicate negative values.

Table 4 .
The accuracy assessment of the height obtained based on archival photos using the dense image matching (DIM) technique (versus stereo measurements).Shades 1 of red indicate positive values, shades of blue indicate negative values.

Table 4 .
The accuracy assessment of the height obtained based on archival photos using the dense image matching (DIM) technique (versus stereo measurements).Shades 1 of red indicate positive values, shades of blue indicate negative values.

Table 6 .
Mean differences between the height of the trees and shrubs obtained using ALS (2016-reference data) and DIM (2015) considering the height of trees and shrubs.

Table 7 .
Mean differences between the height of the trees and shrubs obtained using ALS (2016-reference data) and DIM (2015) considering the crown diameter.

Table 8 .
Mean differences between the height of the trees and shrubs obtained using ALS (2016-reference data) and DIM (2015) considering the crown density.

Table 9 .
Mean differences between the height of the deciduous trees and shrubs obtained using ALS (2016-reference data) and DIM (2015) considering the tree height and crown diameter.

Table 10 .
Mean differences [m]between the height of the deciduous trees and shrubs obtained using ALS (2016-reference data) and DIM (2015) considering the tree height and crown density.

Table 11 .
Mean differences [m]between the height of the deciduous trees and shrubs obtained using ALS (2016-reference data) and DIM (2015) considering the crown diameter and density.