Countrywide Stereo-Image Matching for Updating Digital Surface Models in the Framework of the Swiss National Forest Inventory

Surface models provide key knowledge of the 3-d structure of forests. Aerial stereo imagery acquired during routine mapping campaigns covering the whole of Switzerland (41,285 km2), offers a potential data source to calculate digital surface models (DSMs). We present an automated workflow to generate a nationwide DSM with a resolution of 1 × 1 m based on photogrammetric image matching. A canopy height model (CHM) is derived in combination with an existing digital terrain model (DTM). ADS40/ADS80 summer images from 2007 to 2012 were used for stereo matching, with ground sample distances (GSD) of 0.25 m in lowlands and 0.5 m in high mountain areas. Two different image matching strategies for DSM calculation were applied: one optimized for single features such as trees and for abrupt changes in elevation such as steep rocks, and another optimized for homogeneous areas such as meadows or glaciers. The country was divided into 165,500 blocks, which were matched independently using an automated workflow. The completeness of successfully matched points was high, 97.9%. To test the accuracy of the derived DSM, two reference data sets were used: (1) topographic survey points (n = 198) and (2) stereo measurements (n = 195,784) within the framework of the Swiss National Forest Inventory (NFI), in order to distinguish various land cover types. An overall median accuracy of 0.04 m with a normalized median absolute deviation (NMAD) of 0.32 m was found using the topographic survey points. The agreement between the stereo measurements and the values of the DSM revealed acceptable NMAD values between 1.76 and 3.94 m for forested areas. A good correlation (Pearson’s r = 0.83) was found between terrestrially OPEN ACCESS Remote Sens. 2015, 7 4344 measured tree height (n = 3109) and the height derived from the CHM. Optimized image matching strategies, an automatic workflow and acceptable computation time mean that the presented approach is suitable for operational usage at the nationwide extent. The CHM will be used to reduce estimation errors of different forest characteristics in the Swiss NFI and has high potential for change detection assessments, since an aerial stereo imagery update is available every six years.


Introduction
Knowledge of tree resources and the structure of forests is key for managing various forest functions and services.Information on the third dimension (z component) of 3-d data has become more and more important over the last decade.Forest characteristics are traditionally estimated in forest inventories based on sample plots.The combination of design-based sample inventories and remote sensing data, e.g., satellite data or airborne laser scanning data, reduces estimation errors and makes field-based inventories more efficient [1].Stereo optical imagery or laser scanning data have tremendous potential for obtaining area-wide three-dimensional data over several thousands of hectares.Their usage in the analysis of forested areas over an entire country opens up new possibilities for biomass estimations [2,3], forest structure analysis [4,5], biodiversity assessment [6,7] and the optimization of forest inventories at a national level [8].
Three-dimensional data can either be acquired with an active laser or radar system, or passively with a camera taking optical images [9].In forested areas, airborne laser scanning (ALS) is often the method of choice, because ALS can penetrate the forest canopy and is, therefore, the only means of obtaining sub-canopy forest structure and terrain information [10].ALS data has proven to be the most accurate data for the generation of highly detailed digital elevation models of bare earth, referred to as digital terrain models (DTMs), within forests [11].Several countries now have or are in the process of acquiring national ALS data to provide this crucial terrain information countrywide (for example Sweden, Netherlands, Denmark, Austria, Switzerland, Slovakia, Czech Republic and Finland).For the generation of digital surface models (DSMs) different input data can be used.For instance, surfaces may be modeled using the first-return of the ALS point cloud or using point clouds coming from stereo aerial or satellite images.
Once a detailed DTM is available, canopy height models (CHMs) can be computed by subtracting bare earth heights from the DSM heights.If the terrain does not change significantly, which is the case in most countries, the same ALS DTM can be used for the generation of CHMs at several time steps.Digital airborne imagery over large areas is typically less expensive to acquire than ALS data [10].The acquisition of aerial stereo imagery has a long tradition in many countries and such images are used to update topographic maps and for orthophoto production.There are a variety of studies based on digital aerial images for canopy height generation which have been mainly carried out in Canada [12,13], Germany [14][15][16], Sweden [17] and Finland [18][19][20][21].
Since the stereo-images are acquired during regular national campaigns, there is a constant cycle of image acquisition allowing for continued update of CHMs.This updating is very valuable for studies of environmental processes where change detection is of high importance [22].Since the workflow for the photogrammetric stereo image matching can be highly automated, a constant updating of the CHMs is also feasible.In the case of Switzerland, the national aerial imagery campaign is carried out by the Federal Office of Topography (swisstopo) using Leica ADS 40 and ADS80 sensors.In a cycle of maximum six years, the whole country, with an area of 41,285 km 2 , is updated with aerial image data during the leaf-on season.
A variety of different algorithms and software packages for high density image matching are available today.The availability of high quality digital image data, efficient image matching algorithms and the ever increasing capabilities of CPUs and GPUs have led to major improvements in the generation of dense point clouds from stereo images.We decided to use the SocetSet "Next-Generation Automatic Terrain Extraction" (NGATE) by BAE systems for this study, because we have a long time experience with this software and a large amount of licenses at our institute.NGATE includes area-and feature-based image correlation at different minification levels, while the feature-based method matches, for example, edge segments and contours, the area-based approach attempts to correlate the grey levels of image patches based on the assumption that they represent some similarity [23].The two image-matching approaches of NGATE complement each other.Although the matching performance of NGATE would also improve from higher image overlap, it achieves a good compromise between matching quality and calculation time.
To date, most of the studies using image-based point clouds for forest applications have been carried out in small study areas with relatively flat terrain and even-aged forests [10].In contrast, we present a nationwide assessment over a variety of forests for the first time.Although Switzerland, with an area of 41,285 km 2 , is a relatively small country, it is dominated by a small-scale mixture of different forest types covering one third of its territory, and is consequently an excellent testing ground for our approach.As vegetation height is the most important proxy for a wide range of forestry variables, we focus in our analysis on the accuracy and agreement of modeled tree and canopy heights.The main objectives of the present study are: (1) to demonstrate the highly automated workflow to calculate a nationwide surface model out of regular national aerial images; (2) to evaluate the accuracy and completeness of this digital surface model, which is a prerequisite for its subsequent usage; and (3) to present the potential of this data for National Forest Inventories.

Study Area
Switzerland is a mountainous country of 41,285 km 2 area with a diverse topography covering an altitudinal range between 200 and 4500 m a.s.l.[24].The highest elevations can be found in the Alps in the south and southeastern part of Switzerland.Running from west to east, on the northern side of the Alps, is the Swiss Plateau, where most of the settlement and urban areas are located.The Jura Mountains, which are generally lower than the Alps, are located on the northwestern border of Switzerland and are more rural in character.
Switzerland is characterized by a mixture of different land uses, which are described in the federal land use statistics [25].Four broad designations are distinguished: settlement and urban areas (8% of the surface area), agricultural areas (36%), wooded areas (forest and woods, 31%) and unproductive areas (lakes and rivers, unproductive vegetation, rocks and screes, glaciers and perpetual snow, 25%) [25].Forested areas consist of 43% pure coniferous, 18% mixed coniferous, 12% mixed deciduous, 24% pure deciduous and 3% other forests according to the third NFI, covering the period of 2004-2006 [26].

Image Data
To obtain nationwide coverage, stereo aerial images from regular flight campaigns of the Federal Office of Topography (swisstopo) were used.Swisstopo has two flight programs, one for updating of the topographic landscape model (TLM) [27], usually acquiring images in the leaf-off season in order to also map infrastructure in forests, and one for orthophoto production [28], usually in the leaf-on season (May-September).The strategy is to cover 1/6 of Switzerland for each program annually.This implies that for a complete leaf-on image coverage of the country a maximum six year period has to be considered.Consequently, for the purpose of the present study, image data from 2007 to 2012 were used (Figure 1).Over this time frame, three different sensors were in operation.2007 was the last year swisstopo used the ADS40-SH40 sensor.This sensor collects, at nadir (0°), three spectral bands (blue, green and red), and forward-looking red, green and near infrared (14°, 16° and 18°).Two Pan bands are collected by looking −14° backwards and 28° forwards.The ADS40-SH52 sensor has been used since 2008, with an additional ADS80-SH82 sensor since 2009.These sensors collect four spectral bands (blue, green, red and near infrared) at nadir (0°) and four spectral bands (blue, green, red and near infrared) by looking backward (−16°).Technical details for the different sensor configurations can be found in Appendix A, Table A1.
Swisstopo makes use of two different ground sample distances (GSDs) (Figure 1).In the Jura Mountains, the Swiss Plateau, the pre-alps and along mountain valleys, images have a GSD of ~25 cm.In the Alps, the GSD is ~50 cm.The photogrammetric pre-processing (level 1 generation of raw image data, aerial triangulation (AT) and calculation of RGB and CIR band combinations) is done by swisstopo with the Leica GPro 4.2 software package and image data is provided including camera calibration and orientation data.The orientation accuracy (x, y, z) is reported to be ±1 pixel GSD.

Generation of the Digital Surface Model (DSM)
The entire workflow for the generation of surface models is highly automated (Figure 2).To distribute calculations over a range of different processing machines, the country was divided into 165,500 blocks of 0.5 × 0.5 km.To avoid edge effects, a buffer of 100 m was added to all sides of the blocks, resulting in an area of 0.7 × 0.7 km for the actual image matching.The blocks were intersected with the stereo footprints of the ADS image stripes, and the distances from the center of each block to the flight lines of the ADS image stripes were calculated.Since the stripes overlap by 50%, in the best case, two stripes for one block are available: one stripe with a shorter distance from the center of the block to the nadir line of the stripe (hereafter 1st most nadir stripe) and a possible second stripe, with a larger distance to the nadir line of the stripe (hereafter as 2nd most nadir stripe).The 1st most nadir stereo image stripes had the highest priority for the matching process for each block since obliqueness effects due to the central perspective of the images are reduced.The 2nd most nadir stereo image stripes were used in cases where image matching on the 1st most nadir stereo image stripes was not successful.The list of all image blocks was managed in a central database to maintain an overview of the ongoing image matching processes.When matching started for a particular block, it stopped being available for other processing machines.After the matching procedure was completed, the block was labeled as finished.Sixteen processes (using 16 virtual Windows7 64bit PCs with 2 cores, 2GB RAM and 2.67 GHz on a Citrix Xen environment) were running parallel on two HP Pro Workstations.The image data was available through a 1 Gb intranet from a central image server.
Image matching was performed using SocetSet 5.6 digital mapping software.The combination of the SocetSet "Next-Generation Automatic Terrain Extraction" (NGATE) by BAE Systems and the stereo line scanning data from the ADS sensors gave the best performance with respect to quality and completeness of image matching and correlation speed.As part of the NGATE, matching is performed for every pixel using area-matching and edge-matching algorithms [29,30].Although NGATE computes an elevation for every pixel, a fixed spacing must be defined to resample the internal irregular NGATE point cloud to the desired density.No information is available on the interpolation algorithm.Spacing was set to 1 m for the images with 0.25 m and 0.5 m GSD.Only matched raster points were exported from SocetSet as ASCII and then converted to LAS format.the stripe with the smallest distance from nadir to the centre of the block; and (C) selecting the stripe with the second smallest distance from nadir to the center of the block (1) The result after the matching process for the block using the 1st most nadir stereo stripe and the "ADS steep" image matching strategy.Red areas in the image could not be matched.The completeness for this block is 84%.(2) The result after applying the "ADS flat" matching strategy.At 91%, the completeness threshold is still not fulfilled.(3) The result after matching using the 2nd most nadir stripe and the "ADS steep" strategy (completeness 96%) and (4) the final DSM after "ADS flat" strategy with a matching completeness of 96%.
Within SocetSet several correlation strategies for different image contents (e.g., urban, flat homogeneous, hilly and steep areas) are provided.To reach a high level of completeness of the matched points, some of these matching strategies were adapted.We used two different matching strategies (Table 1): "ADS steep", which is optimized for steep areas and very heterogeneous surfaces with a broad range of different slopes (Figure 3a,b), and "ADS flat", which is optimized for homogenous areas with little contrast (Figure 3c,d).When matching with the 1st most nadir stereo image stripes using the "ADS steep" strategy was finished, the proportion of successfully matched points was calculated.In one block of 500 × 500 m the maximum of matched points at a grid spacing of 1 m is 250,000.We defined the level of completeness to be reached by 99.5%.If this completeness value was not reached, the "ADS flat" strategy was applied.All additionally matched points with this strategy were then added to the point cloud resulting from the initial "ADS steep" strategy.If the proportion of successfully matched points was still below the threshold of 99.5%, the same procedure was followed using the 2nd most nadir stereo image stripe.Information regarding the stereo image stripe and the strategy used for its calculation was then assigned to each point of the point cloud.
Table 1.Parameter sets of the two image matching strategies used in the image correlation workflow.The algorithm iterates seven times using the pyramid levels of the images.The seven values listed in the table correspond to the number of iterations.105 (1) The actual window correlation size corresponds to the difference of parameter (3) and parameter (10).
A smaller window is used in areas with many elevation discontinuities.
Image-matching was done using only one image pair at a given time, which allowed for better control over the amount of matched pixels.Since the green-red-near infrared (CIR) band combination produced the least number of errors, especially in forested regions and in high alpine glacier regions, the two CIR stereo images (backward 16° and nadir 0°) were used for image matching for all of Switzerland.With the exception of the case of the ADS40-SH40 image data from 2007 where the combination of PAN (backward 14°) and RGB (nadir 0°) was used for the matching, because CIR band combination was not possible with this focal plate configuration.

Topographic Survey Points
A total of 198 topographic survey points, planimetic fixpoints of the new Swiss national topographic survey LV95 [31], were used to assess the absolute accuracy of the digital surface model (Figure 4).They cover the whole country on a more or less regular grid, and are reported to have an absolute accuracy of 3 to 5 cm.

Stereo Measurements of the National Forest Inventory
In the framework of the Swiss NFI, continuous landscape variables are interpreted on the same ADS40/ADS80 stereo imagery, which was used for the digital surface modeling.These stereo measurements are carried out on a rectangular sampling grid of 1.4 km × 1.4 km, using interpretation clusters of 50 m × 50 m, for details see [32,33].Each interpretation cluster consists of 25 equally spaced lattice points (see Appendix F).A photo interpreter measures the elevation of each lattice point hitting the surface using a 3D softcopy station and assigns a land cover class to each point [34].For this study, 195,784 measured lattice points, covering 11 land cover types, were used to assess the agreement of the digital surface model.Image matching of water surfaces was not possible since spectral signatures of neighboring pixels were very similar or waves on the water surface led to different reflectance within an image pair; as such this land cover type was excluded from further analyses (Table 3).The values of the four pixels surrounding the lattice points of the stereo measurements were extracted from the digital surface model for comparison with the stereo measurements.For objects with a low object height (herb and grass, rock, sand, stones and visible earth, sealed surfaces, glaciers, firn and snow), the values of the stereo measurements were compared with the minimum value of these four pixels, and for objects with a high object height (trees, shrubs and buildings), the maximum value was used.Only values where the same ADS40/ADS80 image stripe was used for the stereo measurements and the generation of the digital surface model were used in the agreement analysis.

Accuracy Assessment of the Digital Surface Model (DSM)
Knowledge of the accuracy and completeness of the nationwide DSM over different land cover types is crucial for its subsequent usage.Therefore, vertical accuracy assessments were carried out using two reference data sets (topographic survey points and stereo measurements).First, the difference in elevation between the digital surface model and the reference data set was calculated.Because the DSM deviations were not normally distributed, e.g., [35,36] robust accuracy measures described by Höhle and Höhle [37] were used.These included the median, the normalized median absolute deviation (NMAD) and the 68.3% and 95% sample quantiles.For better comparison with other data sets two root mean square errors (RMSE) were calculated, one for the whole sample size and one with outlier correction, where a threshold value of ±3 × RMSE was applied for outlier exclusion.All calculations were done with the open source statistical software R [38].Analyses were conducted only for the sample of successfully matched image points.

Generation of the Canopy Height Model (CHM)
By subtracting the digital terrain model swissALTI3D provided by swisstopo from the calculated digital surface model, a normalized digital surface model (nDSM) was generated.The swissALTI3D of swisstopo is based on ALS data with an overall point density of 0.5 points per m 2 [39], has a resolution of 2 m which was resampled to 1 m.To obtain the canopy height model, the topographic landscape model TLM provided by swisstopo was used to mask out houses, settlement areas and water bodies.

Tree Height of Single Trees
Within the ongoing 4th inventory of the Swiss NFI, the heights of a total of 3109 trees belonging to the upper canopy layer were measured using a Vertex ultrasonic hypsometer.Field measurements were conducted between 2009 and 2013, and are spread over the whole country.The coordinates of the tree stem location at breast height were acquired with a Trimble GeoXH (with an average accuracy after post-processing of ±1 m), enabling a comparison with the values of the CHM.Within a 2.5 m radius buffer around each stem, the maximum height values of the CHM were extracted and compared to the tree height of the corresponding measured trees.Only those trees with at least 15 pixels that could be matched within the 2.5 m buffer (93% of all measured trees) were used for the analysis.In addition, only measurements where the tree height was measured later than the date of the image acquisition were included in the analysis.To know how precisely the trees were measured in the field, double measurements within the NFI of 117 deciduous and 324 coniferous trees were compared.

Canopy Height at Plot Level
Individual tree heights measured in the 4th NFI were used to estimate the maximum canopy height at the circular NFI plot level of an area of 500 m 2 .The maximum values of these terrestrial measurements (only a sub-sample of the trees on the plot is measured-one to three upper layer trees per plot) were then compared to the maximum values of the CHM within the sample plot area.We analyzed the data for different development stages, species mixtures and stand structures.

Generation of Digital Surface Model (DSM)
The nationwide DSM with a resolution of 1 × 1 m was successfully generated.To illustrate this, we provide three examples below: (a) an urban area within the city of Zurich, (b) a typical forest of the Swiss lowlands and (c) a mountain area with glaciers (Figure 5).The examples show how successful the matching in these different land cover areas was.In the example of the city of Zurich, we see that the buildings could only be roughly mapped and their shape is not well represented.Looking at the forested area, we can visually detect the different developmental stages of the forest and some gaps within the forest canopy.Even in the high mountain areas where snow and ice dominate, small structures are visible in the DSM.To cover the entire country, a total of 166,725 image blocks were matched: 54.5% with a GSD of 0.25 m and 45.5% with a GSD of 0.5 m.The different strategies and the combination of the 1st and 2nd most nadir image stripes allowed for an overall completeness of the image matched points of 97.9% (Table 2).Already using only one strategy with one image stripe enabled an overall mean completeness of all image blocks of 94.9%.The greatest improvement was achieved while additionally processing the "ADS flat" strategy, which led to a mean gain in completeness of 1.7%.To achieve this rather minor improvement in completeness, 133,897 blocks were matched with the "ADS flat" strategy, requiring a calculation time of 98 days.For the year 2007, the gain in completeness by adding this second strategy was a very large improvement of 6.5% (see Appendix B, Tables A2 and A3).Manual corrections of the data were done for only 0.14% of the matched blocks, most of which were done because of cloud cover.Manual correction involved visual inspection and replacement of those results showing a high proportion of errors with data from earlier years.
The mean calculation time per image block was 16 min for one strategy.In our study, the matching process, distributed to 16 virtual PCs, was completed within approximately 320 days (including ~10% for database updates and network maintenance).Thus, within two months, the annually acquired image data for one sixth of the country (~7000 km 2 ) can be processed and an updated CHM can be provided.Abandoning the second strategy would reduce the calculation time from 320 to 218 days, corresponding to a reduction of 32%.

Topographic Survey Points
Within the national triangulation network, all 198 topographic survey points could be successfully used for accuracy assessment (Table 3).The median deviations, with values of 0.07 and −0.11 m, are well below one pixel and result in a normalized median absolute deviation (NMAD) of little more than one pixel.The positive median values in the case of 0.25 m GSD indicate lower elevation values of the topographic survey points in comparison to those of the digital surface model.For images with 0.5 m GSD, the opposite was the case.The combination of the DSM and the stereo measurement data enabled an agreement assessment over various land cover types (Table 4).The lowest median and NMAD values were found for sealed surfaces.NMAD values below 1.5 m were found for all land cover classes except forested areas.Forest trees were the most difficult ones to model, as indicated by high NMAD values varying between 1.76 and 3.94 m.Additionally, the influence of the nadir distance on accuracy measurements was evaluated and median errors were found to rise slightly with the increase in distance to the nadir line of the image stripe.This result was independent of the land cover type.Glaciers, firns and snow in the high mountain areas were only covered with images of 0.5 m GSD.We found increasing NMAD values with increasing slope inclination for both categories 0.25 and 0.5 m GSD (Table 5).The same trend was found for the RMSE without outliers and the quantile values.Contrary to the effects of slope inclination, different aspects were found to have no influence on the accuracy of the digital surface model (Appendix C, Table A4).

Tree Height of Single Trees
The comparison of the terrestrially measured tree heights with the tree height calculated from the CHM revealed acceptable results.For areas in the lowlands with 0.25 m GSD, deciduous trees showed a lower median but a higher NMAD than coniferous trees (Table 6).In the mountains where the GSD was 0.5 m, deciduous trees still showed a lower median, but the NMAD was higher for coniferous trees.The higher quantile values for deciduous trees indicate that these trees tend to be more difficult to measure terrestrially.This result is consistent with the comparison of the 441 double tree height measurements (Appendix D, Table A5).There was a high correlation (Pearson's r = 0.83) between the terrestrial measurements and the canopy height model values (Figure 6).There was no difference in correlation between coniferous and deciduous trees (Appendix D, Figure A1).A slight difference in correlation was found between two different elevation classes (Pearson's r = 0.72 for low elevations and Pearson's r = 0.6 for high elevations, Appendix D, Figure A2).Outliers occurred in low numbers.They indicate cases where the forest changed between the date of the image acquisition and the terrestrial height measurement or where the digital surface model was not successful in modeling individual tree crowns.

Canopy Height at Plot Level
Good agreement between the canopy height measured terrestrially and the modeled canopy height at the plot level was found over a wide variety of forest stands (Figure 7, Appendix E, Table A8).If we focus on the different development stages of the forest stand, pole wood was found to have the highest correlation (Pearson's r = 0.91) and mixed stands the lowest (Pearson's r = 0.64).Regarding species mixture, the correlation was dependent on the percentage of coniferous trees within the stand, showing highest values in stands with over 90% coniferous trees (Appendix E, Figure A3, Table A6).Stand structure influenced correlation as well, showing highest values in one-and multi-layered stands and lower values in highly structured stands (Appendix E, Figure A4, Table A7).

Workflow for the Generation of the Models
The highlight of the nationwide digital surface generation approach presented here is the consistency of the input data and matching approaches over an entire country.The workflow used is highly automated and allows for an update of the DSM covering Switzerland every three to six years.The input image data used in this study is routinely acquired for different purposes of countrywide orthophoto production and topographic mapping.Therefore, there are limitations to the spatial and temporal resolution of the data, in comparison to projects with flight campaigns planned solely for a specific research objective.However, the great advantage of the approach lies in the opportunity for model updating that is both effective and cost-efficient.Once a digital terrain model based on ALS data is available, the photogrammetric updating of digital surface data and canopy height models using routinely acquired aerial stereo image data provides a valuable alternative to the acquisition of repeated ALS data.
In the present study, we used the SocetSet "Next-Generation Automatic Terrain Extraction" (NGATE) by BAE Systems for image matching.It enabled consistency in calculations over the different time periods and achieved a good compromise between image-matching quality and calculation time.Other software packages, such as Trimble Match-T and Leica XPro using Semi-global matching were tested in case studies.While good overall image-matching results were achieved, calculation times were at least four times as high as with NGATE implemented in SocetSet.Inaccuracies in the surface models were found for situations such as those listed by Baltsavias et al. [40], for example water surfaces, where matching was not successful due to the movement of the water surface between the two time steps of image acquisition, as well as on steep slopes or in the high mountains where whiteout was an issue.
A high level of completeness was achieved using two different matching strategies ("ADS steep" and "ADS flat") on two pairs of imagery (1st most nadir and 2nd most nadir stripe).The combination of these strategies and input data is important to cope with the heterogeneous mosaic of different land cover types within the country.Although using one strategy applied only to one image pair would reduce the calculation time substantially (by 50%), we maintained the consistent application of the different strategies for the matching of the 1st and 2nd most nadir image stripes for each block for the whole country.The main reason was that we were striving for a very high level of completeness to allow for a broad spectrum of analysis based on this data set.The use of a second image pair was important for areas where clouds or haze led to occlusion of some image parts in the first image pair, which could be matched using the second image pair.In addition, in the case of forested areas, the use of a second image pair allowed for another view angle, for example, for oblique trees at forest edges or within canopy gaps, resulting in a substantial improvement in the image matching for those areas, even if the overall gain in completeness was not high.The combination of the two strategies was important for a successful matching of areas with low texture, such as expanse meadows, grasslands or snow areas.In these areas, the use of the second strategy enabled the matching of areas that could not be matched with the first strategy, thus leading to higher overall matching success.Especially in the case of stereo image pairs from the year 2007, where a combination of RGB and Pan images was used, this second strategy led to a substantial improvement in overall completeness.

Accuracy of the Digital Surface Model (DSM)
The accuracy of the calculated DSM is very good in the case of sealed surfaces (e.g., median difference of −0.04 m and NMAD 0.52 m for 0.25 GSD) and acceptable for most applications in forested areas (e.g., median difference of 0.16 m and NMAD of 2.95 m for deciduous trees, 0.25 GSD).A similar result for a small study area in northeast Switzerland was achieved in previous work using GPS points as reference with a median difference of 0.08 m and a NMAD of 0.21 m in flat terrain and a median difference of −1.12 m with a NMAD of 1.32 m in forested areas in comparison with stereo-measurements [41].The deciduous tree species in this study had a greater median difference and a higher deviation than the coniferous trees.In the present study, the deviation of the values in forested areas was larger than in the above mentioned study, since calculations were made over the entire country spreading over an area of a 41,285 km 2 at an altitudinal range of 200 to 4500 m a.s.l., and including a variety of forest types and tree species.There are no studies outside Switzerland based on ADS stereo aerial images to which the results of the present study could be compared.In a study of a broadleaf-dominated forest in central Europe, average differences of 0.17 m with a standard deviation of 0.29 m between the UltraCamXp image-based heights and the ALS-based heights were achieved for flat areas [15].These values are comparable to the accuracy values of the topographic survey points in Switzerland.
To take into account variation in land cover in Switzerland, a separate analysis was performed distinguishing different land cover types based on stereo measurements.This study revealed median differences between 0.04 and 0.80 m with NMAD values ranging from 0.49 to 3.94 m.Herb and grass, and sealed areas had the highest agreement values with NMAD values of ±2 pixels for both GSD classes.Shrubs, buildings, glaciers, firn and snow, with NMAD values of ±3 pixels, showed good agreement values as well.NMAD values of up to ±5 pixels were found for areas covered by rock, stands, stones and visible earth.As expected and mentioned above, forested areas were the most difficult areas to model, achieving NMAD values of up to ±16 pixels.Looking at the influence of topography, a decrease in accuracy with increasing slope inclination was found.This is in line with a study by Bühler et al. [42] in high alpine terrain comparing ADS80 stereo image elevation values and ALS data and with the findings of Müller et al. [43] in a high mountain environment.Contrary to our expectations, no difference in accuracy was found for different aspect classes.Due to the high radiometric resolution of 12 bit, even in shadowed parts in north-facing mountainous regions, the matching was successful.
The two reference data sets used complemented each other and allowed for comprehensive assessment of the accuracy of the DSM.The topographic survey points are a completely independent data set enabling an estimation of the absolute accuracy.However, this analysis is largely dependent on the orientation of the input imagery.Small shifts in geolocation could lead to higher vertical errors.For stereo measurements, a shift in geolocation can be excluded, since these measurements were carried out on the oriented image data itself.This data set is highly suitable for the assessment of the agreement with the digital surface model over different land cover types, but cannot be used for an absolute accuracy analysis in a strict sense.

Applications within the National Forest Inventory
An important characteristic of this study is the variety of forests analyzed with regards to development stages, species mixture and stand structure.Given the very heterogeneous data set, the comparison of single tree heights and canopy height at the plot level revealed plausible results and demonstrates the potential of the canopy height model.This variation in analyzed forest types is also the reason why the RMSE values are substantially higher than some other relevant studies, for example, a Swedish study that used DMC aerial images in a coniferous forest for DSM generation, where RMSE values between 1.4 and 1.6 m were found [17], or a Finnish study based on UltraCamD images, where RMSE values ranged between 6.77% and 7.55% of tree height [19].
The accuracy of a CHM is largely dependent on the DTM used for its generation.In our case, the ALS DTM is from the years 2001-2003 and the ADS40/ADS80 images from the period 2007-2012.The time differences can be ignored, because in most areas the terrain has not changed significantly over the short period of several years.Major errors, however, can occur if there are differences in the geolocation of the DSM and the DTM, or the two models are not accurately co-registered [18].To avoid the influence of such minor shifts in geolocation, we always used some kind of buffer while comparing different data sets for evaluation purposes.
These results show that the canopy height model is most usable for stand-wise assessment, where the focus is on estimates over the whole canopy, and the influence of individual trees is only minor.The data is less suitable for estimations of single tree parameters and estimations based on individual trees, since the spatial resolution of the input image data is too low for such specific purposes.In the framework of the Swiss National Forest Inventory, the data is used as auxiliary information to post-stratify forest areas according to canopy height for the estimation of different forest characteristics, see also [8,44].This allows for a further reduction of estimation errors, especially when it comes to small area estimations [1,45].

Conclusions
For the first time, a digital surface model with a resolution of 1 × 1 m and a canopy height model covering an entire country have been calculated using an image matched point cloud with a very high completeness of matched points of 97.9%.Switzerland, at just over 41,000 km 2 , is a rather small country, although its mixture of different land cover types, an altitudinal gradient of 200 to 4500 m a.s.l., and the complex topography makes the country particularly interesting for this kind of research.Improvements in matching algorithms, combined with the increasing performance in terms of hardware capacity make nationwide processing of aerial stereo imagery possible, and with the highly automated workflow presented, the manual interaction during the process is reduced to a minimum.This allows for an operational generation of a 3-d data set of homogeneous data quality over a heterogeneous country, even if calculation time is still high.The agreement between the manual stereo measurements and the elevation of the calculated digital surface model at trees (NMAD below +/− 16 pixel) as well as the good correlation (Pearson's r = 0.83, n = 3109) of terrestrially measured trees and the tree height derived from the canopy height model show the high potential of the model for country-wide forest applications.
Using routinely acquired stereo imagery, an update of models can be provided at regular intervals.For 81% of the country at least two time steps are already calculated.Initial studies for change detection are underway.This data set provides high potential to estimate changes in timber volume, gaps within the canopy or development stages of forest stands.A primary future application lies in the estimation of the extent of damage after a large windthrow, fire or insect disturbance within the forest, based on stereo imagery which is acquired immediately after the natural disturbance event.

B. Completeness
For all image data sets of the different acquisition years, completeness is over 90% after the first image strategy used for each block.In the case of the year 2007, values are slightly lower in comparison to other years, which can be related to the use of a different combination of spectral bands (Table A1).The greatest gain in completeness was achieved by using the second strategy on the same image stripe (Tables A2 and A3).A similar improvement was achieved when using a second image stripe and only a minimal gain in completeness was found by using the second strategy on the second image stripe.The overall gain in completeness using the different strategies and two image stripes was 3% in the years from 2008 to 2012 and 8.1% for the year 2007.

C. Accuracy Assessment of DSM
Other than the slope, aspect was found to have only a minor influence on the accuracy of the DSM (Table A4).One could expect that northern exposed slopes affected by shadow effects would show larger errors.This, however, did not occur, indicating that the 12 bit ADS40/ADS80 data allow for good image matching even in shadowy areas.

D. Comparison of Individual Tree Height Measurements
We found no difference in the correlation between tree height extracted from the CHM and terrestrial height measurements when distinguishing between coniferous and deciduous trees (Figure A1).However, as expected, the correlation was higher for trees in the lowlands than for trees at higher elevations (Figure A2).Double measurements of trees within the Swiss NFI were used to estimate how these differences in correlations are influenced by the accuracy of tree height measurement in the field (Table A5).As expected, coniferous trees could be measured more accurately than deciduous trees.

E. Canopy Height at the Plot Level
At the plot level, species mixture had some influence on the agreement of tree height in the CHM and the terrestrially measured tree heights (Figure A3).The best correlation was found for stands dominated by over 90% coniferous trees.Stand structure influenced the correlation as well (Figure A4).Stands with clear layering could be better mapped than stands with a high variation of canopy height at the plot level.

Figure 1 .
Figure 1.Image data used: (a) Coverage of ADS40/ADS80 images acquired in 2007-2012, (b) flight lines for the image acquisition for the entire of Switzerland.The imagery was acquired with a GSD of 0.25 m at lower altitudes (denser flight lines) and a GSD of 0.50 m at higher altitudes in the mountainous regions (less dense flight lines).

Figure 2 .
Figure 2. Image matching workflow.(A) Rasterizing into 0.5 × 0.5 km blocks; (B) selecting the stripe with the smallest distance from nadir to the centre of the block; and (C) selecting the stripe with the second smallest distance from nadir to the center of the block (1)The result after the matching process for the block using the 1st most nadir stereo stripe and the "ADS steep" image matching strategy.Red areas in the image could not be matched.The completeness for this block is 84%.(2) The result after applying the "ADS flat" matching strategy.At 91%, the completeness threshold is still not fulfilled.(3) The result after matching using the 2nd most nadir stripe and the "ADS steep" strategy (completeness 96%) and (4) the final DSM after "ADS flat" strategy with a matching completeness of 96%.

Figure 3 .
Figure 3. Examples of the performance of the two different correlation strategies: "ADS steep", optimized for steep areas and very heterogeneous surfaces such as forest canopies (a) input CIR image and (b) shaded relief of the resulting DSM and "ADS flat", optimized for homogeneous areas with little contrast such as glaciers (c) input CIR image and (d) shaded relief of the resulting DSM.

Figure 4 .
Figure 4. Reference data for accuracy assessment.The 198 topographic survey points of the national triangulation network were used for an independent error analysis.

Figure 5 .
Figure 5. Shaded relief of the generated DSM with a colour gradient for canopy height (0-50 m) on the left and the used ADS80 input image data on the right: (a) buildings in Zurich, (b) forests in the Swiss lowlands and (c) glaciers in the Valais.

Figure 6 .
Figure 6.Correlation between tree heights of the CHM and terrestrial tree height measurements.Total number of measurements is 3109.Tree height in the canopy height model refers to the maximum value within a buffer of 2.5 m radius around each individual tree.Tree height in the field was measured using a Vertex ultrasonic hypsometer.

Figure 7 .
Figure 7.Comparison of canopy height from the CHM with terrestrially measured canopy height at the NFI sample plot level.One of the four stages of stand development was assigned to each sample plot.(a) Young growth refers to regeneration areas with a medium tree height of 8 m; (b) pole wood refers to stands with a high number of stems and a mean tree height of 8-20 m; (c) the timber stage is characterized by closed stands with a mean height > 20 m; and (d) mixed refers to stands where there are different stages of stand development.

Figure A1 .
Figure A1.Correlation of tree height extracted from the CHM with terrestrial tree height measurements distinguishing between (a) deciduous and (b) coniferous trees.The total number of measurements was 3109.Tree height in the canopy height model refers to the maximum value within a buffer of 5 m diameter around each individual tree.Tree height in the field was measured using a Vertex ultrasonic hypsometer.

Figure A2 .
Figure A2.Correlation between tree height extracted from the CHM and terrestrial tree height measurements distinguishing between (a) lower and (b) higher elevations.Total number of measurements was 3109.Tree height in the canopy height model refers to the maximum value within a buffer of 5 m diameter around each individual tree.Tree height in the field was measured using a Vertex ultrasonic hypsometer.

Figure A3 .
Figure A3.Comparison of canopy height of the CHM with terrestrially measured canopy height at the sample plot level of the NFI.One of the four species mixture classes (a-d) depending on the percentage of intermixed coniferous trees was assigned to each sample plot.

Figure A4 .
Figure A4.Comparison of canopy height of the CHM with terrestrially measured canopy height at the sample plot level of the NFI.One of the four stand structure classes (a-d), defined by the layering and structure of the canopy was assigned to each sample plot.

F
Figure A5.Sample plots of the Swiss National Forest Inventory spaced in a regular 1.4 km grid, consisting of 25 lattice points each where height and surface cover information was measured.

Table 2 .
Gain in completeness of matched points with the different strategies and image stripes used.

Table 3 .
Vertical accuracy measurements of the digital surface model based on 198 swisstopo topographic survey points.Differences were calculated by subtracting the elevation values of the topographic survey points from those from the digital surface model.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table 4 .
Vertical agreement measurements of the digital surface model per land cover type based on stereo measurements.Differences were calculated by subtracting the elevation values of the stereo measurements from those of the digital surface model.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table 5 .
Vertical agreement measurements of the digital surface model based on stereo measurements covering five slope classes.Differences were calculated by subtracting the elevation values of the stereo measurements from the values of the digital surface model.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table 6 .
Agreement of tree height calculated from the CHM with the terrestrially measured tree height.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table A1 .
Technical details of the Advanced Digital Sensors (ADS) used.

Table A3 .
Gain in completeness with the different strategies and image stripes used.Comparison of the year 2007 where a combination of Pan and RGB was used and the years 2008-2012 where matching was based on CIR images (S11: 1st nadir stripe/"ADS steep", S12: 1st nadir stripe/"ADS flat", S21: 2nd nadir stripe/"ADS steep", S22: 2nd nadir stripe/"ADS flat").

Table A4 .
Vertical agreement measurements of the digital surface model based on stereo measurements distinguishing eight aspect classes.Differences are calculated by subtracting the elevation values of the stereo measurements from the digital surface model values.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table A5 .
Comparison of double measurements of the heights of 441 trees within the Swiss NFI to estimate measurement error.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table A6 .
Agreement of canopy height calculated from the CHM with the terrestrially measured canopy height according to species mixture.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table A7 .
Agreement of canopy height calculated from the CHM with the terrestrially measured canopy height according to stand structure.A threshold value of ±3 × RMSE was applied for outlier exclusion.

Table A8 .
Agreement of canopy height calculated from the CHM with the terrestrially measured canopy height according to stages of stand development.A threshold value of ±3 × RMSE was applied for outlier exclusion.