Find the information such as human life, natural resource,agriculture,forestry, biotechnology, biodiversity, wood and non-wood materials.
Blog List
Thursday, 29 December 2016
Large-scale road detection in forested mountainous areas using airborne topographic lidar data
Published Date
ISPRS Journal of Photogrammetry and Remote Sensing February 2016, Vol.112:23–36,doi:10.1016/j.isprsjprs.2015.12.002
Author
António Ferraz a,b,c
Clément Mallet a,
Nesrine Chehata a,d
aUniversité Paris-Est, IGN, MATIS, France
bNASA – Jet Propulsion Laboratory, USA
cR&D Unit INESC Coimbra, Portugal
dBordeaux INP, EA 4592, France
Received 15 May 2015. Revised 22 October 2015. Accepted 3 December 2015. Available online 29 December 2015.
Abstract In forested mountainous areas, the road location and characterization are invaluable inputs for various purposes such as forest management, wood harvesting industry, wildfire protection and fighting. Airborne topographic lidar has become an established technique to characterize the Earth surface. Lidar provides 3D point clouds allowing for fine reconstruction of ground topography while preserving high frequencies of the relief: fine Digital Terrain Models (DTMs) is the key product. This paper addresses the problem of road detection and characterization in forested environments over large scales (>1000 km2). For that purpose, an efficient pipeline is proposed, which assumes that main forest roads can be modeled as planar elongated features in the road direction with relief variation in orthogonal direction. DTMs are the only input and no complex 3D point cloud processing methods are involved. First, a restricted but carefully designed set of morphological features is defined as input for a supervised Random Forest classification of potential road patches. Then, a graph is built over these candidate regions: vertices are selected using stochastic geometry tools and edges are created in order to fill gaps in the DTM created by vegetation occlusion. The graph is pruned using morphological criteria derived from the input road model. Finally, once the road is located in 2D, its width and slope are retrieved using an object-based image analysis. We demonstrate that our road model is valid for most forest roads and that roads are correctly retrieved (>80%) with few erroneously detected pathways (10–15%) using fully automatic methods. The full pipeline takes less than 2 min per km2and higher planimetric accuracy than 2D existing topographic databases are achieved. Compared to these databases, additional roads can be detected with the ability of lidar sensors to penetrate the understory. In case of very dense vegetation and insufficient relief in the DTM, gaps may exist in the results resulting in local incompleteness (∼15%). Keywords
Lidar
Airborne
Road extraction
Classification
Mountainous areas
Forests
Large scale mapping
1 Introduction
1.1 Motivation
Accurate road location, use and condition in forested environments is an invaluable input for many ecological, economical purposes, and in particular for evaluating adequate spatial indicators related to public policies (Coffin, 2007 and Robinson et al., 2010). It concerns forest management or biodiversity issues, wildfire protection, animal movement, leisure activities etc. In mountainous areas, such knowledge is highly necessary since it can have a negative impact on water quality, habitats and can cause landslides (Sidle and Ziegler, 2012).
In addition, for a cost-effective wood supply, mobilization conditions (both harvesting and accessibility) is a prerequisite for setting up adequate chains for the timber industry. Road location knowledge is coupled with planning models and decision-making tools implemented in geographic information systems (Sačkov et al., 2014): forest harvesting, logging and investments in forest infrastructures can be optimized, while taking current and scheduled accessibility of forest resources into account. For those purposes, both 2D road location and geometrical characteristics (slope, width, radius of curvature) are required in order to evaluate which kind(s) of vehicules can go down a given road. The related information should also be updated in a yearly basis (Gucinski et al., 2001).
Various solutions exist for forest road database (DB) generation and updating: existing DB integration, human photo-interpretation, GPS field surveys, crowdsourcing web-based solutions, and automatic remote sensing data analysis. The four first solutions exhibit significant advantages in terms of road centerline completeness albeit lacks can be noticed under dense canopy cover. However, limitations exist in terms of 2D location accuracy and geometrical characterization: the road width remains most of the time the unique available feature. This is not sufficient since (i) only coarse values may be provided, (ii) smaller and narrower roads (e.g., trails and haul roads) may not be described, and (iii) no solution has been developed for the updating task yet.
In this paper, the issue for large-scale road detection and characterization in forested mountainous areas is addressed. For that purpose, the airborne topographic lidar technology is adapted. This is due to its ability to provide reliable and accurate terrain altitude even under forest cover. It offers a mature and realistic solution since forest massifs are now acquired over large areas in many countries (Maltamo et al., 2014).
1.2 Related work
A large body of literature has addressed the problem of (semi-) automatic road detection for more than 30 years. However, existing studies are mainly focused on urban and peri-urban areas using geospatial optical imagery (Mena and Malpica, 2005, Amo et al., 2006 and Mayer et al., 2006). A majority of techniques was designed for specific data and environment. Main differences stem from the selected road model, which is highly correlated to the spatial resolution of the input images.
For spatial resolutions superior to 1 m, roads are considered as thin elongated objects. Top-Down methods based on a line extraction procedure are privileged (Lacoste et al., 2005). In structured man-made environments, this allows to directly generate networks by solving graph problems (Hinz and Baumgartner, 2003 and Türetken et al., 2013). Top-Down methods are also adopted in peri-urban areas in order to cope with occlusions and missing information issues, such as vegetation and shadows (Rochery et al., 2004). They are often embedded into a hierarchical multi-scale approach (Laptev et al., 2000). Initial detections are frequently refined with a radiometric profile tracking step (Hu et al., 2007) in order to improve the centerline estimation. These latter methods are all the more relevant than they can be integrated into workflows with already available initial road information (Fujimura et al., 2008 and Miao et al., 2013): existing 2D topographic DBs are the main data source.
Very High Resolution images (<1 m) are considered for urban areas, where roads are seen as homogeneous areas in terms of geometry, topology, texture, and color. Thus, Bottom-Up approaches, based on classification methods, are prefered (Grote et al., 2012, Mnih and Hinton, 2012 and Ziems et al., 2012). The most recent techniques allow the integration of contextual information. This aims to limit the false detection rate, correlated to the high spectral variability of road areas (cars, road marks etc.). Similarly to Top-Down approaches, they are often followed by a network generation and vectorization step (Unsalan and Sirmacek, 2012 and Matkan et al., 2014), which enforces additional regularization constraints.
Finally, Top-Down and Bottom-Up strategies are merged both at the feature and the decision levels in the most recent techniques (Chai et al., 2013, Montoya et al., 2014and Poullis, 2014). This permits to benefit from advantages of both approaches.
Airborne topographic lidar data features three main advantages for the design and the performance of road detection methods. First, the ability to acquire ground points beneath vegetation canopy permits to generate very high resolution Digital Terrain Models (DTMs) in forested areas. Accurate height information can be extracted (Hatger, 2005 and Matkan et al., 2014). DTM morphological analysis gives access to terrain geometrical features unseen with standard optical sensors (Rieger et al., 1999and Sherba et al., 2014) and without any occlusion and shadowing issue (Hu et al., 2014) In mountainous environments, this corresponds to breaklines and road edges (Rieger et al., 1999 and White et al., 2010).
Second, raw 3D lidar point clouds can be directly processed in order to retrieve the vertical distribution of the objects above the ground (Hatger, 2005). In forested areas, this indicates the presence of low vegetation and it helps predicting occlusion areas (Lee et al., 2005). Finally, 3D points are accompanied with an intensity information, more and more often extracted from full-waveform data (Djuricic and Hollaus, 2013). This spectral information is helpful for discriminating many surfaces (roads from buildings in (Clode et al., 2007), roads from bare earth in (David et al., 2009 and Azizi et al., 2014)). Nevertheless, it cannot be assumed that road shoulders always exhibit a distinct spectral behavior from road centers for precisely defining the road areas. Furthermore, classification techniques based on the intensity feature cannot be straightforwardly transferred to forested areas since correction and calibration methods are not mature enough.
In this paper, we aim to develop a fast and efficient method that is able to extract roads in mountainous forested environments at large scales. Consequently, main existing recent methods (David et al., 2009, White et al., 2010 and Djuricic and Hollaus, 2013) do not fit to these requirements since they were not tailored for large-scale processing or/and involve too complex data (namely lidar waveforms).
2 Methods
2.1 Overall strategy
Our road detection and characterization strategy is based on the two-dimensional analysis of a lidar-based 1 m DTM. This allows to provide a solution efficient at large scales, without complex 3D processing techniques. We acknowledge that advanced point cloud classification methods may improve the discrimination accuracy of pathways. However, it would be to the detriment of the computing time and the scalability of the approach.
We do not aim to extract a network but road patches, for two main reasons. First, solutions that exist in the literature are efficient but to the detriment of the processing time, in particular since a priori knowledge on mountainous areas cannot be straightforwardly integrated in such solutions. Second, our pipeline is also designed for topographic database updating and refinement, which do not necessarily need to handle road junctions and correct topology.
Here, the DTM generation process (off-ground point filtering and surface reconstruction) is not considered and no DTM evaluation is performed. We rely on the adaptive Triangulated Irregular Network solution (Axelsson, 2000). It is commonly adopted in the literature, implemented in various softwares, and known to be efficient at large scales. DTM quality may vary with the area of interest: artifacts (outliers or coarse surface interpolation) can exist, but their knowledge is not integrated in our process.
Clode et al., 2007 and David et al., 2009 integrate lidar echo intensity information for road discrimination but, it was decided here to discard these values. First, intensity calibration without sufficient metadata, partly stemming from full-waveform data, is very difficult (Wagner, 2010), in particular for large areas with significant slopes and height range. Second, to the best of our knowledge, there is no generic correction method that can deal with vegetation occlusion(Wagner et al., 2008 and Korpela et al., 2012). Thirdly, it cannot be assumed that road shoulders have a spectral or geometrical behavior distinct from the roadway.
Consequently, our method is purely based on geometrical information. The following mountainous forest road model is proposed:
•
A road is a locally planar elongated feature in the road direction with relief variation in the orthogonal direction.
•
A road width ranges between 2 and 12 m, i.e., relief variations can be retrieved in between these two values. This fits well with the existing literature and with the width information that can be retrieved within various existing road DBs.
The model aims to discriminate main forest roads from other planar areas (including logging and U-turn areas) without assuming that sharp breaklines exist for both sides of the road. Our method is tailored to retrieve roads in slopped environments but not in case of flat terrains (see Fig. 1). Smaller roads may not be discriminated as well.
Fig. 1. Examples of mountainous forest roads. (a and b) roads can be detected with our model but probably not the one in (c) since no relief variation exist in the neighborhood of road edges.
Our workflow is decomposed into three main steps (cf. Fig. 2 and Fig. 3):
1.
Focusing step: three DTM-based features are computed and fed into a supervised Random Forest binary classifier. A road mask is then generated at the same spatial resolution as the initial DTM (namely 1 m).
2.
Road graph generation: an extensive graph is built upon this mask in order to propose the largest number of conceivable road candidates. These patches are then pruned following the classification-based criteria.
3.
3D geometric characterization: for each section, width and slope features are computed from the DTM with a joint profile and Object-Based Image Analysis (OBIA).
Fig. 2. Overview of the proposed approach. See text for more details on the three steps.
Fig. 3. Illustration of the method over an area of interest. (a) Orthoimage with existing road DBs. (b) Slope map. (c) Gradient slope map. (d) Road graph after edge decimation. (e and f) OBIA before and after final filtering. Each section is colored with respect to its width (see Section 2.4). : <3.5 m – : 3.5–4.5 m – : 4.5–5.5 m – : 5.5–6.5 m – : >6.5 m.
2.2 Focusing step
2.2.1 Feature extraction
Maps of slope ( ), roughness ( ) and slope gradient ( ) are calculated from the DTM at the same resolution (1 m) according to a 8-connected neighborhood. The slope is defined as:
equation1
where dX and dY are the derivatives of the altitude with respect to a central pixel s in the East-West and South-North directions, respectively. This attribute is used to detect locally flat areas (Fig. 3a and b). Terrain breaklines are better highlighted by roughness attribute. It evaluates the heterogeneity of the terrain within an area around a given pixel s. In practice, it is computed as the difference between the altitude of the central pixel s and its 8 neighborhoods according to the following formula:
equation2
In order to deal with large flat areas that can be misclassified as road patches (Chen et al., 2014), we propose a multi-scale slope attribute map ( ). It analyzes the variability of the slope of an inner neighborhood ( ) compared to an outer one ( ). In such a way, mountain road contours are emphasized (Fig. 3a–c) without the need for sharp breaklines for both road edges (Brzank et al., 2008). Moreover, contrarily to (David et al., 2009 and Djuricic and Hollaus, 2013), it does not rely on either height dispersion or intensity to filter non-road features such as low vegetation. Let be:
equation3
with:
equation4
equation5
is a mean filter applied to the DTM according to neighborhoods of different sizes ( m and m in our forest road model).
2.2.2 Focusing mask
A supervised classification is carried out at the same resolution of the DTM using the three maps described in the former section. The training dataset was manually obtained (Section 3) and contains two classes of pixels: road and non-road. A Random Forest classifier (RF, (Breiman, 2001)) was chosen as being one of the most reliable methods dealing with high dimensional data as well as with data requiring high level of generalization (Caruana et al., 2008 and Criminisi and Shotton, 2013). Additionally, RF is weakly sensitive to both noisy data and outliers. This highly facilitates the manual establishment of the training data. RF provide a measure of the classifier confidence called the RF margin. An unsupervised margin is used here. It corresponds to the difference of votes for the road and the non-road classes, normalized by the total number of votes (100 is this study). The margin ranges between −1 and +1. High positive and negative values correspond to a high confidence of the classifier, for road and non-road classes, respectively. Margin values around 0 indicates a low classifier confidence.
2.3 Graph generation
It is decomposed into three main steps: vertices and edges initialization, followed by edge pruning.
2.3.1 Vertex initialization
Graph nodes are selected using the focusing mask, according to two criteria. (i) Nodes should be located on pixels with the highest confidence values of the roadclassification; (ii) they should be homogeneously spatially distributed in order to efficiently sample the area of interest. The most exhaustive graph as possible should be built without proposing too many nodes. For this purpose, the Point Process (PP) technique, issued for the stochastic geometry theory (Descombes, 2011), is adopted. This allows to solve the interleaved problem of object counting (the optimal number of nodes) and parameter estimation (their 2D location). Such technique has been widely adopted in the literature in order to detect objects in remote sensing data (Lacoste et al., 2005 and Chai et al., 2013). Traditionally, one of several patterns is attached to the object of interest (e.g., a rectangle) leading to the so-called Marked Point Processes. Our road model would lead us to adopt a rectangular pattern (Lacoste et al., 2005) but it would results in significant computing times. Therefore, here, the objects of interest are 2D points whose parameters are their 2D spatial coordinates . The research space is limited to the focusing classification mask. The most favorable ensemble of objects (called configuration, ) among the set of possibles is given by minimizing the following energy U (Descombes, 2011):
equation6
where constitutes the set of neighboring objects in the configuration . is a unary term, evaluating whether the current object is located on a high confidence road area. is a binary interaction term between neighboring objects which limits their proximity. The neighborhood symmetrical relationship and V terms are defined as follows:
equation7
equation8
equation9
is the characteristic function. and are the two parameters defining the neighborhood size and the minimal distance between two vertices, respectively. In practice, we have m and m (cf. Section 5 for discussions).
Details related to the PP simulation and the iterative optimization technique for U can be found in (Descombes, 2011).
2.3.2 Edge generation
Subsequently, each vertex is linked to its eight nearest neighbors. They are sampled according to the eight cardinal directions and for distances inferior to 50 m to the vertex. Such a value is selected since it allows to link vertices lying on the same road section even in case of vegetation occlusion. This also ensures not to aggregate those from distinct pathways. In practice, vertices may have less than eight neighbors. A simple solution would be to compute a Delaunay triangulation on the vertices. However, it is not adapted here since many vertices can be aligned or lying on a circle for curved roads.
2.3.3 Edge decimation
Graph simplification is based on the initial road model and on the hypothesis that edges featuring significant relief variation cannot be considered as road. For this purpose, candidates are weighted and selected according to two of the three initial DTM-based features. Let be the edge linking and vertices. We have:
equation10
equation11
where and are the pixel sets of and (slope and multi-scale slope maps, respectively), intersected by the edge . ( , edge set), is removed from the graph if:
equation12
The training set defined for the RF classification is used for the automatic and threshold computation. In practice, the mean is sufficient to capture the geometric heterogeneity of roads within the area of interest. We have and .
2.4 3D characterization
This step intends to remove terrain features that show characteristics similar to roads (i.e., locally homogeneous flat areas) such as natural flat zones, agriculture terraces, ground patches under dense canopy. Due to the lack of lidar measurements there, the terrain is roughly represented by a horizontal plane. The filtering is based on the 3D characterization of the road patches and consists of two steps.
2.4.1 Profile analysis
The terrain surrounding the road edges, which are supposed to be a good approximation of the road centerline, is analyzed in order to estimate the area that corresponds to the actual roadway. The analysis of every edge allows calculating a 1 m binary map road/non-road ( ) over which an object-oriented approach will be carried out. For that purpose, a region-growing approach is set up: edges are discretized into points that are 1 m apart from each other. Then, starting from the center line, for every point, the slope profile, defined orthogonally to the edge direction, is evaluated. Such a profile is extended so that either its length is less than 12 m or the slope between two consecutive pixels is less than (the latter is evaluated from the previously established training set). Finally, the value 1 is assigned to all the pixels from that intersect such profiles.
2.4.2 Object-based analysis and geometrical characterization
An Object-Based Image Analysis (OBIA) is applied on the road clusters extracted from using a 8-connected component segmentation (Fig. 3e and f). Statistical features are selected from a Structural Feature Set (SFS, (Huang et al., 2007)) that characterizes objects from the extensions of their direction-lines histogram. An extension in a given direction for a given pixel gathers together all the nearest neighborhood pixels that show a lower difference in the spectral value than a certain threshold. Out of the many attributes provided by the SFS, only length and width are extracted. To reduce the calculation time and analogously to the former steps (Sections 22.2 and 22.3), only the 8 main directions are analyzed. Moreover, the maximal length for the direction-lines extensions is limited to 50 m within (Section 2.3).
To eliminate residual errors corresponding to small flat zones underneath the canopy, objects with a ratio length/width less than 10 are removed (Fig. 3e and f). Finally, the slope of the road clusters is directly estimated using the DTM.
2.5 Quantitative evaluation
The evaluation is carried out comparing the detected segments to a vectorized reference network (Section 3). The performance metrics are based on segment lengths Heipke et al. (1997). The matching process is based on the “buffer method”. Every portion of one dataset within a given distance from the other is considered as matched. The matched detected data are denoted as true positive with length TP1 while the unmatched detected data are denoted as false positive with length FP. The second matching step is performed in the other direction leading to TP2 length of reference data that matches the detected road network. The unmatched reference data are denoted as false negative with length FN. Four usual metrics are then derived; the true positive metric TP which is the average of TP1 and TP2, the false positive FP (commission), the true negative TN (omission) which is the length of the reference data that are not included in the detected network. These metrics are normalized by the length of the reference road network and thus range between 0 and 1. Finally the errors sum is denoted ERR = FP + TN.
3 Data
Our method was tested over a 1425 km2 study site located in the Vosges mountains in North-East of France (Fig. 4). The slope ranges from 200 m to 1365 m. The French mapping agency (IGN, Institut National de l’Information Géographique et Forestière) collected the airborne LiDAR data using an Optech 3100EA device. The angle scan was set to , avoiding many irrelevant areas at each side of the lidar strip. The footprint was 0.8 m in order to increase the probability to have ground points. The point density, which ranges from 2 to 4 pt/m2, agrees with the acquisition parameters used in many countries for large-scale forest mapping purposes. Our study site was characterized by nearly 5.6 billion of 3D points that had been divided into tiles of 1 km × 1 km for processing purposes. The corresponding 1 m DTM was interpolated from a Delaunay triangulation, computed using lidar points corresponding to ground measurements. The latter had been filtered from the point cloud using a routine based on the TerraSolid software (Soininen, 2014). As mentioned in Section 2.1, both results and accuracy are not assessed in this work.
Fig. 4. Overview of the area of interest (yellow). (a) Localization in France. (b) Perimeter of the dataset. (c) 5 m lidar DTM. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The training data used for the RF classification consist in twenty polygons manually delineated over a few tiles. They comprise 16,250 pixels, i.e., only % of the study site area. Although its small size, the training data was effective in calculating both an accurate focusing mask and for the estimation of many parameters required by our method.
Two topographic databases are available for validation purposes. They represent roads as lines in two dimensions. The first one has been maintained by IGN whereas the second one by the French National Forest Office (ONF, Office National des Forêts). Both have been obtained by photo-interpretation, extensive surveying methods, and local authorities DB integration. Their completeness is very high even under dense forest cover. However, they are not spatially accurate in terms of 2D location. Therefore, such DBs are not optimal for comparison with our method. As a result, a ground truth consisting in 320 km of roads was manually outlined over nearly 100 km2 using the lidar DTM and Very High Resolution aerial orthoimages. While both the IGN and ONF DBs will be only used to visually assess the results (Section 4.1), the lidar-based DB allows for an accurate quantitative assessment of results (Section 4.2).
4 Results
4.1 Qualitative assessment
The main advantages of our approach are (i) the ability to detect roads that are not present in the existing databases, (ii) the location accuracy (Fig. 5 and Fig. 6), and (iii) the capability to estimate the road width (Fig. 3).
Fig. 5. Incomplete road detection mainly comes from dense forest canopy which may prevent a clear DTM generation with lidar data. (a) Orthoimage (0.5 m resolution). (b) Slope map . (c) Graph before and after decimation (pink and blue, respectively). (d) Final result with width feature (see Fig. 3 for the color code). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 6. Advantages and drawbacks of the proposed method. (a) Orthoimage of the area of interest, superimposed with two existing road DBs (red and green) and the hydrographic network (blue). (b) Slope map : one can notice the completeness of the road network, the current lacks, and the planimetric shifts with the genuine location of the roads. (c) Our result (see Fig. 3 for the color code). One can see the roads are correctly retrieved (location and completeness), that junctions are rarely detected, conversely to parts of the hydrographic network. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Roads that are not present in the existing sources of information vary from the IGN to the ONF DBs. They mainly correspond either to roads underneath dense forest cover or to roads that cross bare earth and grasslands. In both cases, they cannot be easily extracted with the visual interpretation of aerial images. Moreover, these non-inventoried sections correspond to long patches of more than 100 m. As a result, this feature is a clear advantage of our method despite a low rate of missing roads in the existing DBs (5–10%). However, we should note that the Vosges mountains have been much better inventoried than other French regions.
In addition, our method allows adjusting the location of the roads from 5 to 10 m compared to IGN/ONF DBs. The high spatial accuracy of the lidar-based road maps compared with other techniques have already been proved by White et al., 2010 and Azizi et al., 2014. Here, it is assessed at larger scale. Nearly 40% of the roads are poorly located in both DBs, which also shows the relevance of the method for updating and refinement purposes.
Our workflow can be used to locally calibrate the road location or to assign road characteristics (e.g., width, slope) to road patches (skipping step 2, Section 2.1). Fine details such as the angle or the length of a turn would be refined as well (Fig. 6). This can be crucial for the estimation of the radius of curvature of turns, which are needed for the classification of roads patches according to the vehicles that can drive through them.
The method exhibits two types of drawbacks. First, our method is based on the slope gradient. Therefore, it is only adapted to extract forest roads over certain environments. The discrimination between flat surfaces from roads is only based on an elongation ratio, which is not adapted for the extraction of urban roads, U-turn areas and logging areas (left part of Fig. 3). In such areas as well as over large talwegs, an over-detection is observed. In general, human activities, such as plowing, irrigation, agriculture terraces, bound our method for an over-detection. For instance, machinery activities in recent plantation areas form terrain features similar to forest roads. Indeed, most often, roads between plantation rows are kept clean during the first years for plowing purposes. Unfortunately, these misdetections are difficult to identify without both contextual information and stronger geometric assumptions.
Second, under-detections are observed when the road does not follow our model. This corresponds to roads exhibiting homogeneous slopes compared with the surrounding areas. For instance, high forest cover can prevent lidar to well characterize underneath terrain and to properly extract roads. On the other hand, under-detection can be originated by poor graph decimation. For example, if a graph vertex is located slightly off-road, high and (Eq. (11)) will be assigned to such a vertex which will be erroneously removed (Fig. 5). There is no constrain that forces two vertices to be aligned to each other: a graph vertex can be located either on turns or on road edges. We could set to a lower value to better represent road centerlines and turns. However this would generate a much higher number of vertices and, consequently, higher computation times. Similarly, for terrains occluded by dense vegetation, additional vertices under canopy cover would not change the results since the resulting edges would not fulfill our requirements. They would be removed during the 3D characterization step.
Moreover, road junctions are not taken into account either on the PP or on the OBIA step. In case of multiple intersecting roads or a split of a road into two pathways, the road model do not hold. These areas are poorly retrieved and one direction is privileged over the other ones (Fig. 6c).
4.2 Quantitative assessment
320 km of roads were delineated by photo-interpretation and considered as our ground truth. The buffer width was fixed to 5 m in order to take into account the different planimetric accuracies of reference and detected road networks. One can observe on Fig. 7a that the proposed algorithm leads to high True Positive road length that vary from 90% to 82% without length pruning and with a strong length pruning (L> 3 m), respectively. In the latter case, the True Negative length (TN) increases and many road sections are missed (cf. Section 4.1 and Fig. 8). In the contrary, with a low pruning length, the False Positive length (FP) is higher, leading to many false detections essentially due to small road segments. Performance results are similar to those reported in Sherba et al. (2014) in which a semi-automatic workflow at a smaller scale was proposed. No in-depth critical analysis were proposed in Rieger et al., 1999, David et al., 2009 and Djuricic and Hollaus, 2013 whereas White et al. (2010)only focus on centerline location evaluation.
Fig. 7. Impact of the evaluation parameters. : TP – : FP – : TN – : ERR.
Fig. 8. Road detection evaluation: buffer size = 5 m, minimal length = 2 m. Top: orthoimage. Bottom: results. In green, the well-detected reference roads, in blue overdetections and in red, missing reference roads. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 7b shows the impact of the evaluation buffer size on road detection performances. The higher the buffer size, the better the road detection performances. This is due to different planimetric accuracies between both road networks. Finally, a trade-off is reached with a buffer size of 5 m and a pruning of road parts that are shorter than 2 m.
This quantitative evaluation emphasizes the main performance conclusions stated in the previous section. In Fig. 8, one can note the satisfactory detection of the roads. False positive areas are of two types: punctual misclassification of flat areas, resulting in overdetection spots, and lines orthogonal to the real network, corresponding to talwegs. Similarly, the two main sources of errors can be easily retrieved: vegetation occlusion and inability to handle road junctions. Open areas prevail in the Northern part of the area of interest. This explains why only small gaps and junction problems are noticed there. Conversely, one can see in Fig. 8 that the South-East part concentrates most of the longer underdetected roads. This is due to a more significant terrain slope and denser canopy cover, subsequently leading a poorer terrain description. Since our work is uniquely based on a DTM, no additional constraints can help us to fill such gaps.
4.3 Computing times
The computing times of the method steps were evaluated over more than 70 tiles (1 km × 1 km) over the Grand Ballon area, the highest point of the Vosges mountains (1424 m). It allows to ensure a high diversity of landscapes in terms of slope, canopy cover, road density and complexity. Per-step figures are provided in Table 1 (for a single 3.3 GHz core). They correspond to the mean value of all the processed tiles.
Table 1. Mean computing time per km2 for each step of the method.
Step
Process
Time (s)
1
DTM-based feature computation
2″
Slope gradient computation
6″
Supervised RF computation
16″
2
Point Process
22″
Graph initialization and decimation
19″
3
Profile analysis
1″
OBIA and 3D characterization
36″
Linearisation
11″
Shapefile generation
2″
Total
1′55″
The full workflow takes less than 2 min per km2. This is very satisfactory for operational purposes. It stems from the facts that (1) methods have been adopted and/or tailored to be of low computational demand and (2) most of them are based on well-established and efficient existing libraries: Sherwood (RF Classification, (Criminisi and Shotton, 2013)), LibRJMCMC (Point Process (Brédif and Tournaire, 2012)), Orfeo ToolBox (Image processing (Inglada and Christophe, 2009)) and CGAL (Graph analysis and management, (CGAL, 2015)).
The three most time-consuming tasks are the initial pixel-based RF classification, the vertex detection, and the final object-based analysis.
•
The RF processing time could be reduced with a coarse initial image over-segmentation or by adopting a parallel implementation of the RF algorithm (Sharp, 2008).
•
The vertex detection is already highly optimized since the computing time of 22” corresponds to the 21 million iterations of the underlying stochastic process. PP-based techniques are known to feature significant convergence time and therefore may not be designed for large-scale issues. However, the LibRJMCMC implementation and the non-complex energy formulation allow to cope with such an issue.
•
Eventually, most of the computing times of the OBIA step stem from the Structural Feature Set texture feature extraction. However, it remains one of the most relevant attribute for road detection.
5 Discussion
First, one can note that our workflow is highly dependent of the DTM-based feature extraction step, and, in particular, of the slope-gradient image. Such an image is built over the selected road model, and subsequently, over two thresholds that are the minimal and maximal extents of analysis: m and m. These parameters were selected in order to capture the relief variations between two scales while ensuring that the first one corresponds to a road patch and the second one to a non-road area (for road pixels). Such a behavior is illustrated in Fig. 9. With , the function acts as a mean filter and smooths the slope map (Fig. 9a, b, and c). While large road sections could still be retrieved, smaller elements almost disappear. In addition, flat non-road areas are exacerbated which does not help discriminating roads and trails. In Fig. 9d, e, and f, we can see that m allows a correct detection of road edges albeit we do not focus on locating them accurately here. The variation from m to 12 m shows that we are able to correctly focus on areas with changing topography. m offers the best solution since the resulting image is less noisy and better removes the flat areas. Road sections are also enhanced in terms of magnitude. DTM noise can come from artifacts in the initial point cloud filtering process and from real topographic high frequencies. A large value permits to avoid taking the DTM quality into account in the feature extraction step. In addition, this ensures a general applicability of this parametrization without a per-pixel and per-area neighborhood analysis. Large roads, hauls and trails can be detected with the same parameter set.
Fig. 9. Relevance of the road model for focusing on road areas (1 km × 1 km area). (a) Slope map . (b–f) Slope gradient maps computed with various scales of analysis. Image (f) provides the best discrimination and its parameters correspond to the ones selected for the road model.
Second, the performance of our workflow relies on a correct 2D coarse road detection in order to provide a fast-yet-complete focusing mask. The RF classification method was perfectly tailored for that for two main reasons: (i) its very high generalization performance and (ii) its intrinsic confidence measurement of the labeling process. A very limited training dataset was sufficient to accurately label road areas (see Fig. 10). As discussed in Section 4, false positive and false negative detections come from areas with similar topography (talwegs, see Fig. 10b) and unsufficient terrain description, respectively, but not from the classification method. The RF unsupervised margin provides an additional relevant input to the binary discrimination (Fig. 10c and focus on Fig. 11a). For road pixels, lower margin pixels (<0.8) correspond to non-roadareas, high margin pixels are flat areas while road segments always correspond to very high margin values (the margin mean value for road pixels is 0.94). The separation between these clusters did not require any hard threshold but was smoothly integrated into the vertex detection procedure.
Fig. 10. Illustration of the classification results over a complex area(3 km × 1 km). (a) Slope map . (b) Binary RF classification. (c) RF margin. One can see that the margin is helpful for better discriminating road areas from other topographic features following the road model (talwegs and flat areas), with lower margin.
Fig. 11. Illustration of the graph generation procedure. (a) RF margin values for the road class (white: non-road pixels). (b) Edge decimation : initial graph. : graph after slope-based threshold ( ). : final graph after slope-gradient threshold ( ). See Section 2.3 for more details.
One can see in Fig. 11 that areas without high margins are naturally discarded during the graph initialization. Furthermore, the margin constraint coupled with the spatial regularization inserted in the PP allow to regularly sample the terrain. This avoids a redundant number of vertices in non very high margin areas since the vertex density is correlated to the margin value. This explains why only few vertices are retained in the flat area in the center of Fig. 11b. Flat areas and talwegs correspond to less than 3% of the whole area of interest, this is also why no specific filter was proposed to address such issue.
In Fig. 11, one can see that m and m are suitable values for a fine yet simple network restoration. is selected as the upper bound of analysis of the road model (Section 2.1) so as to favor edges fully lying on high margin regions. A lower value is not necessary for better road curve sampling, similarly to a higher cardinality-based neighbor selection. It would lead to additional vertices and edges and subsequently to multiple road centerlines for large road sections. Indeed, no orthogonal nor parallel constraints were inserted into the regularization term of the Point Process. m is not a critical parameter and allows connecting distant vertices along a given road track (Fig. 11b) without generating too many irrelevant edges. In addition, since it controls the genuine number of neighbors, it acts as an implicit weighting parameter for the prior term of the energetic formulation (Eq. (6)) with respect to the data term. Thus, no trial-and-error tuning is required.
In Fig. 11, one can also notice the good performance of the edge decimation procedure. More than 90% of the generated edges are removed using the slope-based and slope-gradient thresholds ( and , respectively). This figure varies with the complexity and the density of the road network. In practice, it ranges from 77% to 93%, depending on the tile. Both thresholds are relevant and act in a coarse-to-fine flavor. First, approximatively 75% of the deleted edges are directly removed with . Small local relief variations are sufficient to discard such edges for further analysis, corresponding both to significant slopes and some road shoulders (center of Fig. 11, where the flat area in red is eliminated due to a fluctuating topography). Second, with a more global reasoning, removes 25% additional edges located in very flat areas (right part of Fig. 11) and in both sides of road segments (left part). The limited over- and undetection rates discussed in Section 4.2 again show that the thresholds statistically learnt for a limited training set are sufficient for a large-scale detection process.
Finally, the geometrical characterization step is illustrated in Fig. 12. The graph generation procedure, despite the removal of a large number of unrealistic edges, still results in many small edges, disconnected from the main network. This leads to small clusters after profile analysis (Fig. 12a and b). The SFS-based width and length features (Fig. 12c and d), synthesized into the elongation attribute (Fig. 12e), are then particularly suited for efficiently removing these erroneous regions (Fig. 12f). Thus, a 50 m eight-direction analysis provides sufficient reliable values. No additional computational effort is required. An extension of the SFS to 20 neighboring directions and to a distance of 100 m would have lead to an increase of 15–20% in terms of computing times without further hints towards erroneous segment removal. Indeed, one can notice in Fig. 12f, g, and h that the proposed ratio value cannot discard all of wrong segments. This is due to the fact that such areas locally fit to the road model (>2 m width) and have a length superior to 20 m (25 m for those of Fig. 12). The ten-time elongation ratio cannot be increased without removing correct road segments which mainly correspond to trails and hauls lying on unsufficiently-described terrain under dense canopy cover (e.g., bottom-left segment in Fig. 12g and h). A future line of research consists in deriving this elongation ratio from slope and canopy-cover density maps so as to locally better match to landscape characteristics.
Fig. 12. 3D characterization step over a 1 km × 1 km area. (a) Orthoimage with existing IGN road and river DBs. (b) Road clusters after profile analysis. (c and d) Object width and length, respectively, computed from the SFS feature. (e) Object elongation. (f) Final road segments after elongation pruning. (g and h) Vectorized output with width and slope attributes, respectively.
6 Conclusions
The proposed method, the first one on a such scale, allowed a fully automatic road detection in forested mountainous areas with high completeness and few erroneous detections. The developed approach is fast. For this purpose, effective classification and pattern recognition tools were adopted. Our methodology was based on open-source libraries. It is easily reproducible and each step can be enhanced separately. The method is composed of two separate steps; the road detection and the geometrical characterization. The second step can also be considered for the recalibration or updating of existing road databases, that are exhaustive but without satisfactory geometric accuracy.
Only minimal and maximal road bounds parameters have to be tuned manually. We showed that a unique set of parameters is sufficient to retrieve the large variety of forest roads. Additionally, the learning step was defined with a very limited number of manually selected pixels. It allowed to tune almost all parameters automatically and leads to a satisfactory trade-off between the number of true positives and false detections.
In this work, we had few hypothesis about the geometric forest road model. However, using the full 3D point cloud should allow a 3D feature extraction, and therefore a more accurate focus on roads rather than on talwegs. It should also allow to prune unrealistic graph edges known the local off-ground (the presence of low vegetation for instance should eliminate the corresponding graph edges). Future work will be on both methodological improvements and operational use. First, the road detection step, based on a bottom-up approach should be enhanced with a more complex top-down approach for the classification (Wegner et al., 2013 and Montoya et al., 2014), for the graph construction (Chai et al., 2013) or by using a priori knowledge on forest cover, on road material, etc. For an operational use, the road detection has to be qualified automatically to provide a reliability map: it highlights false detections and missing roads to a human operator who will manually refine the results. Such a semi-automatic scheme will allow to reduce the initial learning step impact and to highly limit the false detection rate. It will open the door to active and semi-supervised learning techniques for a higher degree of versatility and scalability. Finally, a specific study will be carried out on the impact of the DTM quality on the detection results.
Acknowledgements
The authors would like to thank Nicolas David (IGN), Alain Munoz (ONF), Sylvain Dupire and Jean-Matthieu Monnet (IRSTEA) for fruitful discussions on the topic, Adrien Gressin for software development, Jean-Stéphane Bailly (UMR LISAH) for providing the code of road evaluation, and Mathieu Brédif for his help on the LibRJMCMC library. This work was supported by the French National Research Agency (ANR) through the FORESEE project (ANR-2010-BIOE-008) and by the Portuguese Foundation for Science and Technology under project grant PEst-OE/EEI/UI308/2014.
Caruana, R., Karampatziakis, N., Yessenalina, A., 2008. An empirical evaluation of supervised learning in high dimensions. In: International Conference on Machine Learning. IMLS, Helsinki, Finland, 5–9 July, 2008.
Djuricic, A., Hollaus, M., 2013. Extraction of forest roads from full-waveform airborne laser scanning data. In: EGU. 7-12 April 2013, Vienna, Austria.
Gucinski, H., Furniss, M., Ziemer, R., Brookes, M., 2001. Forest roads: a synthesis of scientific information. Tech. Rep. PNW-GTR-509, USDA Forest Service, Pacific Northwest Research Station: Corvallis, OR, USA.
Sharp, T., 2008. Implementing decision trees and forests on a GPU. In: European Conference on Computer Vision, 12–18 October 2008, Marseille, France. pp. 595–608.
is a mean filter applied to the DTM according to neighborhoods of different sizes ( 
m and 
m in our forest road model).
is the characteristic function. 
and 
are the two parameters defining the neighborhood size and the minimal distance between two vertices, respectively. In practice, we have 
m and 
m (cf. Section 5 for discussions).
: <3.5 m – 
: 3.5–4.5 m – 
: 4.5–5.5 m – 
: 5.5–6.5 m – 
: >6.5 m.
), roughness ( 
) and slope gradient ( 
) are calculated from the DTM at the same resolution (1 m) according to a 8-connected neighborhood. The slope 
) compared to an outer one ( 
). In such a way, mountain road contours are emphasized (Fig. 3a–c) without the need for sharp breaklines for both road edges (Brzank et al., 2008). Moreover, contrarily to (David et al., 2009 and Djuricic and Hollaus, 2013), it does not rely on either height dispersion or intensity to filter non-road features such as low vegetation. Let 
. The research space is limited to the focusing classification mask. The most favorable ensemble of objects (called configuration, 
) among the set of possibles is given by minimizing the following energy U (Descombes, 2011):
constitutes the set of neighboring objects in the configuration 
is a unary term, evaluating whether the current object 
is located on a high confidence road area. 
is a binary interaction term between neighboring objects which limits their proximity. The neighborhood symmetrical relationship 
and V terms are defined as follows:
be the edge linking 
vertices. We have:
and 
are the pixel sets of 
( 
, edge set), 
and 
threshold computation. In practice, the mean is sufficient to capture the geometric heterogeneity of roads within the area of interest. We have 
and 
.
) over which an object-oriented approach will be carried out. For that purpose, a region-growing approach is set up: edges are discretized into points that are 1 m apart from each other. Then, starting from the center line, for every point, the slope profile, defined orthogonally to the edge direction, is evaluated. Such a profile is extended so that either its length is less than 12 m or the slope between two consecutive pixels is less than 
(the latter is evaluated from the previously established training set). Finally, the value 1 is assigned to all the pixels from 
, avoiding many irrelevant areas at each side of the lidar strip. The footprint was 0.8 m in order to increase the probability to have ground points. The point density, which ranges from 2 to 4 pt/m2, agrees with the acquisition parameters used in many countries for large-scale forest mapping purposes. Our study site was characterized by nearly 5.6 billion of 3D points that had been divided into tiles of 1 km × 1 km for processing purposes. The corresponding 1 m DTM was interpolated from a Delaunay triangulation, computed using lidar points corresponding to ground measurements. The latter had been filtered from the point cloud using a routine based on the TerraSolid software (Soininen, 2014). As mentioned in Section 2.1, both results and accuracy are not assessed in this work.
% of the study site area. Although its small size, the training data was effective in calculating both an accurate focusing mask and for the estimation of many parameters required by our method.
and 
(Eq. (11)) will be assigned to such a vertex which will be erroneously removed (Fig. 5). There is no constrain that forces two vertices to be aligned to each other: a graph vertex can be located either on turns or on road edges. We could set 
: TP – 
: FP – 
: TN – 
: ERR.
, the 
m to 12 m shows that we are able to correctly focus on areas with changing topography. 
: initial graph. 
: graph after slope-based threshold ( 
: final graph after slope-gradient threshold ( 
25% additional edges located in very flat areas (right part of Fig. 11) and in both sides of road segments (left part). The limited over- and undetection rates discussed in Section 4.2 again show that the thresholds statistically learnt for a limited training set are sufficient for a large-scale detection process.
No comments:
Post a Comment