Statistical Analysis for Tidal Flat Classification and Topography Using Multitemporal SAR Backscattering Coefficients

Abstract

Coastal zones are very dynamic natural systems that experience short-term and long-term morphological changes. Their highly dynamic behavior requires frequent monitoring. Tidal flat topography for a large spatial coverage has been generated mainly by the waterline extraction method from multitemporal remote sensing observations. Despite the efficiency and robustness of the waterline extraction method, the waterline-based digital elevation model (DEM) is limited to representing small scale topographic features, such as localized tidal tributaries. Tidal flats show a rapid increase in SAR backscattering coefficients when the tide height is lower than the tidal flat topography compared to when the tidal flat is covered by water. This leads to a tidal flat with a distinct statistical behavior on the temporal variability of our multitemporal SAR backscattering coefficients. Therefore, this study aims to suggest a new method that can overcome the constraints of the waterline-based method by using a pixel-based DEM generation algorithm. Jenks Natural Break (JNB) optimization was applied to distinguish the tidal flat from land and ocean using multitemporal Senitnel-1 SAR data for the years 2014–2020. We also implemented a logistic model to characterize the temporal evolution of the SAR backscattering coefficients along with the tide heights and estimated intertidal topography. The Sentinel-1 DEM from the JNB classification and logistic function was evaluated by an airborne Lidar DEM. Our pixel-based DEM outperformed the waterline-based Landsat DEM. This study demonstrates that our statistical approach to intertidal classification and topography serves to monitor the near real-time spatiotemporal distribution changes of tidal flats through continuous and stable SAR data collection on local and regional scales.

1. Introduction

Tidal flats are extensive coastal ecosystems that form intertidal areas and are periodically inundated due to the ebb and flow of tides, which cover an area of 128 K km² of the Earth’s surface [1,2]. They play a vital role in removing inorganic and organic matter from the land, stabilizing the coastline, and providing protection against coastal flooding caused by hurricanes and storm surges [3,4,5]. Despite the geographical and ecological importance of tidal flats, the global area of tidal flats has decreased by 16% over the past 30 years [1,6,7]. In particular, the Korean coastal environment has been undergoing rapid changes due to coastal reclamation and development, thus the total area of tidal flats in Korea has decreased by 30–40% compared to the 1980s [8]. Tidal flats are constantly changed by waves, wind, and biological and anthropological effects. Therefore, when considering the physical, biological, geological, and morphological aspects of tidal flats, it is very important to accurately investigate extent and exposure of the tidal flats along with the tide level and the spatiotemporal change of the tidal flat topography [1,9].

Ground-based point measurements provide information to improve our understanding of complex tidal flat environments on a local scale, despite involving expensive, labor intensive, and time consuming processes. Recently, the enhanced operational stability of unmanned aerial vehicle (UAV) platforms has increased the use of rapid and automated intertidal mapping. The high spatial resolution of UAV images can be very useful in detecting small changes in tidal flats. Long et al. [10] attempted topographic monitoring over tidal flats using a fixed-wing mapping drone combined with a digital camera and showed high accuracy in the result data. Brunetta et al. [11] presented a method to detect detailed topographic changes in tidal flats with various ground control point configurations. However, this method requires on-site visits to obtain the data and the operating environment is readily affected by weather conditions, such as extreme wind, rain, and storms. Similar to in situ measurement, UAV platforms are challenging to use for the characterization of large scale tidal flat environments due to their limited flight time and coverage.

Several studies have demonstrated that satellite-based remote sensing represents an effective tool for monitoring morphologic change and sedimentary facies distribution over wide tidal flats, on regional and continental scales [9,12]. Specifically, the Normalized Difference Water Index (NDWI) from Landsat TM/ETM+ data was applied for the large scale intertidal classification and coastline extraction of East Asia [8] and Australia [13]. However, optical remote sensing imageries are always limited by weather conditions such as sunlight and cloud coverage. Synthetic aperture radar (SAR) images have also been applied to tidal flats with less cloud contamination [14], but SARs are readily affected by high speckle noise. The contrast between ocean and land is not always clear on a pixel-by-pixel basis in high resolution SAR images, which makes the waterline extraction of SAR images a challenging problem [15]. Thus, the integration of multitemporal and multiband SAR data can be very useful for the study of tidal flats. For example, Zhang et al. [16] used multiband SAR images, including C-band ENVISAT ASAR and L-band ALOS PALSAR, to analyze the characteristics of seasonal changes in reed marches and rice fields. Salameh et al. [17] showed that the frequent sampling of the Sentinel constellation helps to overcome the limitations induced by a shortage of images. The combined use of Sentinel-1 and -2 can provide an adequate number of images that covers the whole tidal cycle in a relatively short period.

Intertidal topography has been generated through waterline extraction using multitemporal remote sensing data. Optical satellite images have been largely used to generate topographic maps [7,18], and SAR images have also been moderately applied in the waterline extraction method [19]. This uses a time series of satellite images taken during different tidal periods to map intertidal topography. The heights of the waterlines can be measured from in situ data or marine hydraulic models, and a gridded DEM is created using spatial interpolation techniques [20]. Although waterline extraction is simple and fast, it is sensitive to noise and the edge location is not accurate with only a limited number of remote sensing images [15,21]. Although this method is still commonly used, DEM generation is constrained between the highest and lowest level of acquired images and provides contour lines without the representation of small scale topographic features, such as localized tidal tributaries [7,22]. Using remote sensing data to generate tidal flat topography only provides information about the selected regions and times when the satellite images were acquired [23]. Thus, the use of an insufficient number of images in the waterline extraction method may lead to inaccuracies [2].

In this paper, we propose the pixel-based statistical analysis of intertidal classification and topography using multitemporal Sentinel-1 backscattering coefficients. Specifically, this study (1) distinguishes the tidal flats from land and ocean via the Jenks Natural Breaks (JNB) classification using Sentinle-1 SAR data; (2) estimates intertidal topography from a multitemporal SAR backscattering logistic model on a pixel-by-pixel basis; and (3) evaluates our Sentinel-1 intertidal topography against airborne Lidar and waterline-based Landsat TM/ETM+ Digital Elevation Models (DEMs). This study attempts to determine the feasibility of our innovative statistical analysis for intertidal classification and topography to monitor the near real-time spatiotemporal distribution changes of tidal flats through continuous and stable SAR data collection on local and regional scales.

2. Materials and Methods

2.1. Study Area

Our study area, Hwangdo (35°35′05.48″N, 126°22′39.92″E), is located in the innermost part of Cheonsu Bay on the central western coast of Korea and has a width of 1.65 km and length of 5.15 km (Figure 1). The Cheonsu Bay is characterized as a semi-closed bay with a shallow depth of less than 25 m. The sediment facies of this study area comprise mud flats, mixed flats, and sand flats, from the high tide waterline to the low tide waterline [9,24]. The total area of Cheonsu Bay has decreased from 380 to 180 km², which was caused by the embankment construction for land reclamation in the late 1980s and early 1990s [25,26]. The water depth around the Cheonsu Bay has been maintained at over 10 m since the dike was constructed. The extensive loss of tidal flats requires the continuous monitoring and assessing of the spatiotemporal change of the tide flat topography.

2.2. Tide

This study area predominantly has semi-diurnal tides. Tide data from Boryeong station were collected from the Korea Hydrographic and Oceanographic Agency (http://www.khoa.go.kr; accessed on 1 September 2021) for the years 2014–2020. Figure 2 shows the time series of the tide height and the acquisition times of the Sentinel-1 SAR images used in this study. The tide heights ranged from 6.68 m in mean higher high water (MHHW) to 1.22 m in mean lower low water (MLLW) for our study period. The highest astronomical tide (HAT) was 8.16 m and the lowest astronomical tide (LAT) was −0.19 m. The tide heights associated with the Sentinel-1 images were mainly located between 2 and 6 m in tidal range, which accounts for the lower and upper limits used when estimating our tidal flat topography.

2.3. SAR Image Processing

The Sentinel-1 mission is a two satellite constellation (i.e., Sentinel-1A and Sentinel-1B) with the same reference orbit [27]. Sentinel-1A and -1B were launched on 3 April 2014 and on 25 April 2016, respectively. In this study, 120 images were acquired from an ascending orbit on path 127 for the years 2014–2020. The data were Ground Range Multi-Look Detected (GRD) level-1 products acquired in Interferometric Wide Swath (IW) mode with High Resolution (HR) and VV polarization. The image processing was automated through multiple batch programs with a Graph Processing Tool (GPT). In general, the SAR data are provided as a digital number (DN), which are relative to the incidence angle and the highest intensity value. Radiometric correction was required to facilitate the comparative analysis with other studies. In this study, we converted the DN values to gamma naught (${\gamma }^0$) as dB values. Multi-look processing was not applied any further so as to preserve the intrinsic spatial resolution of the Sentinel-1 GRD HR products, which were generated by averaging five looks in range and one look in azimuth in the Singe Look Complex (SLC) product. A Lee Sigma filter with a 3 × 3 window size was adopted to reduce the speckle noise without losing spatial resolution. Terrain correction was applied to the Sentinel-1 images using the Shuttle Radar Topography Mission (SRTM) DEM data. The radar backscattering coefficient (${\gamma }^0$) images were converted into geographic coordinates with respect to the World Geodetic System (WGS84) reference system.

2.4. Statistical Analysis

2.4.1. JNB Optimization

The JNB method was adopted to separate the tidal flats from the ocean and land in our study area. The JNB method is a statistical analysis which can sectionalize data by minimizing each class’s average deviation from the class mean while maximizing each class’s deviation from the means of the other classes [28]. In this study, the algorithm began the process of dividing the ordered ${\gamma }^0$ values in terms of the corresponding tide heights. Secondly, the algorithm calculated the sum of squared deviations from the class means (SDCM) and the sum of squared deviations from the mean of the complete array data set (SDAM). The calculations had to be repeated using different breaks in the dataset to determine which set of breaks had the smallest in-class variance. The goodness of variance fit (GVF) index was indicative of the class’s fitness obtained from the JNB classification.

GVF = 1 − SDCM/SDAM

(1)

The GVF ranges from 0 to 1, where 0 indicates the worst and 1 is the best case. Therefore, tidal flats tend to be close to 1 whereas land and ocean are close to 0. This is supported by the fact that the radar backscattering coefficients in tidal flats are readily divided into two classes. Generally, the ${\gamma }^0$ shows high values when tidal flat topography is higher than tide height and low values when tidal flat is covered by water.

2.4.2. Logistic Function

After distinguishing the tidal flats from ocean and land in the SAR data with a certain threshold of GVF, the logistic model was applied to characterize the temporal evolution of the SAR backscattering coefficients and estimate the intertidal topography in this study. The temporal variation of the radar backscattering coefficients can be formulated by the logistic function [29] as below.

$${\gamma }^0\left(x, y\right)=\frac{k1}{1+e^{-a\left(h_i-h_t\right)}}+k2, i=1, \dots , M$$

(2)

where ${\gamma }^0\left(x,y\right)$ is the corresponding gamma naught at the tidal flat pixel (x, y), h_t is the tidal flat topography, h_i is the tide height and M is the number of SAR images used in the time series. The parameter a is related to the steepness of the logistic function, which ranged from −7 to −10 in the unit of m⁻¹ for our study area. The two parameters k1 and k2 are related to the average radar backscattering coefficients for the tidal flat exposure period at low tide and inundation period at high tide, respectively. In this study, our JNB optimization algorithm estimated both parameters k1 and k2 on a pixel-by-pixel basis. The JNB algorithm divides the data into two groups: high ${\gamma }^0$ values for the tidal flat exposure period and low ${\gamma }^0$ values for the inundation period with respect to the tide heights. The parameter k1 was calculated as the difference between the average ${\gamma }^0$ values of the tidal flat exposure and inundation periods. The parameter k2 was the average ${\gamma }^0$ for the inundation period. After the tidal flat topography h_t was initialized from our JNB algorithm, h_t was optimized in an iterative manner for the best fitting logistic model to derive our final Senintel-1 DEM.

2.5. DEM Evaluation

The Sentinel-1 DEM from the JNB optimization and logistic model was evaluated by the 2015 airborne Lidar campaign from the Korea Hydrographic and Oceanographic Agency (http://www.khoa.go.kr; accessed on 1 September 2021). The Sentinel-1 DEM was also compared to the Landsat DEM image generated by the waterline extraction method [9]. For comparison analysis, the 1 m airborne Lidar and 30 m Landsat DEMs on the Universal Transverse Mercator (UTM) projected coordinates were resampled to the 10 m spatial resolution of our Sentinel-1 DEM using bilinear interpolation. They were then converted into the WGS84 geographic coordinate system. The DEM evaluation was performed using common statistical indices, mean absolute error (MAE), and root mean square error (RMSE). The MAE is the average absolute values of the differences between the reference and the corresponding observation. The Root Mean Square Error (RMSE) is a quadratic scoring rule which measures the average magnitude of the error and gives a relatively high weight to large errors.

$$\mathrm{MAE} =\frac{\left|\left(y_r-y_o\right)\right|}{n}$$

(3)

$$\mathrm{RMSE} =\sqrt{\frac{\sum {\left(y_r-y_o\right)}^2}{n}}$$

(4)

where n is the number of observations and y_r and y_o are the reference and observation values, respectively. The RMSE tends to be larger or equal to the MAE. Our evaluation might have been influenced by different DEM resolutions and resampling methods, but their sensitivity to the DEM error was beyond the scope of this study.

3. Results and Discussion

3.1. Statistical Analysis

The multitemporal radar backscattering coefficients were processed to estimate the intertidal topography through JNB optimization and logistic function. In Figure 3, the mean and standard deviation of the multitemporal radar backscattering coefficients are illustrated. The order of the mean radar backscattering coefficients from highest to lowest was land, tidal flat, and ocean. Ocean exhibited low ${\gamma }^0$ values due to the specular scattering of the water surface. Interestingly, the order of the standard deviation of ${\gamma }^0$ from highest to lowest was tidal flat, ocean, and land. Tidal flats were characterized by high standard deviation, which represents high temporal variation over different tide heights. In contrast, the ocean and land areas had lower standard deviations due to the smaller temporal variability of the corresponding backscattering coefficients. This is related to the fact that tidal flat topography is generally located between land and ocean bottom topography, and tidal flats are exposed as land at low tide and covered by water at high tide.

For the intertidal classification, JNB optimization was applied to ocean, land, and tidal flats, considering that tidal flats tend to show a rapid change of Senitnel-1 SAR backscattering coefficients along with a steady change in tide height. Figure 3c shows the GVF map from the JNB optimization of the backscattering temporal variability. The tidal flat showed higher GVF values compared to the lower GVF values for ocean and land. This supports the fact that tidal flats show high variability in the backscattering coefficient due to the periodic in and out flux of the tide water, thus they can be separated from land and ocean using our JNB optimization. For the intertidal topography, the logistic function was implemented to characterize the temporal evolution of the SAR backscattering coefficients on a pixel-by-pixel basis with the model parameters of local terrain height, steepness, and tide height. Figure 3d shows the result of our intertidal topography, estimated from the best fitting logistic model for areas with high GVF values.

Figure 4 illustrates the temporal variability of the radar backscattering coefficients of ocean, land, salt pond, and tidal flats for more detailed results of our regression models. Tidal flats were divided into high, intermediate, and low tidal flats, in terms of elevation, to illustrate each of the temporal variabilities in Figure 4. Regardless of the tide height, ocean was characterized by low backscatter coefficients of around −15 to −20 dB, and the mean and standard deviation values were −17 ± 1.8 dB. Land showed relatively high gamma naught values from −10 to −5 dB, and the mean and standard deviation values were −7.9 ± 0.9 dB. Interestingly, salt ponds were clearly divided into two groups of high and low backscatter coefficients, with high standard deviation of ±4.4 dB. However, unlike tidal flats, this was not related to the tide heights and could not be formulated by our logistic function. This implies that the backscatter characteristics of salt ponds is closer to ocean when seawater is trapped and closer to land when seawater is left to evaporate in the sun. Figure 4d–f shows that each pixel of tidal flats is characterized by its own distinct temporal variability of the radar backscattering coefficients with different model parameters. Tidal flats show the mean and standard deviation values of −9.5 ± 3.0 dB, −14.3 ± 3.9 dB, and −15.6 ± 4.0 dB for high, intermediate, and low elevation terrains, respectively. This supports the notion that intertidal classification can be improved in an adaptive manner when considering that tidal flats comprise a wide variety of tidal flat sedimentary facies, textures, and exposure times. The error term from the best fitting logistic function may be explained by environmental factors, such as drainage rate, exposure time, and soil porosity [30]. Additionally, our method reduces the problem of choosing an appropriate fixed threshold to distinguish tidal flats from ocean and land in the Sentinel-1 SAR images.

3.2. DEM Evaluation

The accuracy of the Sentinel-1 DEM was evaluated against the airborne Lidar DEM from the Korea Hydrographic and Oceanographic Agency. In Figure 5, the mean RMSE and MAE were averaged into 0.1 bins of incremental GVF. The results show RMSE values, ranging from 0.1 to 2.6 m, and MAE values from 0.1 to 1.7 m. A drastic decrease in the RMSEs and MAEs is observed in GVF values higher than 0.2. Therefore, areas with GVF values higher than 0.2 were defined as tidal flats in this study and entered into our logical function to estimate intertidal topography.

We also compared our Sentinel-1 DEM to the Landsat DEM from waterline extraction, as well as the airborne Lidar DEM (Figure 6). The Senitnel-1 DEM outperformed the Landsat DEM in terms of the evaluation indices of RMSE and MAE. The Sentinel-1 DEM better captured high topographic reliefs (>5 m) compared to the Landsat DEM. Our pixel-based method allowed more accurate measurement of small morphologic changes than the waterline-based method. This can also be explained by the difference of intrinsic spatial resolutions for the 10 m Sentinel-1 images and the 30 m Landsat images.

Figure 7 illustrates the topographic profiles of the airborne Lidar, Sentinel-1, and Landsat DEMs demonstrated by the solid white line in Figure 6a. Overall, the topography profiles show a decreasing trend from north to south and the height difference amongst the three DEMs on the y axis from 2 to 5 m. The Seninel-1 DEM better captured small scale topographic reliefs in heights ranging from 3.5 to 4.5 m on the northern latitudinal profile (>36.57 in degree) whereas the Landsat DEM seemed to be smoothly interpolated with low reliefs. Also, our Sentinel-1 DEM well-represented the steep slopes on the southern latitudinal profile which was close to the boundary between tidal flats and ocean. This can be supported by the fact that the waterline-based method tends to produce low vertical precision on steep slopes due to the edge detection error of the waterline.

4. Conclusions

Tidal flats are very dynamic natural systems that experience short-term and long-term morphological and hydrological changes. Their highly dynamic behavior requires frequent monitoring. Ground-based measurements and UAV-based methods provide very accurate measurements, but they lack the ability to provide repeated observations. Satellite-based methods have been useful for accomplishing repeated observations for a large spatial coverage. This study presents our innovative statistical approach for intertidal classification and topography on a pixel-by-pixel basis using multitemporal SAR backscattering coefficients.

We demonstrated the use of JNB optimization to separate tidal flats from ocean and land using multitemporal Sentinel-1 backscattering coefficients. This method uses an iterative approach to find a distinct statistical behavior of the tidal flats along with tide heights. Tidal flats show a high variability of the SAR backscattering coefficients due to the periodic in and out flux of tidal water. In this study, areas with GVF values higher than 0.2 were classified as tidal flats, which was supported by the DEM evaluation against the airborne Lidar observations. We also successfully generated Sentinel-1 DEMs with 10 m spatial resolution and a height accuracy of 0.20 m in MAE and 0.28 in RMSE from the logistic function. The Sentinel-1 DEM outperformed the waterline-based optical image from the Landsat DEM. Our logistic function relies on a pixel-based algorithm to model the temporal variation of the multitemporal SAR backscattering coefficients, where three unknown parameters (i.e., k1, k2, h_t) were initialized from the JNB optimization. The height accuracy of our logistic model can be improved when considering the effect of tidal flat exposure time, sedimentary facies, and textures into the temporal variability of the SAR backscattering coefficients.

This study can also be used to generate important data as a means of managing the marine environment and maritime activities more coherently at the marine spatial planning (MSP) [31]. These results are essential to continue making potential improvements in quickly monitoring extensive tidal flats and coastal changes through continuous and stable SAR data collection on a local and regional scale. Future development could focus on more efficient models to estimate and forecast tidal flat topography and exposure with environmental factors, which could also lead to a significant improvement in investigating the spatiotemporal evolution of tidal flat productivity and habitat.