Sinkholes in Winona
County, MN have been mapped four times since 1985 using
different techniques including field observations, topographic
maps, air photos and Global Positioning System (GPS)
measurements. As of early 2009, these efforts had identified
and inventoried 672 sinkholes in Winona County that are
recorded in the Minnesota Karst Feature Database (KFDB) (See
the KFDB at: http://deli.dnr.state.mn.us/). The acquisition of
one-meter resolution Light Detection and Ranging (LiDAR) images
has significantly increased the speed and accuracy of sinkhole
mapping. One meter shaded relief LiDAR Digital Elevation Models
(DEMs) for Winona County were visually scanned to compare
sinkhole locations in the KFDB with the LiDAR images and to
find new sinkholes in the LiDAR DEMs. The results of this
method indicate that the number of actual sinkholes in Winona
County could be as many as four times more sinkholes than
identified by the pre-LiDAR surveys. To automate sinkhole
detection from LiDAR data at a regional scale, an algorithm was
developed in MATLAB based on image processing techniques. The
algorithm has three steps. The first part detects potential
sinkhole locations as depressions in the DEM using a
morphological operation (erosion). The second part of the
algorithm delineates sinkhole boundaries by automatically
fitting an active contour (snake) around the potential sinkhole
locations. In the last step, a pruning process, based on the
relationship between depth and area of depressions, was applied
to discard shallow depressions. The proposed method was
evaluated on selected parts of Winona County. Evaluations of
precision and recall returned positive results at 82% and 91%
levels, respectively, which are sufficiently accurate to permit
regional-scale, reconnaissance sinkhole mapping in complex
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2 Introduction quality and quantity in underlying carbonate aquifers, as part of the hydrological cycle. Sinkholes have become convenient (but inadequate) indicators of the presence of resource management decisions by regulators. Complete, accurate inventories of sinkholes are therefore needed, but Various techniques and methods are used to map sinkholes including topographic maps, air photo interpretation, methods at a regional scale. For example, depending on the contour interval (map scale) on topographic maps, sinkholes under forest often cannot be seen on the aerial (elevation) resolution of DEMs derived from LiDAR completeness of sinkhole mapping at the regional scale. resolution of LiDAR data is to create hillshade images visually scan the hillshade image at varying resolutions to identify sinkholes. They can also be compared to air Bing Maps. Although visually scanning is simple and for large regions. Also, sinkhole characteristics like area, perimeter and depth can only be measured or determined manually using visual techniques, which is very tedious and can be prone to accuracy problems. An automated method to locate and measure sinkholes from Abstract Sinkholes in Winona County, MN have been mapped four times since 1985 using different techniques including 672 sinkholes in Winona County that are recorded in the Minnesota Karst Feature Database (KFDB) (See the and accuracy of sinkhole mapping. One meter shaded relief LiDAR Digital Elevation Models (DEMs) for Winona County were visually scanned to compare sinkhole locations in the KFDB with the LiDAR The results of this method indicate that the number of actual sinkholes in Winona County could be as many LiDAR surveys. To automate sinkhole detection from LiDAR data at a regional scale, an algorithm was developed in potential sinkhole locations as depressions in the DEM using a morphological operation (erosion). The second part of the algorithm delineates sinkhole boundaries by potential sinkhole locations. In the last step, a pruning process, based on the relationship between depth and area of depressions, was applied to discard shallow depressions. The proposed method was evaluated on selected parts of Winona County. Evaluations of precision and recall returned positive results at 82% and mapping in complex landscapes. LOCATING SINKHOLES IN LIDAR COVERAGE OF A Mina Rahimi Water Resources Science, University of Minnesota, 173 McNeal Hall, 1985 Buford Ave., St. Paul, MN 55108, USA, email@example.com E. Calvin Alexander, Jr. Earth Sciences, University of Minnesota, 310 Pillsbury Dr. SE, Minneapolis, MN 55455, USA, firstname.lastname@example.org 469
NCKRI SYMPOSIUM 2 13TH SINKHOLE CONFERENCE Paleozoic carbonates and siliciclastics. As shown in containing sandy dolomite and quartz sandstone, forms a karst plateau across much of Winona County. Most found in the areas where the sedimentary cover bedrock surface is less than 15 m (50 ft) thick (Figure 2). The mapping of sinkholes in Winona County in southeastern Minnesota began in the early 1980s. sinkholes in Winona County as part of the Minnesota Atlas of Winona County (Balaban and Olsen, 1984). interpretation. The sinkhole locations were compiled surveyed sinkholes to update the sinkhole database in Filin and Baruch (2010) proposed a method to automatically detect sinkholes and associated characteristics on a large scale. They detected the inner where Z is the elevation from LiDAR DEM data. Then, they applied the active contour method (Kass et al., 1988) to delineate sinkhole boundaries. They used the sinkholes for comparing the relative depth of inner shallow depressions. Study Area Winona County in southeastern Minnesota is part of the 1985). Karst lands in Minnesota are developed in 470 = Figure 1. Bedrock geology and distribution of sinkholes in Winona County.
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2 LiDAR DEM images; and 2) to map new sinkholes using the LiDAR DEM images. The goals of the second method are: 1) to apply an algorithm to identify sinkholes automatically in some parts of Winona County; 2) to characteristics like depths, areas and perimeters; 4) to prune depressions which may not be true sinkholes from the list; and 5) to compare the results from processing the algorithm with the visually scanned datasets in the KFDB in order to evaluate the accuracy of the algorithm. Methods Visual Scanning of LiDAR DEMs November 18 and November 28, 2008. The vertical accuracy is 0.161 m root mean square error (RMSE) at a derived from LiDAR were visually scanned at varying resolutions to identify sinkholes. As many as possible of the sinkholes in the early 2009 KFDB dataset have and Alexander (2002) mapped additional sinkholes in Data Base (KFDB) for Southeastern Minnesota in a includes sinkholes, springs, seeps, sinking streams and sinkholes to the KFDB, and 672 sinkholes had been inventoried in Winona County by early 2009. This paper presents and compares two different methods to map sinkholes: 1) to visually scan LiDAR DEM images and 2) to develop an algorithm to automatically detect, delineate, characterize and validate potential compare sinkhole distribution in Winona County that had been mapped during previous decades with the new 471 Figure 2. Map of Minnesota karst lands.
NCKRI SYMPOSIUM 2 13TH SINKHOLE CONFERENCE 472 is compared with the highlighted cells that are covered by the kernel window. After that, the minimum value of next cell and this procedure continues until it reaches are assigned 1 and those which have different values sinkhole has the value of 1 while its surroundings have 0. Thus, the minimum point of the depression is located. is too small, only a few cells are contributed and many local number of cells included, the local minimum calculation Sinkhole depressions have various sizes and shapes, and they can sometimes be compound sinkholes: smaller sinkholes within a larger closed depression. Thus, to locate all of these depressions different sizes of kernel windows are needed; small kernel windows are optimal for small depressions and larger windows are better for larger depressions. Figure 4 shows the impact of the kernel window size on the number of seed points detected in LiDAR DEMs. Comparing kernel windows of 25 with 55 pixels illustrates that small depressions are detected with kernel size 25 while they are missed by kernel size 55. sinkholes locations. In this process, additional sinkholes that were previously missed and new sinkholes which and mapped. dates for some locations. Birds eye view feature resolution, pictometric photos from several directions for particular locations. Both types of coverage can be used visually to inspect the locations of sinkholes. Erosion and Active Contour Algorithm To automatically detect sinkholes and their boundaries, image processing techniques. This algorithm has several steps: 1) detect local minimum points (seed points); 2) delineate depression outlines around each seed point; potential sinkhole; and 4) prune the list of potential sinkholes to differentiate sinkholes from shallow depressions that may not be true sinkholes. Finally, the remaining potential sinkholes were tested for validity checked and entered into the KFDB. points or the lowest point of depressions in LiDAR through their geometric characterization using a morphological tool in MATLAB called erosion. This tool processes images based on their shape. It compares the value of each pixel in the input image with its neighbors and assigns the value on a corresponding cell in the output image. The morphological operation uses neighbors. It can be a matrix with any size. The erosion operation compares the cell value with its neighbors in the kernel window and returns the minimum value in it for that cell in the output image. kernel window as it moves across the input image. As Figure 3. The procedure of erosion function to find local minimum points.
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2 473 process may not converge to the actual boundary of the depression in many cases. In the second step of the algorithm, an active contour, a was used to identify the depression boundary of each of the seed points. The boundary is a closed curve that is gradient in the surrounding region around the seed point (Kass et al., 1988). The gradient is directly derived from the elevation map shown in Figure 5 (top Figure). The magnitude of the gradient corresponds to the slope of the depression (i.e. white cells in the edge map, gradient map, shows the curve around the seed point passing through cells, each with a maximum gradient corresponding to maximum slope. This method, however, is known to be sensitive to initial conditions, such as initial radius, and the Figure 5. In the EdgeMap, top figure, white cells correspond to the maximum slope of a sinkhole. In the bottom, the green vectors are determined by Gradient Vector Flow. These vectors point toward the edge of the sinkhole boundary where there is maximum slope. The red contours show initialization and iterative processing until the contours converge to the sinkhole boundary. Figure 4. Effect of kernel window size on detecting seed points on LiDAR DEMs.
NCKRI SYMPOSIUM 2 13TH SINKHOLE CONFERENCE 474 counting the number of pixels which are located inside the perimeter. Another parameter, depth, is determined by subtracting the median of all of the pixel values along the perimeter from the pixel value of the seed point. The formula is as follows: where z is the elevation value derived from LiDAR data. Pruning DEMs. Filin and Baruch (2010) suggest different validity tests to separate local and shallow depressions from true sinkholes. One test is compactness. For example, as sinkholes often follow a circular shape, the only candidates accepted as sinkholes are those contour lines whose compactness is nearly 1 (i.e., close to a circle). Winona County sinkholes due to their irregular shapes. if the compactness test is used in Winona County. So, another method is required to prune these shallow depressions. Winona County which contains the most representative topography and sinkhole shapes was selected. In the selected area, typical sinkholes were manually area and depth. For each sinkhole, the perimeter was marked by drawing a polygon. Based on the polygon, the area of the sinkhole was calculated. Then, the depth of the sinkhole was obtained by subtracting the elevation of the deepest point within the polygon and This training dataset was used to identify extreme sinkholes in terms of their size and depth. Two types of of at least 90%, compared to the depth of the shallowest Using these two extreme types of sinkholes (see Figure passing a line through the two extremes. In the pruning boundary. Figure 5 shows an example of an EdgeMap, depression boundary where the slope is maximized. Also, it is shown that in homogeneous regions where the where a provided energy function is minimized. In this application, the energy function is (partially) chosen the slope is at its maximum. As presented below, two from the curve itself and b) external forces extracted (Eq.1) Where E is the total energy function, the external energy the interal energy function, controls the behavior of the of internal energy) determine the tension and rigidity of the curve. The tension parameters control how much force is exerted on the contour to make it smaller. The rigidity parameter controls the smoothness and bending of the active contour is solved iteratively, and therefore it needs initialization. The bottom image in Figure 5 shows the iterative process to delineate a sinkhole boundary. Also active contour function, compared with the EdgeMap. Sinkhole Characterizations and perimeter can be calculated for each individual depression automatically. To calculate the perimeter, the distance formula is used: (Eq. 2) where P is perimeter for the individual depression, n is the number of boundary points, and and are the coordinates of the boundary points. As the LiDAR data = ( ) + ( ) + ( ) = ( )+ ( ) =
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2 475 To evaluate this threshold, a smooth region with no sinkholes was selected and the algorithm was run (Figure 7). As expected, many depressions were were eliminated. The three remaining depressions, dams. This example clearly shows that the threshold works well. Results and Discussion Visual Scanning of LiDAR DEMs The previous mapping of Winona County sinkholes had recorded 672 in the KFDB through 2009. Table 2 compares the Winona County sinkhole data in the 2009 KFDB and the results of visual scanning of the Winona County LiDAR data set. The data produced four distinct groupings. as their KFDB locations. These sinkholes served as step, those candidates whose depth vs. area falls below In order to increase the recall rate (Table 1), (possibly at the expense of decreasing the Precision rate), the all the extreme sinkholes (in the training dataset) are located above the test line. This is particularly important for sinkholes with small areas that are in the early stages of development, and thus their shallow depth may place them below the test line. To accommodate decreased by 0.1 meters for sinkholes whose areas were the test line. But, in the pruning, sinkholes with areas 0.26 meters will not be included in the inventory. Figure 6. Two extreme sinkholes in terms of depth and area are identified. The first extreme, in the lower left, are the sinkholes with depths of at least 90% of the shallowest sinkhole. The second extreme, in the lower right, are sinkholes with depth-to-area ratios of at least 90% of sinkholes with the smallest depth-to-area ratio.
NCKRI SYMPOSIUM 2 13TH SINKHOLE CONFERENCE 476 important learning tools. They helped to illustrate what Winona County sinkholes look like in LiDAR DEMs in terms of shape and size. DEMs, but at slightly different locations than were recorded in the KFDB. The difference in locations was attributed to location errors in the KFDB. The old data was explicitly known to have location errors up to hundreds of meters. LiDAR allowed determination of more accurate locations for those sinkholes and to quantify the location uncertainty in the earlier data. to 180 meters. Most of the location corrections were Table 2. Comparison of the original KFDB with the LiDAR Sinkhole Data. Figure 7. Yellow points are identified as shallow depressions by pruning so they are removed from inventory. These points are located below the test line in Figure 6. Points in blue are depressions near the road. They are removed by the buffer tool in ArcGIS. The red points are False Positive (FP) points. They are located above the test line but they are not sinkholes. They are ponds behind dams or in ditches. True Positive (TP): Corrected results False Positive (FP): Unexpected results False Negative (FN): Missing results = + = + Table 1. The definition of recall and precision. Sinkholes in 2009 KFDB Sinkholes visible in LiDAR but not in KDFB Visible in the LiDAR DEMs Sinkholes not visible in LiDAR DEMs Location not Location 66 168 672 651
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2 477 of the entire county. In Winona County 651 potential new sinkholes, not listed in the KFBD, have now been mapped, as shown in Figure 1. Field checks are necessary to verify which LiDAR features are sinkholes and which are other features. If all of these features are sinkholes, they will nearly double the number of mapped sinkholes each currently open sinkhole holds, then Winona County may have up to four times as many sinkholes as are listed in the KFDB, based on visual mapping. Erosion and Active Contour Algorithm A small region of southwestern of Winona County, Minnesota was selected to evaluate the best parameters for the active contour method including examining the initial radius. As mentioned in the method section, the active contour is solved iteratively and then it needs boundary around the seed point. As seen in Figure 10, the sizes and depths of depressions range from very small ones with depths of less than 0.21 meter to very large ones with depths of 1.5 meter and greater. With this variety of sizes and shapes, it is impossible to identify all of the depressions with only one parameter. Therefore, different sets of parameters were examined and three of them were selected. The the largest initial radius (15 m) for the active contour. This parameter set detects large and deep depressions. The second set, with the same kernel window size but shallower depressions. The third set with the smallest kernel window size and initial radius (5 m) is able were not visible on the hillshade derived from LiDAR not visible on LiDAR, can be seen on aerial images. to the surrounding, visible soil moisture contrasts are detectable on aerial images under the right moisture stress conditions. As illustrated in Figure 9, sinkhole D0018 and is seen on the Bing map but it is not visible on LiDAR. DEMs facilitates the mapping of sinkholes with high accuracy and precision. The LiDAR covers the entire region, including many areas previously unsearched Figure 9. Comparison of sinkholes that is visible on Bing map and on LiDAR. Figure 8. Histogram of the relocation distance.
NCKRI SYMPOSIUM 2 13TH SINKHOLE CONFERENCE 478 of shallow depressions; however, some of them have remained. The remaining points after pruning are not true sinkholes; they are ponds behind dams, depressions in ditches, local depressions in quarries and points near stream beds or roads (Figure 12). Note that points near are not counted in calculating precision and recall. As seen in Figure 7, most of the depressions are shallow local depressions (less than 0.15 meter depth). Such shallow depressions are farmed across and are typically to identify very small and shallow depressions. shallow depressions. The problem of small depressions is that gradient changes are very smooth so that they cannot be easily words, there is no sharp transition from the minimum point of depressions toward their surroundings. are more distinct, so the active contour method can identify the boundary for more shallow and small depressions. Figure 11 illustrates an example of three parameter sets for the active contour function. As depression boundaries better than parameter sets 1 and 2. This example clearly shows how the larger initial radius produces a better match with the depression boundary where the depression is large. Validity test To assess the precision of these methods, including erosion, the active contour, and this threshold Winona County with sinkholes of various sizes and shapes were selected and the procedures were run. Figure 10. Sinkholes in Winona County have various size and shape range from very shallow and small to very large and deep. Figure 11. Three different parameter sets for active contour. The red has the initial radius (5 m), the blue has the initial radius (10 m) and the yellow has the initial radius (15 m).
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2 479 0.46 meter: so it plots below the test line in Figure 6 and is are discarded by this pruning is very low compare to the After pruning, the results show out of 127 initial These are called true positives (TP). Of that sample, 21 of them are false positives (FP), which mean they are not sinkholes but have remained after pruning. quarries. The remaining 9 sinkholes were not detected by these methods or were discarded by pruning. These are called false negatives (FN). Consequently, the precision and recall results were calculated for the algorithm. The precision for the selected region in southwestern Winona County is 82%. This means that 82% of the detected sinkholes are true sinkholes, and the remainders are false positives. The recall is 91%, which indicates this method only misses 9% of sinkholes. Considering the heterogeneity of Winona County (complex topography, woods, quarries, natural method works wells to detect sinkholes. This automatic increase the precision and recall. to isolate these subtle local depressions was needed. Although pruning discards most of shallow depressions, true sinkhole may also be removed. Based on the threshold, Figure 13. Sinkhole with an area of 1800 square meters and a depth of 0.46 meters is eliminated through pruning. Figure 12. Local depressions have remained after pruning that called False Positive (FP). They are located in ditches and quarries.
NCKRI SYMPOSIUM 2 13TH SINKHOLE CONFERENCE 480 References System Research Institute. Available from: Survey. 8 p. Bing Map available at: Dalgleish JB. 1985. Sinkhole distribution in Winona County, Minnesota [master's thesis] Minneapolis (MN): University of Minnesota. 95 p. + map. Filin S, Baruch A. 2010. Detection of sinkhole hazards using airborne laser scanning data. Photogrammetric development of a karst feature database for southeastern Minnesota. Journal of Cave and Karst of the Upper Mississippi Valley Region. Studies in Kass M, Witkin A, Terzopoulous, D. 1988. Snakes: active contour models. International Journal of Magdalene S, Alexander EC Jr. 1995. Sinkhole distribution in Winona County, Minnesota Engineering and Environmental Problems in Karst Terrane. Proceedings of the Fifth Multidisciplinary Conference on Sinkholes and the Engineering and Conclusions The advent of high resolution LiDAR DEMs facilitates accurate and thorough sinkhole mapping. In the visual scanning process, comparing LiDAR data with KFDB locations which are the same as LiDAR locations; KFDB sinkhole locations slightly different from LiDAR data; KFDB sinkhole locations that are not visible on LiDAR; and additional sinkholes which are not listed in the KFDB. Comparison of these two data sets indicates that Winona County probably contains up to four times as many additional sinkholes as are indicated in the KFDB. mapping, an algorithm was developed to detect sinkholes automatically. To assess this method, selected regions in southwestern Winona County were analyzed. First, seed points on LiDAR DEMs. Then, the active contour method was applied to identify depression boundaries based on seed points. Next, the list of potential sinkholes was characterized. Finally, a threshold was set, using the relationship between area and depth, to distinguish true sinkholes from other local depressions. After this pruning, the precision shows that 82% of detected sinkholes are true sinkholes and the remainders are false positives, located along natural watercourses, ditches or roads or in of sinkholes correctly, and misses only 9% of sinkholes Considering the region to which the method was applied, with a variety of features (such as wetlands, reasonable to map sinkholes. In future work, this method will be applied for other areas of Winona County, the results will be compared with the KFDB, the LiDAR DEMs will be visually scanned and Acknowledgments Environment and Natural Resources Trust Fund as on Minnesota Resources (LCCMR).
13TH SINKHOLE CONFERENCE NCKRI SYMPOSIUM 2