Automatic Surface Water Mapping Using Polarimetric SAR Data for Long-Term Change Detection

Abstract

Mapping the distribution and persistence of surface water in a timely fashion has broad value for tracking dynamic events like flooding, and for monitoring the effects of climate and human activities on natural resource values and biodiversity. Traditionally, surface water is mapped from optical imagery using semi-automatic approaches. However, this process is time-consuming and the accuracy of results can vary among image interpreters. In recent years, Synthetic Aperture Radar (SAR) images have been increasingly used. Microwave signals sensitive to water content make SAR systems useful for mapping surface water, saturated soils, and flooded vegetation. In this study, a fully automatic method based on robust stepwise thresholding was developed to map and track the change in the extent of surface water using Polarimetric SAR data. The application of this method in both Radarsat-2 and Sentinel-1 data in central Ontario, Canada demonstrates that the developed robust stepwise thresholding approach could facilitate rapid mapping of open water areas with a promising accuracy of over 95%. In addition, the time-series extent of surface water extracted from May 2008 to August 2016 reveals the dynamic nature of surface inundation, and the trend was consistent with the local precipitation data.

1. Introduction

Surface water bodies are the most significant resources on our planet. Human activities, climate, and other environmental changes have strong impacts on the locations and sizes of open water bodies, which in turn can affect surrounding environments, biodiversity, and human wellbeing [1,2,3,4]. Advanced remote sensing technologies provide an excellent opportunity to obtain spatially extensive and temporally repeated data relevant to mapping surface water [5,6,7,8]. With optical imagery, surface water is usually extracted based on the difference in the spectral reflectance between land and water. Water absorbs most of the radiation in near-infrared (NIR) and mid-infrared (MIR), whereas other cover types, such as vegetation, soil, and human infrastructure have a higher reflectance in these spectral ranges. In addition to the direct use of the surface reflectance in multispectral satellite imagery, several indices are commonly employed, including the normalized difference water index (NDWI) that uses the visible wavelength (0.52–0.60 micrometers) and NIR (0.77–0.90 micrometers) and modified normalized difference water index (MNDWI) that uses the visible and Mid-infrared wavelengths (1.55–1.75 micrometers) [9,10]. Despite the values of such metrics, the acquisition of passive optical data depends on weather conditions and solar illumination, which make it a challenge for them to be used in surface water monitoring. For example, flooding usually happens during the heavy rain, which makes it difficult for the optical sensors to obtain data of the earth surface covered by clouds [11].

As an active sensor, Synthetic Aperture Radar (SAR) can collect images day or night, and its microwaves can penetrate through most clouds, rain fields, water vapor, and aerosol layers. In particular, as water has a unique signature in SAR imagery, a calm water body behaves as a specular reflector with respect to the radar wavelength and thus appears as a low-intensity area in SAR images. This contrasts with the higher intensity of the rougher surrounding terrain, characterized by diffuse scattering [12,13,14,15]. However, the contrast is dependent on factors such as polarization and incidence angle of the emitted microwave radiation. For the current Radarsat-2 sensor with an incident angle ranging from 20 degrees to 54 degrees, calm water is considered to be specular reflection, and the backscatter is close to zero [16].

Several studies have been conducted to map surface water with SAR data through the thresholding techniques utilizing the low return signal behavior of open water. Specifically, methods include manually-defined histogram thresholding [17,18], the Otsu’s algorithm [19], radiometric threshold based on a gamma distribution [20,21], and threshold sets based on statistical analysis of the trainings [22,23]. Those approaches are not fully automatic, and some of them require additional processes to generate high accuracy surface water map, such as supervised classification approaches [22,24,25,26], region growing segmentation [20,23,27,28], and machine learning approaches [29,30,31]. Even with their success, these methods require time and human resources and thus are not operational [18,32]. In addition, some methods are incident angle-dependent, as higher incident angle data are selected for their high contrast between surface water and land [27,33]. Although higher accuracy has been achieved, it is not suitable for long term surface water mapping and emergency flood map production as the amount of the data is reduced. A fully automatic and incident angle independent process is thus much needed to support the timely production of surface water mapping.

In this study, a fully automated algorithm is developed to characterize effectively water bodies using Polarimetric C-band SAR data. Surface water extent was mapped through a robust stepwise automatic thresholding (SAT) approach and validated with a manual digitized surface water inventory. The SAT method was also compared with other commonly used approaches for surface water mapping using NDWI and Otsu’s method [34]. To further explore the potential of mapping surface water at a global scale, the algorithm was tested on both Radarsat-2 and Sentinel-1 SAR data. The algorithm was applied on a time series of Radarsat-2 data to map surface inundation dynamics.

2. Materials and Methods

The study area is located in Ontario, Canada (77°28’21" W, 45°59’10" N)(Figure 1) within the Boreal Shield ecoregion, where the climate is generally continental with long cold winters and short warm summers [35]. The historical weather data show the monthly precipitation in the site ranges between 29 mm and 111 mm with the annual temperatures averaged at approximately 4.5 °C. Figure 2 shows the specific precipitation data in the area. Also, the area is within the western portion of black duck breeding range with the southern area covered by deciduous forest and transitioning to the northern boreal forest. The site was selected due to the availability of existing SAR and waterfowl data for individual wetlands, as needed to support initiatives of resources managers conducting wildlife-habitat assessments.

This study utilized a total of 51 scenes of C-band (5.405 GHz) Radarsat-2 Single look complex (SLC) images acquired from May 2008 to August 2016. The images all were acquired from Fine Quad (FQ) polarized mode, with a resolution of 5.2 m and 7.7 m in slant range and azimuth, respectively. Also, images were collected in both descending (acquisition time around 11:30 Coordinated Universal Time (UTC)) and ascending (acquisition time around 23:15 UTC) orbits, and from varies beam mode, including FQ4, FQ9, FQ14, FQ18, FQ23 and FQ28 to maximize the number of scenes for the temporal analysis. (Detailed information about the images is listed in Appendix A.)

Another C-band imager Sentinel-1A with a swath width of 250-km at fine spatial resolution was explored in this study. For algorithm testing, one scene of Dual-polarized (HH/VH) Sentinel-1 SAR data under interferometric wide-swath (IW) mode acquired on 9 June 2016, was first processed to Level-1 ground range detected (GRD), and then downloaded via python download from National Aeronautics and Space Administration Alaska Satellite Facility (NASA/ASF). The spatial resolution was 20.3 m and 22.6 m in slant range and azimuth, respectively, and the incident angle was between 32.9–43.1 degrees.

A digital elevation model (DEM) was employed to perform the terrain correction on the images, as the images were aligned and corrected for elevation interference to ensure that all SAR images were overlaid with each other and with the cartographic files. The DEM used for this purpose was the South Central Ontario Ortho-photography Project (SCOOP) grid (2 m), generated from the stereo imagery collected between April to May 2013.

To cross-validate the surface water mapped by SAR data, one scene of SPOT-7 imagery acquired on 03 May 2016 was also employed. The multispectral satellite imagery has four spectral bands, blue (450–520 nm), green (530–590 nm), red (625–695 nm), and near-infrared (760–890 nm), along with a spatial resolution of 6.6 m. The imagery has been ortho-rectified to be the same georeferenced as the SAR imagery.

The proposed automatic algorithm consists of three steps: (1) pre-processing of SAR data to generate a backscattering coefficient, georeferenced with a provincial 2m LiDAR derived DEM to correct the relief displacement, and boxcar de-speckle filtering to remove the noise within the data, (2) SAT for generating the water extent map, and (3) post-processing of the image using morphological operations to remove an island within a lake and small misclassified areas (Figure 3). The HH polarization was used in this study as it provides the highest contrast between calming open water and land.

The pixel value of the SAR image was determined by the strength of the radar signal reflected from a unit area on the corresponding location in the scene, the backscatter coefficient β₀ was used in calibration to convert the values from the digital number to the reflectivity of the surface objects. Parameter β₀ was the radar cross-section per surface unit of the target with respect to the local incident angle. All SAR data were then converted from raw units to power units (decibels).

Due to the fact that the speckle effect caused by the coherent radiation used by radar systems was present, the de-speckle filter was employed to remove the salt and pepper noise while preserving edges and textural structures prior to data analysis. In this study, a boxcar filter with a 5-pixel by 5-pixel window was adopted, resulting in a unique valley-hill pattern in the histogram that represented a better distinction between water and non-water bodies, shown in Figure 4. The strategy of detecting the bottom of this valley has been proven to be a good threshold that separates water bodies (low backscattering coefficient) and non-water bodies [36]. Also, this pattern exists in all incidence angle of Radarsat-2 SAR data, not only the higher ones. Thus, the normalization between incident angles is omitted in this study. All the pre-processing procedures were designed as a workflow done automatically with python script utilizing the functions provided by PCI Geomatica.

To automatically determine the threshold value located at the valley (shown in Figure 4b), a step-wise method was developed. Specifically, a set of third-order polynomials was employed to fit the histogram in a manner of moving steps. The reason for this is that the third-order polynomial has the shape that best describes the histogram of backscattering coefficient after boxcar de-speckle, and easier to solve the turning points compared with the higher order polynomials. Also, as the image histogram was rough with notches throughout the curve, using local minimum detection would detect the bottom between any two adjacent notches, and yield multiple results that do not represent the true minimum, i.e., the threshold value. Moreover, due to the roughness of the histogram, using one polynomial to fit the entire histogram would yield a threshold value with significant uncertainty. Therefore, a step-wise polynomial fitting was the best choice to solve the threshold from the bottom of the valley.

The iterative process started at the minimum value of the backscatter coefficient value. First, a partial histogram was extracted. Second, a third-order polynomial was fitted to the selected histogram. Last, the local minimum (turning point) was solved from the first derivative of the polynomial equal to zero. As the third-order polynomial was used, the smaller value of the solution was extracted and recorded as the threshold candidate (Figure 5a). Meanwhile, the difference between the two turning points are calculated, and 1/10 of the result value is stored as the step size. The process continued, and moved with the calculated step size to the next part of the histogram, another polynomial was fitted, and a second threshold candidate was solved (Figure 5b). The iterative process continued until the solved threshold value moved outside the fitting area of the selected histogram (Figure 5c). Finally, the optimal threshold value was obtained from fitting through the pre-determined candidates (Figure 5d). It is worth mentioning that when the slope of the two sides of the valley is different, i.e., very asymmetric shape, the candidates were scattered along with the histogram rather than clustered at the bottom of the valley. Therefore, fitting through all the candidates instead of simply calculating the average assures the optimal result.

After the threshold was determined, each pixel in the SAR image was classified into either water or land, depending on whether its value was smaller or greater than the threshold, respectively. Post-processing was used to refine the classification result. First, morphologic dilation [37] was performed by filling any holes or gaps that appeared as “islands” inside the big water bodies. A size filter was applied to remove any isolated water regions that were smaller than 3 pixels.

To further test the proposed automatic thresholding method on different sensor types, one scene of Sentinel-1 data from the same area was employed. For sentinel-1 data pre-processing, the Science Toolbox Exploitation Platform (SNAP) Toolkit developed by the European Space Agency was used. The workflow was built to generate the HH intensity image from high-resolution Level-1 ground range detected products, and then the image was pre-processed for geocoding, terrain-flatten, and the signal was corrected to sigma naught backscatter coefficients (β₀). The boxcar de-speckle filter with 5 × 5 kernels was used to reduce speckle noise. Terrain flattening and terrain correction was conducted using the recently released Shuttle Radar Topography Mission (STRM) 1 arc-second (approximately 30 m) digital elevation model (DEM) (SRTMGL1).

In addition, the Otsu thresholding method [38], one of the most commonly used automatic image thresholding techniques, was utilized to obtain another surface water map and compared with the one generated with the proposed SAT method. The threshold value was determined through an iterative process that minimizes the within-class variance at the same time it maximizes between-class variances.

To access the binary classification generated by the SAT, both qualitative and quantitative analyses were conducted based on the manually digitized shapefile and 7-m SPOT true color imagery. As the shape and the presence of the wetland can be different between the creation dates, the reference data were visually interpreted using the SPOT imagery.

The last step was to process a total of 51 Radarsat-2 scenes of the overlapping area with the above-mentioned procedure. As the study area has a flat terrain, the displacement introduced by different incident angles was less than one pixel (~7-m) after geo-referencing with the 2-m DEM. Surface water bodies were extracted with the proposed SAT method, and the total area of the delineated surface water extent was calculated and plotted against the precipitation data to conduct the time series analysis.

3. Results

3.1. Surface Water Classification

The surface water map generated from proposed SAT method with Radarsat-2 data, acquired on 5 August 2016, is shown in Figure 6a. For comparison, the corresponding manual digitized polygons are shown in Figure 6b. The majority of surface water bodies were successfully identified and extracted, with the kappa coefficient k = 0.9822. A few linear water streams are missing in the classification results as omission errors, and a few commission errors are evident near the bottom of the site, which appears to be inundated water or flooded vegetation falsely identified as open surface water.

The surface water map generated from Sentinel-1 data with proposed SAT (Figure 7a) is consistent with the manually digitized polygons over large water bodies. Similar to the Radarsat-2 result (Figure 6a), the majority of surface water bodies were successfully identified and extracted with k = 0.9797, while part of the water streams was missing in the classification results as omission errors. On the other hand, the water/land class derived from the same Sentinel-1 data with Otsu thresholding is illustrated in Figure 7b, and the surface water extent was not successfully extracted, with k = −0.7181. Although all water bodies are identified, the majority of the land is falsely classified as water.

3.2. Accuracy Assessment

One scene of a high-resolution SPOT multispectral image was chosen to be within the closest date with the Radarsat-2 data for comparison and accuracy assessment. Figure 8a shows the water bodies derived from SPOT NDWI image using Otsu thresholding and illustrates that large size open water bodies, such as the lakes, have been successfully delineated with k = 0.9586. Also, the linear streams which connect the lakes have been differentiated from the land. However, the built-up land features (top left) are falsely identified as surface water objects. The built-up land features have the same positive values as water features in the NDWI image. This is due to the fact that the index is able to efficiently suppress vegetation and soil to negative values, but fails to separate the built-up land background from the positive value range.

To quantitatively evaluate the classification result, a manually digitized shapefile was used to validate the accuracy of classifications and the results are shown in

Table 1. The matrix of classification accuracy.

Dataset (Method)	False Negative (%)	False Positive (%)	Total Accuracy (%)
Radarsat-2 (SAT)	1.15	0.22	98.62
Sentinel-1 (SAT)	1.27	0.29	99.70
Sentinel-1 (Otsu)	0.29	39.35	60.34
NDWI (Otsu)	3.09	0.17	96.73

. As illustrated in Table 1, the overall accuracy for both Radarsat-2 and Sentinel-1 SAR data exceeds over 95%, and SAR water class has higher total accuracy than the one generated from the NDWI image. The water classification of Sentinel-1 data using Otsu’s method has a total accuracy of 60.34%.

For the accuracy of SAR derived water class, the false identification rate is higher than that derived from NDWI image, indicating more land features may have been falsely identified as water (i.e., commission errors were high).

3.3. Time-Series Analysis

The time-series extent of surface water was extracted using the Radarsat-2 data from May 2008 to August 2016 shows the dynamics of surface water. As shown in Figure 9, the surface water extent generally varies in relation to monthly accumulated precipitation as measured at the local weather station. Decreases in precipitation were associated with a lower area of surface water for several of the observed months (August 2008, May 2014, August 2014, and June 2015). The area of water bodies follows the rate of 3-Day accumulated precipitation reasonably well. In general, the area of water bodies increases after an increase of 3-Day accumulated precipitation and the change in surface water agrees with the change in precipitation. The area of water bodies follows the rate of 3-Day accumulated precipitation reasonably well. In general, the area of water bodies increases after an increase of 3-Day accumulated precipitation and the change in surface water agrees with the change in precipitation.

4. Discussion

For the entire Radarsat-2 dataset (51 scenes), with the average imagery size of 2700 pixels by 3800 pixels, the proposed method computes the series of water maps in 559.26 s on a computer equipped with Intel^® i7-6700 3.4 GHz and 16 GB RAM. Among the whole processing time, SAT thresholding takes 248.92 s, post-processing takes 12.95 s, and writing the image file takes 297.39 s.

The proposed SAT thresholding method successfully extracted the majority of surface water bodies. This can also be observed from the histogram (Figure 10), as the proposed method detects the threshold close to the bottom of the valley (indicated by the blue line), while the Otsu’s method locates the threshold near the peak (indicated by the red line). It can be observed that the surface water bodies with the flooded areas only cover a minor portion in the study area, producing a different histogram from the water dominated studies as shown in Martinis et al., 2009 [27], and the Otsu method failed to properly classify the surface water from the SAR data within our study area.

It is worth mentioning that Otsu thresholding is more effective to separate water and non-water classes when the histogram exhibited as a bi-modal histogram, i.e., two distinct objects presented as one maximum for water and one maximum for land. This can be observed in the NDWI image, and Otsu method is able to separate the two classes with an accuracy of 96%. On the other hand, the histogram of boxcar filtered SAR data only has a single modal presented, as shown in Figure 10, using Otsu yield low classification accuracy, around 40%. To better align with the literature, the de-speckle filter employed in Otsu thresholding for SAR data is Lee filter, instead of the boxcar filter used by the proposed method, and the classification accuracy raised to a reasonable 60%.

A variety of factors might complicate the derivation of water surfaces from SAR data. The low backscattering from some bare agriculture fields and radar shadows led to the overestimation of the water extent. As well, some inundated water resulting from heavy precipitation appeared in the SAR data, but was not evident as water in the reference data. For the water surface derived from the NDWI image, the higher commission errors than the omission errors are possibly due to flooded vegetation being classified as water. For the water surface derived from the SAR data, omission errors were observed from smooth objects like pavement and bare soil. Future studies can include additional road masks to reduce the overestimation introduced by roads in urban areas. Omission errors were observed from flooded vegetation close to rivers and lakes (left side of Figure 5b). This finding is consistent with previous studies that have found fully Polarimetric data were preferred for detecting flooded vegetation [25,39]. Furthermore, the Radarsat-2 detected inundated water areas unidentified in the provincial inventory are typically smaller water bodies. Those areas can be a potential breeding area for water birds.

Notably, part of the commission errors may also result from the difference in acquisition time between the spectral image and the SAR image. The reference manual digitized shapefile was generated from interpretation from the recent spectral image, and the actual boundary observed from SAR data did not precisely follow the same edges from the spectral image. Also, the influence of roughening of the water due to wind on the radar backscatter results in diminished water/land contrast, which can affect the classification accuracy. Future study will cooperate with multi-temporal filters [40] and non-local speckle filters [41] to potentially reduce the speckle-noise while the time-sensitive edges are adapted.

5. Conclusions

An automatic surface water extraction algorithm was developed in this study. The result of developed fully automatic SAT method demonstrates the ability to map the dynamics of surface water with promising accuracy. The preprocessing of SAR data and classification is fully automatic utilizing the Geomaticas and the SNAP toolbox and Python packages. The whole process is fully automatic, and the default values that have been encoded are found to be efficient in extracting surface water extent for different scenes. Our findings demonstrate that using an SAT technique can automatically generate the highly accurate surface water extent from C-band SAR data (Radarsat-2 and Sentinel-1) in an efficient and timely manner. Based on the accuracy assessment, 98% of the water bodies have been successfully extracted, and the surface water extent closely follows the boundaries of the main water bodies, such as the lakes and rivers. Omission errors observed from small, linear targets, i.e., streams and river, and commission errors are identified from roads and inundated vegetation. Future study can include local road masks to reduce the overestimation, and test the effects of non-local speckle filter for minimizing the influence of roughening of the water due to wind. Also, it worth mentioning that the proposed method has only been applied to flat terrain, and additional assessments should consider outcomes for rough terrain. Another valuable finding of this study is that time-series results reveal seasonal dynamics of inundation, and show a strong relationship of the response in surface water extent to the changes in precipitation.

Using SAR data thus permits a detailed timely fashioned surface water map for wetland habitat monitoring that was previously missed from inventory data and land classification map. The fully automatic algorithm developed in this study can be implemented in an operational system to generate provincial or national scale, and benefit the long-term water monitoring frameworks for wetland and management. Its ability to track the changes on the surface water will facilitate the habitat assessments for biodiversity dependent on wetlands, such as the effect of temporal variation in wetland habitat availability on breeding American black duck populations and standardized monitoring of change in a wetland habitat.