Automatic Semi-Global Artificial Shoreline Subpixel Localization Algorithm for Landsat Imagery

Abstract

Shoreline mapping using satellite remote sensing images has the advantages of large-scale surveys and high efficiency. However, low spatial resolution, various geometric morphologies and complex offshore environments prevent accurate positioning of the shoreline. This article proposes a semi-global subpixel shoreline localization method that considers utilizing morphological control points to divide the initial artificial shoreline into segments of relatively simple morphology and analyzing the local intensity homogeneity to calculate the intensity integral error. Combined with the segmentation-merge-fitting method, the algorithm determines the subpixel location accurately. In experiments, we select five artificial shorelines with various geometric morphologies from Landsat 8 Operational Land Imager (OLI) data. The five subpixel artificial shoreline RMSE results lie in the range of 3.02 m to 4.77 m, with line matching results varying from 2.51 m to 3.72 m. Thus, it can be concluded that the proposed subpixel localization algorithm is effective and applicable to artificial shoreline in various geometric morphologies and is robust to complex offshore environments, to some extent.

1. Introduction

The coastline, the boundary of land and sea, is one of the 27 most important land surface features, and is vulnerable to natural processes such as coastal erosion/accretion, sea level changes and human activities [1]. Coastline mapping is, therefore, becoming a fundamental work for coastal erosion monitoring, coastal resource management, coastal environmental protection and coastal sustainable development [2,3,4,5,6]. In reality, the shoreline accurate position is difficult to be localized, as the position changes continually through time, because of cross-shore and alongshore sediment movement in the littoral zone and especially because of the dynamic nature of water levels at the coastal boundary (e.g., waves, tides, groundwater, storm surges, setups, runups, etc.) [7].

With the advantages of cost-effectiveness and large spatial and temporal scales, satellite remote sensing data have been used widely for coastline mapping [1,7,8,9]. When shoreline changes are sufficiently large (several tens of meters), satellite remote sensing can enable semi-automated comparison of large-scale areas by providing a common protocol for all sites [10], thus making comparisons consistent [11]. However, most observed shoreline changes are presently much smaller [12], so that the coarse spatial resolution of pixels prevent the accurate determination of shoreline positions when monitoring shoreline changes [13]. In this case, shoreline change observations can only be obtained by means of repeated in situ surveys, analysis of aerial or satellite high-resolution photographs at several time intervals, or a combination of both approaches [11,14]. However, the expensive price and shortage of historical data cannot meet large-scale shoreline supervision demands and the requirements of shoreline change analysis. Thus, it is important to conduct research on how to accurately determine the shoreline’s position from long-term sequences of medium spatial resolution satellite images.

In recent years, many articles have appeared on how to use super resolution mapping (SRM) or subpixel edge localization (SEL) algorithms to extract the shoreline accurately. In these articles, positioning accuracy is quantitatively evaluated by four indicators: mean absolute error (MAE), standard deviation (SD), root mean squared error (RMSE), and line matching (LM) [15].

SRM has been applied to low- or medium-resolution satellite remote sensing images to overcome the limitation of the image spatial resolution of the original image. Li et al. [16] proposed that SRM can be categorized into two groups. The first group [17,18,19,20] is directly applied to satellite images instead of the intermediate spectral unmixing result, whereas the second group [21,22,23,24,25,26,27,28,29] is expected to improve the result’s accuracy when highly-accurate fraction images are available from a spectral unmixing model [30]. According to shoreline SRM, Foody et al. [31] presented a soft fuzzy classification utilizing a geostatistical approach to obtain accurate waterline locations. Muslim et al. [32] proposed a localized soft classification approach to predict the shoreline location by a two-point histogram and pixel-swapping algorithms. Muslim et al. [33] proposed a contouring and geostatistical method to geographically position the coastline within image pixels. Zhang et al. [34] integrated a geostatistical approach and the high-resolution spatial structure prior model to undertake super-resolution mapping, which can properly illustrate the spatial distribution of the coastline at a fine scale. Comparisons [35] have been made using three soft classification methods and three subpixel mapping methods for coastal area classification.

SEL algorithms are often designed as follows: first, the initial position is obtained by edge detection; second, a local edge model is adopted to refine the initial edge position to the subpixel level. Subpixel detection techniques can be grouped into three categories [36]: moment-based; least-squares-error-based; and interpolation-based. Concerned with subpixel shoreline localization, Pardo-Pascual et al. [37] extracted subpixel shorelines utilizing local spatial structures from Landsat TM and ETM+, where the RMSE obtained ranged from 4.69 to 5.47 m. Almonacid-Caballer et al. [38] determined the annual mean shoreline subpixel position from Landsat images, and the extracted shorelines were biased from the seaward direction by approximately 4–5 m. Qingxiang Liu [39] presented a subpixel vector-based shoreline method to monitor shoreline changes at Narrabeen–Collaroy Beach, Australia, over 29 years. The experimental results show that after the correction of tidal effects, the RMSEs of annual mean shorelines are within 5.7 m. Pardo-Pascual et al. [40] evaluated the accuracy of shoreline positions obtained from the infrared (IR) bands of Landsat 7, Landsat 8, and Sentinel-2 imagery on natural beaches, where the mean error reached 3.06 m (± 5.79 m) from Landsat 8 and Sentinel-2 images.

In fact, there are different types and various geometric morphologies of shorelines in complex offshore environments. From the point of view of different shoreline types, there are artificial shorelines and natural shorelines. A shoreline may also be considered over a slightly longer timescale, such as a tidal cycle, where the horizontal/vertical position of the shoreline could vary anywhere between centimeters and tens of meters (or more), depending on the beach slope, tidal range, and prevailing wave/weather conditions [7]. It is, therefore, more difficult and challenging to evaluate the shoreline location accurately, especially natural shorelines.

From a geometric morphology point of view, shorelines include simple straight, quasi-straight, and high curvature shorelines, or combinations of these. It is difficult for traditional subpixel shoreline algorithms to determine various geometric morphological shoreline subpixel positions. Most algorithms mentioned above obtain the most accurate shoreline position for simple straight shorelines, but fail for high curvature shorelines in which the positional error increases [37].

Owing to the complex offshore environment, which includes suspended sediment, foam, different land-cover types etc., there are many mixed pixels and noise along the shoreline. Thus, it is difficult to find pure pixels along the natural or artificial shoreline for spectral unmixing, and the pure pixels obtained by global spectral analysis may not represent shoreline local spectra finely. Meanwhile, these conditions may also lead to difficulties in modeling the local edge for SEL algorithms.

Compared with a natural shoreline, an artificial shoreline is stable and its position is not affected by tidal effects or other factors, thus the reference shoreline can be extracted from high-resolution satellite images from different imaging times. Thus, artificial shorelines have the advantage of being more easily validated than natural shorelines. In this study, we focus on how to determine the artificial shoreline position accurately. In addition, we propose a method called the semi-global shoreline subpixel localization (SGSSL) algorithm. The main thoughts underlying SGSSL are simplifying the shoreline subpixel localization problem to a segmented shoreline subpixel fitting problem, expressing a shoreline segment geometric morphology perfectly, and minimizing the intensity integral error in local windows. To express various geometric morphological shorelines, we utilize multi-scale corner points to divide the initial shoreline into relatively simpler shoreline segments. To prevent offshore environment interference on subpixel localization, we analyze the water index intensity homogeneity for designing local windows. In designed local windows, intensity integral errors are minimized to obtain the subpixel shoreline positions. The entire method is dependent not only on shoreline geometric morphology and global spectral features but also local window intensity analysis and segmented shoreline geometric morphology. Thus, the proposed method is named semi-global subpixel shoreline localization (SGSSL).

2. Study Areas & Datasets

2.1. Study Areas

With urbanized development, there are increasingly more artificial shorelines located along Chinese coastal areas. Caofeidian Port and the Xiamen coastal area were selected as the study areas. As shown in Figure 1a, the Caofeidian Port located at 118.5°E, 39°N, is adjacent to China’s Beijing Tianjin Hebei urban agglomeration and is one of China’s important ore transportation ports. As shown in Figure 1b, the Xiamen coastal area, located between 118°E–118.5°E and 24.35°N–24.6°N, is next to the Taiwan Strait, and Xiamen is an important port for international economic and cultural exchange. Five artificial shorelines of various geometric morphologies as experimental areas were chosen to evaluate SGSSL, and their key characteristic parameters are listed in

Table 1. Key characteristics of the five experimental areas.

		Experimental Areas 1/2	Experimental Area 3	Experimental Areas 4/5
Location		Caofeidian Port	Caofeidian Port	Xiamen coastal area
Shoreline type		artificial	artificial	artificial
Geometric morphology		simple straight	combination of quasi-straight and curved shape	high curvature/combination of quasi-straight and curved shape
Experimental image	Data	Landsat-8 OLI images (Path 122, Row 033)	Landsat-8 OLI images (Path 122, Row 033)	Landsat-8 OLI images (Path 119, Row 043)
	Date	04/25/2015	04/25/2015	10/13/2015
	Resolution	15 m fusion image	15 m fusion image	15 m fusion image
Reference image	Data	GF-2 image	GF-2 image	GF-2 image
	Date	05/31/2015	05/31/2015	02/06/2015
	Resolution	1 m fusion image	1 m fusion image	1 m fusion image

. The Gaofen-2 (GF-2) satellite data is selected as our reference data. GF-2 is the first civil optical remote sensing satellite independently developed by China with a spatial resolution better than 1 meter [41].

Our experiments verify the proposed algorithm from the following aspects. First, the subpixel shoreline localization results are superimposed on the original data to evaluate the visual effect of the proposed algorithm. Subsequently, compared with reference shoreline from GF-2, the four error indicators of the subpixel shoreline are calculated to verify the correctness and adaptability of the algorithm to different geometric morphology shorelines. Finally, the differences between the subpixel shoreline length and the reference shoreline length are calculated, which could illustrate ability of the proposed method to preserve details.

2.2. Data Pre-Processing

First, using the rational polynomial coefficient (RPC) of the GF-2 image, ortho-rectification of the GF-2 multi-spectral (MS) data and panchromatic (PAN) data were performed separately.

Then, the Gram Schmidt [42] pan sharpening algorithm was used to fuse the MS data with the PAN data; then, the spatial resolution of the Landsat 8 OLI fusion images was 15 m and that of the GF-2 fusion images was 1 m.

The registration parameters were estimated by correspondence feature points that were selected manually and by the polynomial model. The maximum registration error between the two fusion images (15 m/pixel fused Landsat8 OLI image, 1 m/pixel fused GF-2 image) is less than 3 m. The registration error will bring uncertainty to the accuracy assessment and is discussed in Section 5.1.

3. Materials and Methods

According to the main thoughts underlying SGSSL, the overall process is shown in Figure 2. First, global spectral and geometric morphology analysis are conducted: the initial shoreline is extracted using the Otsu [43] automatic threshold method from water index images; and geometric morphology control points, abbreviated as morphology control points (MCPs), are extracted using the multi-scale Harris algorithm [44]. The primary MCPs are utilized to divide the initial shoreline. Then, semi-global analysis is performed by the segmentation-merge-fitting (SMF) method. In the SMF process, the segmented shoreline subpixel location is determined by minimizing the intensity integral error, finally obtaining a continuous subpixel shoreline vector.

In Section 3.1, an ideal image is taken as an example to illustrate the basic principle of SGSSL. In Section 3.2, the challenges faced when the basic principle is applied to real satellite images are explained. In Section 3.3, the method of conducting the global analysis for SGSSL is introduced, including obtaining the initial shoreline position and extracting the MCPs. In Section 3.4, the process of performing a semi-global shoreline analysis for SGSSL is proposed, including designing local windows and details of the SMF processes.

3.1. Basic Principles of Subpixel Shoreline Localization

The entire shoreline can be regarded as consisting of many shoreline segments. The proposed subpixel shoreline localization algorithm is based on the following two assumptions:

Assumption 1: Any shoreline segment can be approximated by the polynomial function y = f(x).

Assumption 2: The shoreline segment divides the image into two homogeneous regions with intensities A and B (A < B).

The ideal binary image is built in the image coordinate system O-xy, as shown in Figure 3, in which there is one shoreline segment.

A local window (the purple box in Figure 3) centered on the ith shoreline pixel (x_i,y_i) is set, so the sum of intensity in the ith window is:

$$SU{M}_i={\displaystyle \sum _{j=y_i-m_2}^{y_i+m_1}{G}_{x_i,j}}{\begin{array}{c}\end{array}}_{(i=1,2,\dots n)},$$

(1)

where (x_i,y_i) is the current shoreline point’s pixel coordinate, G the pixel’s intensity, n the number of pixels in the shoreline segment, and m₁, m₂ are the pixels’ number above or under the shoreline csegment in the local window, respectively. According to Assumption 2, the integral of the intensity in the ith window is:

$$SU{M}_i{}^{*}=A\times S_{A_i}{}^{*}+B\times S_{B_i}{{}^{*}}_{(i=1,2,\dots n)},$$

(2)

where ${{\mathrm{S}}_{{\mathrm{A}}_{\mathrm{i}}}}^{*}$ and ${{\mathrm{S}}_{{\mathrm{B}}_{\mathrm{i}}}}^{*}$ are areas covered by A or B in the local window, respectively:

$$S_{A_i}{}^{*}+S_{B_i}{}^{*}{=m}_1+m_2+1.$$

(3)

According to Assumption 1, as the cubic function can express more geometric details and has superior morphological adaptability, the shoreline segment’s polynomial function is:

$$f\left(x\right)=a+bx+c{x}^2+d{x}^3.$$

(4)

So, the area under the shoreline segment in the ith local window ($S_{A_i}{}^{*}$) can be calculated as:

$$\begin{array}{l}{S}_{A_i}{}^{*}={\displaystyle {\int }_{x_i-1/2}^{x_i+1/2}\left(a+bx+c{x}^2+d{x}^3+0.5+m_2-y_i\right)}dx\\ =0.5+m_2-y_i+a+x_{i}b+\left(x_i{}^2+\frac{1}{12}\right)c+\left(x_i{}^3+\frac{1}{4}{x}_i\right)d\begin{array}{cc}& {}_{\left(i=1,2,\dots n\right)}\end{array}\end{array},$$

(5)

With the above derivations, the intensity integral of the ith local window’s (SUM_i^*) can be described as:

$$\begin{array}{l}SU{M}_i{}^{*}=A\times S_{A_i}{}^{*}+B\times S_{B_i}{}^{*}\\ =\left(1/2+m_2-y_i\right)A+\left(1/2+m_1+y_i\right)B+(A-B)a+(A-B)x_{i}b\\ +(A-B)\left(x_i{}^2+1/12\right)c+(A-B)\left(x_i{}^3+x_i/4\right)d\end{array},$$

(6)

In ideal conditions,

$$\{\begin{array}{c}SU{M}_i=SU{M}_i{}^{*} \\ {\displaystyle \sum _{j=y_i-m_2}^{y_i+m_1}{G}_{x_i,j}}=\left(1/2+m_2-y_i\right)A+\left(1/2+m_1+y_i\right)B+(A-B)a+(A-B)x_{i}b+(A-B)\left(x_i{}^2+1/12\right)c+(A-B)\left(x_i{}^3+x_i/4\right)d\end{array},$$

(7)

A similar idea has been researched by Trujillo-Pino et al. [36] for medical and indoor images, in which a subpixel edge location algorithm based on the partial area effect (PAE) was proposed. However, in that work, the algorithm neither considered the various geometric morphologies of the real shoreline/contour nor proved its application to actual satellite remote sensing images.

To solve Equation (7), the equation can be represented simply as:

$$p_i\beta =q_i{}_{\begin{array}{cc}& \left(i=1,2,\dots n\right)\end{array}}$$

(8)

where:

$$\begin{array}{l}{p}_i=(A-B)[\begin{array}{cccc}1& x_i& x_i{}^2+1/12& x_i{}^3+x_i/4\end{array}]\\ \beta ={\left[\begin{array}{cccc}a& b& c& d\end{array}\right]}^T\\ R_i=\left(1/2+m_2-y_i\right)A+\left(1/2+m_1+y_i\right)B\\ q_i=SU{M}_i-R_i.\end{array}$$

The following non-homogeneous equation can be obtained:

$$\begin{array}{l}\left[\begin{array}{cccc}{p}_{11}& p_{12}& p_{13}& p_{14}\\ p_{21}& p_{22}& p_{23}& p_{24}\\ ⋮& ⋮& ⋮& ⋮\\ p_{n1}& p_{n2}& p_{n3}& p_{n4}\end{array}\right]\beta ={\left[\begin{array}{c}{q}_1\\ q_2\\ ⋮\\ q_n\end{array}\right]}_{(i=1,2,\dots n)}\\ \beta =(P^{T}P{)}^{-1}{P}^{T}Q\end{array}$$

(9)

The intensity integral error in the ith local window (e_i) is defined as:

$$e_i=\left|SU{M}_i-SU{M}_i^{*}\right|={\left|q_i-p_i\beta \right|}_{(i=1,2,\dots n)}$$

(10)

By the least squares methodology, the cubic polynomial coefficients vector $\beta$ = (a,b,c,d)^T can be solved with intensity integral error minimization, and then the shoreline subpixel localization is determined.

3.2. Challenges for Subpixel Shoreline Localization in Remote Sensing Images

In Section 3.1, the derivation was based on the ideal image in Figure 3. When dealing with real shorelines, as Figure 4a shows, the initial shoreline (colored in purple) with various geometric morphologies does not satisfy Assumption 1. Therefore, the initial shoreline should be divided into segments of relatively simple morphologies, which can be approximated by the cubic polynomial function. For a shoreline to be divided appropriately, the multi-scale Harris corner algorithm [44] should be utilized to extract the MCPs.

When processing actual satellite images, the water index result is affected by sensor imaging noise and the interaction between adjacent classes. Meanwhile, heterogeneous pixels exist in the local window with intensity changes (Figure 4b, green box). Then, Assumption 2 is not satisfied, which leads to the intensity integral error expressed in Equation (10) not equaling zero.

It is therefore necessary to design the local window to guarantee the intensity integral error approaches zero to ensure the performance of SGSSL correctly.

3.3. Shoreline Global Analysis

3.3.1. Determination of Initial Shoreline Position

To make full use of satellite image global information, we will conduct a shoreline global analysis, which includes a spectral feature analysis and geometric morphology analysis to determine the initial shoreline position at the pixel level and extract appropriate MCPs.

3.3.2. Determination Initial Shoreline Position

Before extracting MCPs, we should determine the initial pixel level shoreline from the original fused satellite images. This procedure includes the following steps. First, the water index image is calculated. The modified normalized difference water index (MNDWI) [45] is preferred, and the reason why MNDWI is preferred is discussed in Section 5.1. Subsequently, the Otsu method [43] is applied to the water index image, in which an optimal threshold T* is selected automatically by maximizing the inter-class discrepancy. Third, using the optimal threshold value T*, the water index image is divided into a binary image, which includes non-water and water classes. A series of points representing the pixel level shoreline is obtained.

3.3.3. MCP Extraction

To satisfy Assumption 1, the initial pixel level shoreline with various geometric morphologies should be divided into segments with relatively simple morphology by MCPs. As multi-scale Harris detection [44] is sensitive to corners, MCPs are extracted by multi-scale Harris detection [44] from a binarized water index image:

$$M=\mathsf{\mu }\left(x,{\sigma }_I,{\sigma }_D\right)={\sigma }_D^{2}g\left({\sigma }_I\right)\otimes \left[\begin{array}{cc}{L}_x^2\left(x,{\sigma }_D\right)& L_x{L}_y\left(x,{\sigma }_D\right)\\ L_x{L}_y\left(x,{\sigma }_D\right)& L_y^2\left(x,{\sigma }_D\right)\end{array}\right],$$

(11)

where σ_I is the integral scale, g(σ_I) the Gaussian convolution kernel with integral scale σ_I, and L_a the derivative computed in the a direction. The multi-scale Harris cornerness measure combines the trace and the determinant of the scale-adapted second moment matrix:

$$cornerness=\det \left(\mu \left(X,{\sigma }_I,{\sigma }_D\right)\right)-\alpha trac{e}^2\left(\mu \left(X,{\sigma }_I,{\sigma }_D\right)\right)$$

(12)

where $\alpha$ is an empirical coefficient, $\alpha \in \left(0.04,0.06\right)$. It should be noted that decreasing the value of $\alpha$ will increase the cornerness. Since our aim is to extract MCPs from the water/non-water binary image, the value of $\alpha$ can be set to 0.06 (the maximum empirical value). The local maximum of cornerness at each scale determines the scales’ corner positions.

Second, each corner is verified depending on whether the Laplacian of Guassian (LOG) attains the maximum at the scale, and the LoG values are calculated by:

$$|LoG\left(X,{\sigma }_n\right)|={\sigma }_n^2|L_x{}_x\left(X,{\sigma }_n\right)+L_y{}_y\left(X,{\sigma }_n\right)|,$$

(13)

Comparing the LoG values with the adjacent two scale space images at the same position, if:

$$F\left(x,{\sigma }_n\right)>F\left(x,{\sigma }_l\right),l\in \left\{n-1,n+1\right\},$$

(14)

these corners would be reserved as multi-scale Harris corners.

Considering shoreline geometric morphological changes can be classified as dramatic variations and minor variations, the MCPs should include primary MCPs and supplementary MCPs. As the primary MCPs locate the positions at which the shoreline’s morphology changes drastically, theoretically, the multi-scale Harris corner points can be directly viewed as primary MCPs. However, during the process of building multi-scale image pyramids, images would be blurred and image structure details could be missed. Therefore, the initial scale’s corner subset is also preserved from multi-scale Harris detection [44] as supplementary MCPs. The partial area effect (PAE) subpixel algorithm [36] should be conducted for MCP positions to obtain their accurate subpixel position. These two types of MCPs and subpixel positions will be used in the SMF process.

3.4. Shoreline Semi-Global Analysis

To utilize the semi-global information of a shoreline segment correctly, we conduct shoreline semi-global analysis, which includes a designed local window and the SMF method to determine the local homogeneous intensity and obtain proper shoreline segments for subpixel localization.

3.4.1. Designing Local Window

In this subsection, we introduce how to design the local window [36] and estimate homogeneous intensities A and B [36] to ensure that the intensity integral error e_i in Equation (10) is close to zero, so as to satisfy Assumption 2.

First, the maximum gradient direction of each shoreline point is calculated, and for every shoreline point, Sobel edge detection is utilized to calculate the gradient (G_x, G_y). Then, the larger gradient is preserved as the points’ maximum gradient direction.

It should be noted that the shoreline segment consists of numerous points. Therefore, if shoreline points of a certain maximum gradient direction (G_x or G_y) have a larger proportion in the segment, then that direction will be used as the main direction for the segment.

If the main direction for the segment is G_y, the (m₁+m₂+1) ×1 local window is designed, and the cubic polynomial function of a segment is:

$$y=a+bx+c{x}^2+d{x}^3.$$

(15)

If the main direction for the segment is G_x, the 1×(m₁+m₂+1) window is designed and the cubic polynomial function of segment is:

$$x=a+by+c{y}^2+d{y}^3.$$

(16)

Second, since the homogeneous pixels’ intensities are stable, they have a minimum gradient in the local window’s direction. The algorithm adjusts m₁, m₂ to find the minimum gradient pixels (pixels in the blue box in Figure 4b). Once we find the minimum gradient pixels, their water index intensities are used to estimate the intensity A, B and their coordinates are used to set the window size [36].

In the ith window, the intensities A_i and B_i are the farthest pixels from the shoreline. To ensure the correlation between pixels in the local window, we limit $m_1\le 4,m_2\le 4$.

$$A_i=G_{i,j-m_2},B_i=G_{i,j+m_1},$$

(17)

Since there is a correlation between the intensities of adjacent points in the shoreline segment, to ensure that the homogeneity of A or B further prevents isolated noises, according to the adjacent shoreline points’ relative positions, we determine slope k of this shoreline point and calculate the more homogeneous intensity estimation values A_i^* and B_i^*:

$$\{or\begin{array}{c}{A}_i{}^{*}=\frac{A_i+A_{i+1}}{2},\hspace{1em}{B}_i{}^{*}=\frac{B_i+B_{i-1}}{2},\begin{array}{cc}& if\hspace{1em}k\ge 0\end{array}\\ A_i{}^{*}=\frac{A_i+A_{i-1}}{2},\hspace{1em}{B}_i{}^{*}=\frac{B_i+B_{i+1}}{2},\begin{array}{cc}& if\hspace{1em}k<0\end{array}\end{array},$$

(18)

where $k=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$, (x_i,y_i) and (x_i+1,y_i+1) are the adjacent shoreline points’ coordinates.

3.4.2. Segmentation-Merge-Fitting Method

The MCPs determined in Section 3.3.2 include the primary MCPs and the supplementary MCPs, which have been derived from multi-scale Harris corners to represent the major and minor shoreline geometric morphologies. However, utilizing them all without any selection it will lead to over-segmented shorelines and a high computational burden for the subpixel localization algorithm. The segmentation-merge-fitting (SMF) method is proposed to obtain appropriate shoreline segments adaptively, according to the residuals of the least-squares process and the main directions of the adjacent shoreline segments. In the SMF, the least-squares fitting process and the subpixel positions of the MCPs are added as constraints (Equation (19)) to connect adjacent shoreline segments and obtain the continuous subpixel shoreline.

$$\{\begin{array}{c}{y}_{mc{p}_i}-f(x_{mc{p}_i})=0,\\ y_{mc{p}_j}-f(x_{mc{p}_j})=0,\end{array}$$

(19)

The detailed steps are shown in Figure 5.

The detailed explanation for all of the above steps follows:

0. Preparatory work: Build the shoreline morphological control point set (SMCPS) and add all primary MCPs to the SMCPS; set the threshold t of the least squares residual to 0.08.

1. Segmentation: Utilize the SMCPS to divide shoreline to obtain shoreline segments and calculate the main direction of each shoreline segment.

2. Merge: Judge whether the main direction of the current shoreline segment is the same as the main direction of the adjacent segment:

If TRUE, merge these two adjacent segments and go to the Step 3 ‘Delete Point’; else, go to Step 4.

3. Delete point: Remove the current morphological control point connecting the two adjacent segments from the SMCPS.

4. Fitting: Polynomial coefficients are calculated by the constrained least squares methodology, and the least squares residuals are computed for each point. If there are four shoreline points adjacent to the shoreline MCP, the residuals of which are larger than the threshold t, this shoreline MCP should be removed from the least square constraints. Then, the current shoreline segment must be recalculated. Else, go to Step 5.

5. Judge the ‘re-segmentation’ condition: If, in a shoreline segment, shoreline points with least-squares residuals larger than the threshold t continuously appear and these points’ number is larger than 4, the shoreline segment must be re-segmented. Go to Step 6; else go to Step 7.

6. Add Point:

① If the segment is a merged shoreline in Step 2, the MCPs removed in Step 2 should be restored in the SMCPS.

② If the segment is not a merged shoreline in Step 2, the supplementary MCPs located in the segment are added into the SMCPS.

③ If there are no supplementary MCPs in the segment, the point with the largest least squares residual should be selected and added into the SMCPS.

Update the SMCPS, and go to Step 1.

7. Traversing: If all shoreline segments have been traversed, keep the subpixel results and go to Step 8; else, select the next shoreline segment and go to Step 2.

8. End: If the SMCPS is constant during this traversal, end the loop and go to Step 9; else, go to the Step 2 and re-traverse.

9. Output: The subpixel shoreline results are output.

3.5. Verification Method

The reference shorelines are extracted manually from GF-2 fusion images, which satisfy the standards of the “Technical Regulations for Satellite Remote Sensing Survey on Island & Coastal Zones” [14]. Using the SGSSL algorithm, a subpixel shoreline can be obtained. The reference shoreline and the shoreline determined by SGSSL can be compared.

We choose four error indicators to assess the SGSSL performance: the MAE, RMSE, SD and LM. The MAE (Equation (20)) is obtained by averaging all distance errors, because all the errors are obtained by calculating the absolute value of the distance from the SGSSL shoreline to the GF-2 reference shoreline. The MAE and RMSE describe the SGSSL result bias towards the reference shoreline. The SD indicates the variability around the MAE (Equation (21)):

$$MAE=\frac{{\displaystyle {\sum }_{i=1}^N\left|d_i\right|}}{N}$$

(20)

$$SD=\sqrt{\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{(d_i-MAE)}^2}}$$

(21)

$$RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^N{d}_i^2}}{N}}$$

(22)

where $\left|d_i\right|$ is the distance from the subpixel shoreline point to the reference shoreline.

As Figure 6 shows, for the calculation of the LM [15], S_Δ is the sum of the area enclosed by the SGSSL shoreline and the reference shoreline, and L_real is the length of the reference shoreline. In Figure 6, the black dotted line represents the reference shoreline and the solid line represents the shoreline determined by SGSSL.

$$LM=\frac{S_{\mathsf{\Delta }}}{L_{real}}$$

(23)

where $S_{\mathsf{\Delta }}=S_1+S_2+S_3+S_4$.

4. Results

Here, we evaluate the results of the proposed SGSSL in terms of visual comparison, quantitative assessment, and shoreline detail preservation ability.

4.1. Visual Comparison

There are various geometric morphological shorelines in the selected coastal experimental areas. In the first two experimental areas (Figure 7a,e), the initial shoreline colored in purple has been extracted by Otsu [43] in Figure 7b,f. Owing to the fact that shoreline morphology is simple, MCPs labeled by yellow crosses can divide the initial shoreline into relatively simpler segments (Figure 7b,f). The proposed SGSSL algorithm determines the subpixel shoreline, which is represented by the red line in Figure 7c,g and which coincides with the real shoreline well. Although in the local zoomed image (Figure 7d,h) the initial pixel level shoreline points in yellow are located slightly landward, the final subpixel shoreline results (red line) still locate accurately.

In the latter three experimental areas (Figure 7i,m,q), the shoreline morphology is relatively complex, and the initial shoreline in purple has also been extracted by Otsu [43] in Figure 7 j,n,r. With the SMF method, we keep the selected MCPs labeled with yellow crosses (Figure 7j,n,r) and using them, the shoreline can be divided into relatively simpler segments to be perfectly expressed by red line in Figure 7k,o,s. Although in the local zoomed image (Figure 7l,p,t) the initial pixel level shoreline points in yellow are located slightly landward, the final subpixel shoreline results also coincide with the actual position well.

4.2. Quantitative Assessments

Table 2. The quantitative assessments in experimental areas.

Experimental Area	MAE (m)	SD (m)	RMSE (m)	LM (m)
1	2.94	1.93	3.51	2.87
2	3.34	2.16	3.97	3.30
3	3.67	3.06	4.77	3.72
4	2.72	2.61	3.77	2.77
5	2.48	1.72	3.02	2.51

summarizes the quantitative assessment results. In all experimental areas, the MAE at the subpixel level lies in the range of 2.48–3.34 m with an average of 3.03 m; the RMSE varies from 3.02–4.77 m with an average of 3.80 m; while the LM lies in the range of 2.51–3.72 m, with an average of 3.03 m. All quantitative assessments prove that the proposed SGSSL is reliable.

It should be noted that these assessment results may be affected by registration errors. How the registration errors influence the quantitative assessments is discussed in Section 5.1.

4.3. Shoreline Detail Preservation Ability

Some shoreline details observed in the high-resolution images would be blurred or missed in low-resolution images. Hence, it is necessary to verify the detail preservation ability of SGSSL for shoreline details by comparing the lengths between subpixel results and high-resolution images.

In

Table 3. Subpixel shoreline length and reference shoreline length.

Experimental Area	1	2	3	4	5
Subpixel Shoreline Length (m)	3206.22	3000.19	2739.38	6324.41	3572.84
Reference Shoreline Length (m)	3205.25	2994.95	2675.57	6236.77	3571.32
Length Difference Ratio	0.03%	0.17%	2.33%	1.39%	0.04%

, the length difference ratios between subpixel results and reference shorelines are calculated, and the maximum value is less than 2.5%, which indicates that the proposed SGSSL can effectively preserve shoreline details.

5. Discussion

5.1. Registration Error Influence on Quantitative Assessment

There are unavoidable registration errors when the accuracy assessment is conducted between the reference shorelines and SGSSL results. The registration errors will bring uncertainty to the quantitative assessment of SGSSL.

For objective analysis, three typical shoreline geometric morphologies are chosen for the registration error effect analysis.

First, the upper part of experimental area 1 is chosen. Because the upper part of experimental area 1 is a nearly vertical shoreline, the displacement of the horizontal direction will bring an apparent influence to the final MAE. Using 1 m as the displacement interval (0.067 pixels for Landsat8 OLI data) and with the maximum displacement limited to 5 m, the registration error influence in the horizontal direction across the range of [–5 m, +5 m] can be calculated and observed in Figure 8a,b.

Second, the middle part of experimental area 5 is chosen. Because the middle part of experimental area 5 is a nearly horizontal shoreline, displacement in the vertical direction will bring an apparent effect on the final MAE. Using 1 m as the displacement interval (0.067 pixels for Landsat8 OLI data) and with the maximum displacement limited to 5 m, the registration error effect at the vertical direction in the range of [–5 m, +5 m] can be calculated and observed in Figure 8c,d.

As Figure 8a,b show, displacements in the right direction will compensate for the registration error influence on the MAE. The MAE without registration error compensation is 5.73 m, which is 2.59 m larger than the entire experimental area 1 MAE result. With registration error compensation in the right direction, the MAE decreases. At a displacement of 3 m in the right direction, the MAE is 2.73 m; at the same direction of displacement of 4 m, the MAE reaches its most accurate value, 1.75 m. After that, MAE results cannot be compensated.

As Figure 8c,d show, for the middle part of experimental area 5, displacement will also have an effect on the final MAE. The MAE is 1.69 m without any registration error compensation, which is smaller by 0.59 m than the entire MAE of experimental area 5. With displacements in the up or down directions, the worse MAE will be obtained. Additionally, at a displacement of 3 m in the up direction, the MAE is 3.39 m; at the displacement 3 m in the down direction, the MAE is 2.9 m.

Third, experimental area 4 is chosen. Because the artificial shoreline in experimental area 4 is almost circular, the registration error will not influence the SGSSL result in certain directions. So the reference shoreline is displaced in four different directions, namely up, down, left and right, with 1 m of the displacement interval (0.067 pixels for Landsat8 OLI data), and the maximum displacement is limited to 5 m. In Figure 9, the MAE results of SGSSL with different displacements are shown.

From Figure 9, it can be observed that the MAE becomes more accurate in the 270° direction, namely in the down direction. The MAE first becomes most accurate at 2.60 m, which is the minimum MAE. The MAE then increases to 2.63 m in the 270° direction with a displacement of 2 m, and increases to 2.77 m in the same direction with a displacement of 3 m. For other directions, the MAE results are not better. The worst subpixel accuracy appears in the direction of 90°, the up direction, and the worst MAE is 3.74 m.

From the above analysis, it can be concluded that the registration error brings uncertainty to the SGSSL quantitative assessment results. With the same registration method, different geometric morphological shorelines are influenced by registration errors to different extents. The lowest accuracy of the MAE appears at the vertical shoreline, and the MAE reaches at 5.73 m, which is still better than 0.5 pixels (15 m of the fused Landsat8 OLI image). Hence, the registration error with a maximum value of less than 3 m is acceptable for SGSSL.

5.2. Water Index

To select the water index with optimal positioning accuracy, the positioning errors of the SGSSL algorithm under three different water indices—the normalized difference water index (NDWI) [46], MNDWI, and automated water extraction index (AWEI) [47]—are calculated and compared.

As Figure 10 shows, the accuracies of shoreline positioning under three different water indices have all reached the subpixel level, indicating that the proposed algorithm is applicable to all water indices. It is obvious that the MNDWI is best in the selected experimental areas. One of the reasons is that the short-wave infrared 1 (SWIR₁,1566.50 – 1651.22 nm) band is used in the calculation of the MNDWI, and the most accurate and robust sub-pixel shoreline positioning results are often obtained using the SWIR₁ band [40]. Therefore, this paper prefers to use the MNDWI to enhance the differences between land and water, but, considering complicated offshore environments and data sources, in other coastal areas utilizing other water indices is acceptable.

5.3. Intensity Integral Error Analysis

Owing to sensor imaging noise and the interaction between adjacent classes, the intensity integral error in the ith local window is probably not equal to zero.

As Figure 11 shows, the initial shoreline pixel coordinate is (x_i,y_i), and the red line is the real shoreline that crosses the pixel (x_i,y_i). The subpixel level coordinates of the point in the real shoreline are (x₀,y₀), (x₀∈[x_i−1/2, x_i+1/2], y₀∈[y_i−1/2, y_i+1/2]). Once the local window size is determined by finding the minimum gradient pixels along the window direction, the window sizes m_1, m₂ and the homogeneous intensity A_i, B_i are all obtained. At shoreline point (x₀,y₀), the intensity profile is drawn along the window direction, presuming the direction lies in the y axes in Figure 11.

If the shoreline segment can be expressed by a cubic polynomial. S₁,S₂,S₃ are areas enclosed by the intensity profile and the y axis (window direction),

$$\begin{array}{l}{S}_1(x_0)=A_i\times \left(f\left(x_0\right)-\left(y_i-0.5-m_2\right)\right)\\ S_3(x_0)+S_4(x_0)=B_i\times \left(\left(y_i+0.5+m_1\right)-f\left(x_0\right)\right)\end{array},$$

(24)

Therefore, the sum of the ith local window intensity is:

$$SU{M}_i={\displaystyle \sum _{j=y_i-m_2}^{y_i+m_1}{G}_{x_i,j}}={\displaystyle {\int }_{x_i-1/2}^{x_i+1/2}\left(S_1(x_0)+S_2(x_0)+S_3(x_0)\right)}d{x}_0,$$

(25)

where, in the local window, G is the pixel intensity. The approximation of the sum of the ith local window intensity is:

$$\begin{array}{l}SU{M}_i{}^{*}=A_i\times S_{A_i}{}^{*}+B_i\times S_{B_i}{}^{*}\\ =A_i\times {\displaystyle {\int }_{x_i-1/2}^{x_i+1/2}\left(f\left(x_0\right)-\left(y_i-0.5-m_2\right)\right)}d{x}_0\\ +B_i\times {\displaystyle {\int }_{x_i-1/2}^{x_i+1/2}\left(\left(y_i+0.5+m_1\right)-f\left(x_0\right)\right)}d{x}_0\\ ={\displaystyle {\int }_{x_i-1/2}^{x_i+1/2}(S_1(x_0)+S_3(x_0)+S_4(x_0))}d{x}_0\end{array},$$

(26)

Thus, in the local window the intensity integral error e_i can be described as:

$$e_i=|SU{M}_i-SU{M_i}^{*}|=|{\displaystyle {\int }_{x_i-1/2}^{x_i+1/2}(S_2(x_0)-S_4(x_0))}d{x}_0|$$

(27)

where the intensity integral error e_i is related to three factors: the window size (m₁+m₂+1), the homogeneous intensity difference (B-A), and the intensity slope at the shoreline point (x₀,y₀). Regarding the three factors, the intensity slope is determined by image information. The other two factors are determined by appropriate homogeneous intensity estimation A, B and the window’s size (m₁+m₂+1), which can ensure that the e_i approaches zero.

Furthermore, to verify the correctness of the local window design method and to verify whether or not the intensity integral error e_i approaches zero, the relative error δ is calculated within the local window in Landsat OLI8 MNDWI images.

$$\delta =\frac{SU{M_i}^{*}-SU{M}_i}{SU{M}_i}.$$

(28)

The shorelines extracted manually from GF-2 images are viewed as the reference shorelines. In Equation (28), SUM_i^* is the integral of intensity in the ith window and calculated according to Equation (6). The reference shoreline coordinates are used to calculate S^*_Ai and S^*_Bi. SUM_i is the sum of intensity in the ith window and calculated according to Equation (3).

As Figure 12a shows, 6743 local windows covering different offshore environments over different periods of time in the study areas were sampled. The relative error distribution is shown in Figure 12b, the mean of δ is 0.0174 and the variance is 0.0278. The probability of δ less than 5% is 87.57%, and that of δ being less than 10% is 99.85%. Thus, it can be concluded that in most local windows, the relative error δ can be seen as a small number, whose absolute value approaches zero.

5.4. Segmentation-Merge-Fitting Process

As mentioned in Section 3.4.2, the initial shoreline may not be perfectly expressed by the cubic polynomial function. In Figure 13a, it can be observed that the initial shoreline could be divided into relatively simpler segments by primary MCPs (colored blue). As shown in Figure 13b,c (a magnified version of Figure 13b), after using primary MCPs and constrained least squares solving, one shoreline segment containing many shoreline points of larger fitting residuals still exists (marked by the yellow box in Figure 13b), which must be segmented further. In Figure 13d, all supplementary MCPs in green would be used to re-segment this problematic segment. Then, as the Figure 13e shows, some segments would be merged. Finally, appropriate shoreline segments that perfectly agree with polynomial functions are obtained using the SMF method, and are named seg1–seg3 in Figure 13f, the final SMCPS selected from MCPs are labeled by yellow crosses. As shown in Figure 13g,h (a magnified version of Figure 13g), the final subpixel positioning results coincide well with the real shoreline position.

Table 4. Quantitative assessments of subpixel shoreline segments.

Error Indicator	Primary Shoreline	Final Shoreline Segments Result
		Total	Seg 1	Seg 2	Seg 3
MAE (m)	10.53	2.58	3.22	2.12	2.66
SD (m)	12.11	1.84	2.20	1.26	1.97

lists the subpixel assessment indicators (MAE and SD) of shoreline segments marked by yellow boxes in Figure 13 during the SMF process. For the initial longer shoreline segment divided by primary MCP, the MAE and SD of subpixel localization results are 10.53 m and 12.11 m, respectively, which indicate that the initial subpixel localization result is problematic. With the SMF process, the final selected MCPs distribute appropriately along the shoreline and three correct segments are preserved. In addition, the subpixel localization MAE results lie in the range of 2.12 m to 3.22 m and the SD results from 1.84 m to 2.20 m.

5.5. Robustness to Complex Offshore Environment and Salt-And-Pepper Noises

In Figure 14, due to the complex offshore environment, for example, the existence of suspended sediment (Figure 14a), the local window is difficult to obtain, which will directly affect the subpixel localization accuracy.

In Figure 14a,b, suspended sediment situations can be grouped into three different categories depending on the concentration extent and accumulated area.

When the suspended sediment concentration is low (the regions in the blue boxes in Figure 14), the influence will be suppressed or even eliminated in MNDWI images, regardless of the size of the accumulated area.

When the suspended sediment concentration is high and the accumulated area is small with a width less than four pixels, this region (green boxes in Figure 14) will be regarded as the intensity variation region in the local window. In this situation, with the local window design method, the minimum gradient pixels would be found in water, whose intensity is homogeneous.

When the suspended sediment concentration is high and the accumulated area is large (in yellow boxes in Figure 14), with the local window design method, the homogeneous pixels will be selected directly in the suspended sediment region.

In conclusion, regarding the above three forms of suspended sediments, the homogeneous intensity estimations are dealt with effectively and will not reduce the subpixel localization accuracy.

Table 5. MAE of local regions in suspended sediment environments.

	low Concentration	High Concentration
		Large Area	Small Area
MAE (m)	0.96	2.28	3.55

summarizes the MAE results of the selected local suspended sediment region, which lie in the range of 0.96–3.55 m, proving that our proposed SGSSL is robust to suspended sediments to some degree.

Furthermore, we evaluate the proposed SGSSL performance under the influence of salt-and-pepper noises. In Figure 15, the yellow points are determined by the PAE subpixel algorithm proposed by Trujillo-Pino [36], and the red lines are determined from SGSSL, where the white and black points are salt-and-pepper noise. With increasing noise percentage, the PAE subpixel results become increasingly problematic. However, the results determined by SGSSL are always accurate and are not affected by increasing salt-and-pepper noise.

As shown in Figure 16, the positioning accuracy results of the two methods under different percentages of salt-and-pepper noise are quantitatively evaluated. With increasing noise percentage, it is obvious that the SGSSL algorithm exhibits better accuracy, proving that the proposed SGSSL is robust to salt-and-pepper noise to some extent.

6. Conclusions

With the merits of efficient, large-scale investigational capability, satellite remote sensing shoreline mapping plays an important role in the monitoring of coastal resource management. However, low spatial resolution, various shoreline geometric morphologies, and complex offshore environments prevent the accurate positioning of shorelines. In this study, therefore, we proposed a semi-global subpixel shoreline localization (SGSSL) algorithm for accurately determining artificial shorelines.

The proposed method utilizes not only global spectral information and shoreline morphological features, but also local water index homogeneity features and simplifies the entire shoreline subpixel positioning problem with a segmented shoreline fitting solution. The method considers the following factors: (1) MCPs are utilized to divide the initial shoreline into segments of relatively simple geometric morphology; (2) minimum gradient pixels are found to design a local window; (3) the intensity integral error is minimized in every local window within a segment to initially determine the subpixel location; and (4) the SMF process is presented to obtain the shoreline segments that can be perfectly expressed by a cubic polynomial function and to determine the final subpixel results.

In experiments, five artificial shorelines of various geometric morphologies from Landsat 8 OLI images were selected. The accuracy of the proposed method was evaluated using four error indicators: the MAE, SD, RMSE, and LM. The subpixel RMSE results are all less than 5 m, ranging from 3.02–4.77 m; and the LM results are all less than 4 m, ranging from 2.51–3.72 m, proving that subpixel shoreline accuracy obtained by the proposed method is stable over different experimental areas with various morphologies.

It can be concluded that the proposed algorithm is applicable to the various geometric morphologies of artificial shorelines and is robust to complex offshore environments and salt-and-pepper noise, to some extent.

Limitations of the proposed algorithm include the fact that its performance heavily depends on MCP distribution. Although the SMF process helps in obtaining optimum segments, in some experimental images a lack of MCPs will lead to irreparable subpixel accuracy loss. Another issue worth mentioning is that the proposed algorithm relies on the initial pixel level shoreline, which is the local window determination basis. More adaptation thresholding methodology should be applied to guarantee the initial pixel level shoreline’s correct position. Finally, the proposed SGSSL has only been verified on a selected artificial shoreline, other types of shoreline, for example, sandy shorelines and mangrove shorelines, have not been evaluated.

In future research, the performance of the method should be improved by a more flexible MCP extraction algorithm and a more reliable initial shoreline determination method. In terms of application prospects, the method will be combined with multi-source satellite images or ground truth data in the continuous monitoring of shoreline dynamics, coastal terrain mapping and other related research topics.