Urban Riverway Extraction from High-Resolution SAR Image Based on Blocking Segmentation and Discontinuity Connection

: An urban riverway extraction method is proposed for high-resolution synthetic aperture radar (SAR) images. First, the original image is partitioned into overlapping sub-image blocks, in which the sub-image blocks that do not cover riverways are regarded as background. Sub-image blocks covering riverways are then ﬁltered using the iterative adaptive speckle reduction anisotropic di ﬀ usion (SRAD) that introduces the relative signal-to-noise ratio (RSNR). The ﬁltered images are segmented quickly by the Sauvola algorithm, and the false riverway fragments are removed by the area and aspect ratio of the connected component in the segmentation results. Using the minimum convex hull of each riverway segment as the connection object, the seeds are automatically determined by the di ﬀ erence between adjacent pyramid layers, and the sub-image block riverway extraction result is used as the bottom layer. The discontinuity connection between river segments is achieved by multi-layer region growth. Finally, the processed sub-image blocks are stitched to get the riverway extraction results for the entire image. To verify the applicability and usefulness of the proposed approach, high-resolution SAR imagery obtained by the Gaofen-3 (GF-3) satellite was used in the assessment. The qualitative and quantitative evaluations of the experimental results show that the proposed method can e ﬀ ectively and completely extract complex urban riverways from high-resolution SAR images.


Introduction
As an important carrier of urban water resources, the urban riverway is a type of ground object formed by natural and anthropogenic factors and affects various important functions, including urban ecological landscapes, flood control, and the regulation of heat island effects [1]. In recent years, the use of remote sensing satellite images to extract urban riverways quickly and accurately has significantly improved the rational utilization, macro-monitoring, real-time planning, and future protection of urban water resources.
In remote sensing, the water body's spectral feature serves as the primary indicator in the identification and extraction of the riverway [2]. In the visible waveband, the water body has less light absorption and higher reflectivity [3]. However, physical and chemical characteristics of different water body segments may cause measurement variations, thereby increasing the difficulty of water body extraction. As the wavelength increases, the reflectivity of the water body decreases; when the wavelength exceeds 0.75 µm, it almost becomes an absorber. The wavelength used for Synthetic Aperture Radar (SAR) imaging is usually at centimeter-level, which is much larger than the wavelength. As a result, the water body behaves like a region with lower overall brightness in the SAR image [4].
As an active microwave sensor, SAR has pronounced advantages over optical sensors, providing wider coverage in all weather conditions [5]. SAR images have become the main data source for water body extraction.
Few studies have been conducted on SAR image urban riverway extraction, usually part of larger research delineating water bodies [6]. SAR-based water body extraction methods based on SAR image mainly include texture analysis, data fusion, and threshold-based segmentation. The texture analysis method uses local binary patterns and a gray-level co-occurrence matrix to describe the texture features of water bodies [7,8]. This method is a classical approach for water body extraction, but many calculations and uneven texture features in water bodies reduce its accuracy and applicability. In the data fusion method, water body extraction is achieved mainly through the fusion of optical images, thematic maps, digital elevation models, and other information [9][10][11]. Due to the introduction of varying data types, this approach has large random errors and can sometimes be impractical and unfeasible due to difficulties in obtaining fusion data. For the threshold-based segmentation method, the threshold used for delineating water features is obtained through the water body's spectral features in the SAR image. Due to its simplicity and efficiency, this approach has been widely used in water body extraction studies. For example, in [12][13][14][15][16], histograms, Otsu, maximum entropy, and clustering methods have been used to select the threshold for water body extraction. However, due to the inherent speckle noise in SAR images, the results yielded many false extractions. Some postprocessing methods are used to eliminate false extractions. However, if these false extractions are connected to real water bodies, the postprocessing methods would be invalid, seriously reducing the extraction accuracy. The most simple and effective way to solve this problem is to filter the SAR image in advance. For this purpose, filters, such as Lee, Frost, Kuan, and speckle reduction anisotropic diffusion (SRAD) in [17][18][19][20], can be used prior to the water body extraction.
Compared with lakes, reservoirs, wetlands, and other water bodies, urban riverways in high-resolution SAR images have meandering geometries and are usually accompanied by natural obstacles (such as tree shadows) and human-made infrastructure (such as bridges and dams), appearing in the image as curved riverway segments with certain discontinuities [21]. The urban riverway spans over a large image region, while the ratio of the riverway object to the background is very small, resulting in the whole image being extracted with a larger scale and more complex background. These properties of urban riverways in high-resolution SAR images make the above methods unsuitable for extraction. To quickly and effectively extract riverways from large-scale urban SAR images, an extraction method based on blocking segmentation and discontinuity connection is proposed. The method first partitions the original image into a collection of overlapping sub-image blocks and extracts riverway segments in sub-image blocks by iterative adaptive SRAD and the Sauvola local threshold algorithm. The discontinuity between the riverway segments is then automatically connected, and the extraction results of each sub-image block are stitched to obtain the final extraction result.
The main methodological improvements and contributions of this paper are as follows.
(1) The original image is partitioned into sub-image blocks with a shared boundary region. These sub-image blocks are then classified into two categories according to whether they cover the riverway. As not all sub-image blocks cover the riverway, two schemes can be implemented on the classified sub-image blocks: the sub-image blocks that do not cover the riverway can be processed as background, while the sub-image blocks that cover the riverway are the processing units to be extracted. By introducing image blocking and a classification processing mechanism, the limited computing resources can be concentrated on specific regions of interest, which can improve the efficiency of riverway extraction. This can also effectively reduce background information not relevant in riverway extraction, thereby improving accuracy. Due to the shared boundary region of adjacent sub-image blocks, the riverway segment occupying these blocks has a buffer that can be better stitched during segment extraction.
(2) For the sub-image blocks covering the riverway, the iterative adaptive SRAD that introduces RSNR is used for filtering as a preprocessing step. The riverway segments are quickly extracted by the Sauvola algorithm [22], and the false riverway fragments are removed by the area and aspect ratio of connected components. The filtering procedure can greatly reduce the effect of speckle noise in SAR images during segmentation. The introduction of the RSNR quantitative index in SRAD to control the number of filtering iterations in different sub-image blocks can reduce artificially specified errors. Using the Sauvola algorithm in segmentation, the filtered image can overcome the uneven brightness of the image. Likewise, the use of morphological features of riverways in postprocessing the segmentation results can effectively improve the extraction accuracy.
(3) For each extracted sub-image block, the minimum convex hull in each riverway segment is calculated to form the search sub-image block with a convex hull as the connection object used as the bottom stratum of the multi-layer pyramid [23]. The seeds between the riverway segments in the pyramid's top layer are located by the search strategy. The seeds are then grown by the growth strategy, and the growth results are mapped in the next layer. The process of seed location and growth is repeated layer by layer until the bottom layer realizes the automatic connection of the discontinuities. By taking the minimum convex hull that envelops each riverway segment as the connection object, the meandering riverway can be taken as a whole to avoid a riverway being connected at the meanders. By constructing pyramids for multi-layer region growth, the discontinuity between riverway segments can be connected automatically.
The rest of this paper is organized as follows. Section 2 gives a detailed description of the high-resolution SAR image urban riverway extraction process. In Section 3, the experimental results are shown. Finally, a short conclusion is given in Section 4.
The proposed urban riverway extraction method consists of five parts: (1) image block partition and classification, (2) SRAD filtering preprocessing, (3) riverway segments extraction, (4) discontinuity connection, and (5) output of extraction result. The framework of this method is shown in Figure 1.

Overlapping Block Partition and Classification
The original image is first partitioned into a collection of non-overlapping, fully covering image blocks. The image blocks are then extended to adjacent blocks along the X-and Y-axes of the

Overlapping Block Partition and Classification
The original image is first partitioned into a collection of non-overlapping, fully covering image blocks. The image blocks are then extended to adjacent blocks along the Xand Y-axes of the image space with overlapping degrees (v 0 and h 0 ) to obtain a set of sub-image blocks. The overlapping block partition process can be expressed as where B = {B l , l = 1, . . . , h × v} and B' = {B' l , l = 1, . . . , h × v} represent the partitioned non-overlapping and overlapping image block collections, l is the index of the image block, and v and h are the number of image blocks along the Xand Y-axes.
where · is the down rounding operator, and V and H are the sizes of the non-overlapping block set along the Xand Y-axes.
The sub-image blocks in B' can be classified into blocks covering riverways and blocks not covering riverways by visual interpretation. Suppose that the number of sub-image blocks B' l 1 covering riverway and B' l 2 not covering riverway in B' is n 1 and n 2 , where l 1 ∈ {1, . . . , n 1 }, l 2 ∈ {1, . . . , n 2 }. For B' l 2 , which is processed as background, the result is denoted as J l 2 = {J i, j , (i, j) ∈ Ω l 2 }, where J i, j = 0 means background, and Ω l 2 is the spatial domain corresponding to B' l 2 .

SRAD Filtering Preprocessing
For B' l 1 , the spatial domain is denoted as Ω l 1 , l 1 ∈ {1, . . . , n 1 }. Taking I i, j , (i, j) ∈ Ω l 1 , as input, SRAD uses the following update function for the iterative filtering, where ξ and ξ are time and spatial steps, respectively. The parameter c t−1 i,j is the diffusion coefficient at (i, j) for the (t−1)-th iteration and is expressed as where and 2 represent the gradient and Laplacian operator, respectively, and c(q) is defined as where q i, j is the instantaneous diffusion coefficient, which can be calculated by the equation The smoothing scale function, q 0 (ξ), can be approximated by where ρ is the speckle reduction exponential, and 0 < q 0 ≤ 1 is the initial diffusion coefficient. In order to adaptively control the number of filtering iteration, RSNR is introduced as a termination condition, defined as where I t i,j and I t−1 i,j are the results of the t− and (t−1)-th filtering of I at (i, j), respectively. When the filtering iteration satisfies the Equation (9), the iteration is terminated, and the corresponding filtered image where ε is the preselected threshold, which is generally adopted as ε ≥ 0.01.

Extraction of Riverway Segments
Considering the spectral features of the riverway and the uneven brightness of the filtered image, the Sauvola algorithm is used to directly binarize the filtered image. The false riverway fragments are removed to realize the extraction of riverway segments.
In the optimized Sauvola algorithm, the integral image of F l 1 is used to quickly obtain the average pixel intensity µ i, j and the standard deviation σ i, j in a w × w window centered on (i, j). To calculate the threshold T i, j at (i, j), the following equation is used, where (i, j) ∈ Ω l 1 , κ ∈ [0.2, 0.5] is the empirical adjustable parameter. To facilitate the subsequent processing, the T i, j is used to binarize F l 1 : where R i, j = 0 means background, and R i, j = 1 means riverway, R l 1 = {R i, j , (i, j) ∈ Ω l 1 }. The connected component search is performed on the obtained R l 1 , and the false riverway fragments are removed by setting the area and aspect ratio thresholds of the connected component: where z is the index of the connected component and also serves as the riverway segment index, z ∈ {1, . . . , Z}; Z is the total number of connected components in R l 1 ; T a is the area threshold; and T r is the aspect ratio threshold. If Equation (12) is satisfied, the z-th connected component is considered as a true riverway segment. After traversing all connected components, the number of true riverway segments is denoted as Z', and the extraction result is R'

Discontinuity Connection between Riverway Segments
As other objects often overshadow the upper parts of urban riverways, a certain distance discontinuity is generated between the extracted riverway segments. In order to obtain complete riverway information, it is necessary to automatically connect the discontinuities between riverway segments.

Construction of the Convex Hull
To better correspond to the complex shapes of urban riverways, the minimum convex hull for each riverway segment is first calculated and regarded as the connecting object.
Remote Sens. 2020, 12, 4014 6 of 23 Take a riverway segment and its nine edges (E 0 , E 1 , . . . , E 8 ) as example. E 0 , the point with the smallest Y-coordinate value, is taken as the reference point. The other points E 1 -E 8 are connected to E 0 and sorted according to their angles with the X-coordinate, as shown in Figure 2a. The coordinates of three adjacent edge points form a third-order square matrix, where the determinant can be calculated. For example, given E 0 (x 0 , y 0 ), E 1 (x 1 , y 1 ), and E 2 (x 2 , y 2 ), the determinant of the square matrix formed by these points can be calculated as where | · | is the determinant operator. If P < 0, then E 1 is determined as a convex hull point. As shown in Figure 2, all the edge points are processed (the green lines indicate the intermediate process, and the red lines are determined convex hull borders) to obtain the convex hull of the riverway segment.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 27 where | · | is the determinant operator. If P < 0, then E1 is determined as a convex hull point. As shown in Figure 2, all the edge points are processed (the green lines indicate the intermediate process, and the red lines are determined convex hull borders) to obtain the convex hull of the riverway segment. Using this method to process the R'l 1 , the results are denoted as Pl 1 = {Pl 1 (z), z = 1, …, Z'}, where l1 ∈ {1, …, n1}, Pl 1 (z) is the convex hull of the z-th riverway segment in R'l 1 .

Pyramid Representation of Convex Hull Image
Taking Pl 1 as input, a Gaussian Pyramid (GP) is constructed by Gaussian smoothing and sub-sampling Pl 1 .
where Gk is the k-th layer of GP; k ∈ {1, …, K}, K is the total number of layers; Hr is the 2-D Gaussian low-pass filter for reducing resolution; ⊗ is the convolution operation; and S is the step size, such that (↓S)[ · ] means downsampling with step size S. The Connection Pyramid (CP) is constructed by upsampling the Gk and superimposing it with  Using this method to process the R' l 1 , the results are denoted as is the convex hull of the z-th riverway segment in R' l 1 .

Pyramid Representation of Convex Hull Image
Taking P l 1 as input, a Gaussian Pyramid (GP) is constructed by Gaussian smoothing and sub-sampling P l 1 .
where G k is the k-th layer of GP; k ∈ {1, . . . , K}, K is the total number of layers; H r is the 2-D Gaussian low-pass filter for reducing resolution; ⊗ is the convolution operation; and S is the step size, such that (↓S)[·] means downsampling with step size S. The Connection Pyramid (CP) is constructed by upsampling the G k and superimposing it with the G k-1 , k ∈ {1, . . . , K}. When k = K, where C k represents the k-th layer of CP; C k is the extension of C k to match the size of G k-1 , H e is the 2-D Gaussian high-pass filter for resolution expansion, and (↑S)[·] means upsampling with step size S.

Multi-Layer Region Growth
The resolution difference of adjacent GP layers is used to automatically determine the seeds of each CP layer. The discontinuity connection is realized by combining the region growth results of seeds in each CP layer.

(a) Automatic location of seeds
For any layer C k of CP, three kinds of pixel intensity values can be obtained based on its construction process.
, indicating that the position (i, j) may belong to the discontinuity between riverway segments. Therefore, the search strategy for the seeds in C k , k ∈ {1, . . . , K} is to select the points with C k (i, j) = 1 as the seeds.

(b) Growth range of seeds
With the increase of k, the growth range of seeds should be gradually reduced, and the maximum growth range of the k-th layer can be expressed as where S is the sampling step size for constructing GP.
(c) Growth strategy of seeds For any seed in C k , k ∈ {1, . . . , K}, the seed can grow in four directions along in eight angles with its position as the center, as shown in Figure 3. The four directions are horizontal (1↔5), right diagonal (2↔6), vertical (3↔7), and left diagonal (4↔8). The seed's growth strategy is to expand along each growth angle from seed until the pixel with intensity value 2 is found in its maximum growth range and there is no pixel with intensity 0 on the growth path. On each growth angle, the length of the growth path in pixels is denoted as If the growth results in any of the four directions meet the discontinuity condition, the seeds are connected to the riverway segment. Taking the horizontal direction as example, when the length of the growth path r 1 k + r 5 k ≤ R k , the seed can be determined to belong to the discontinuity between riverway segments. The connection process is employed to change the seed into a riverway segment, such that Starting from C K , the determined seeds are used for region growth, and the growth results are passed to the next layer to connect the discontinuities, layer by layer. When k = 1, C 0 is obtained, which is the connection result of P l 1 . For the final riverway extraction image , J i, j = 1 means riverway, and J i, j = 0 means background. If there is a little background noise inside the riverway extracted in J l 1 , the connected component method of Section 2.3 can be used to optimize it.
where S is the sampling step size for constructing GP.
(c) Growth strategy of seeds For any seed in Ck, k ∈ {1, …, K}, the seed can grow in four directions along in eight angles with its position as the center, as shown in Figure 3. The four directions are horizontal (1↔5), right diagonal (2↔6), vertical (3↔7), and left diagonal (4↔8). The seed's growth strategy is to expand along each growth angle from seed until the pixel with intensity value 2 is found in its maximum growth range and there is no pixel with intensity 0 on the growth path. On each growth angle, the length of the growth path in pixels is denoted as r o k , where o ∈ {1, …, 8} is the growth angle index.

Riverway Extraction Result Output
In order to obtain the output image in the riverway extraction, the resulting sub-image blocks have to be stitched. Considering the accuracy and computational efficiency, the output operation can be simply regarded as the inverse process of image-overlap partitioning.
First, for the riverway extraction results of each sub-image block, the overlapping degrees of h 0 and v 0 pixels are removed along the Xand Y-axes of the image space. The non-overlapping riverway extraction results are then stitched to form the output image. The process of riverway extraction can be expressed as where J l and J' l are the l-th overlapping and non-overlapping image blocks, respectively, and l ∈ {1, . . . , h × v}, J = {J i, j , (i, j) ∈ Ω} is the final riverway extraction output image, where J i, j ∈ {0, 1}, J i, j = 1 means riverway and J i, j = 0 means background.

Experiment and Analysis
For this study, a personal computer with Intel (R) Core (TM) CPU 32G, Windows 10 32-bit professional operating system, and MATLAB R2018b was used in the experiments. The experimental results were then evaluated quantitatively and quantitatively.

Data
The GaoFen-3 (GF-3) satellite is China's first C-band multi-polarization SAR satellite with a 1-m resolution. It was launched on 10 August 2016 and can receive and transmit fully polarized electromagnetic waves with 12 imaging modes, such as Spotlight (SL), Fine Strip (FS), Ultra Fine Strip (UFS), and Standard Strip (SS).
As a case study, a partial region of a GF3-UFS intensity image was selected, which was taken in Beijing, China in September 2017, with a size of 2800 × 4000 pixels and a 3-m spatial resolution. Figure 4 shows the original image of the study area, which contains objects, such as riverways, roads, buildings, and bridges that cross the riverway and connect to the road. Compared with other objects, the riverway shows a lower gray value and has a meandering shape.

Experiment and Results
Taking the partition parameters of Figure 4

Experiment and Results
Taking the partition parameters of Figure 4  Three sub-image blocks-B' 1,3 , B' 1,4 , and B' 3,4 in B' l 1 -with different riverway morphology were selected to display the results of each processing step.
The parameters of SRAD for each sub-image block in B' l 1 were set as ξ = 0.5, ξ = 1, ρ = 0.1, and q 0 = 0.5. The filtering results of the selected three sub-image blocks are shown in Figure 6a1-c1. The speckle noise was effectively filtered out, and the water and land boundaries had been well preserved. Figure 6a2-c2 are the results of segmenting the filtered outputs using the Sauvola algorithm, where w = 50 and κ = 0.3. From the segmentation results, there were some false riverway fragments. The false riverway fragments can be removed by taking the thresholds T a = 400 and T r = 1.5 of the connected components to obtain the extraction results of the riverway segments, as shown in Figure 6a3-c3.
As shown by the extraction results of the riverway segments in Figure 6a3-c3, there were discontinuities between the riverway segments. Some riverway segments showed branching and annular distribution (see green ellipses in Figure 6b3), which can interfere with the connection of discontinuities. To accommodate the complex riverway morphology, each extracted riverway segment is tightly enveloped with a minimum convex hull, as shown in Figure 7a1-c1. As shown in the figures, the convex hull can effectively take the annular riverway segments as a whole. Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 27 The parameters of SRAD for each sub-image block in B'l 1 were set as ξ = 0.5, ξ = 1, ρ = 0.1, and q0 = 0.5. The filtering results of the selected three sub-image blocks are shown in Figure 6(a1)-(c1). The speckle noise was effectively filtered out, and the water and land boundaries had been well preserved. Figure 6(a2)-(c2) are the results of segmenting the filtered outputs using the Sauvola algorithm, where w = 50 and κ = 0.3. From the segmentation results, there were some false riverway fragments. The false riverway fragments can be removed by taking the thresholds Ta = 400 and Tr = 1.5 of the connected components to obtain the extraction results of the riverway segments, as shown in Figure 6(a3)-(c3).
As shown by the extraction results of the riverway segments in Figure 6(a3)-(c3), there were discontinuities between the riverway segments. Some riverway segments showed branching and annular distribution (see green ellipses in Figure 6(b3)), which can interfere with the connection of discontinuities. To accommodate the complex riverway morphology, each extracted riverway segment is tightly enveloped with a minimum convex hull, as shown in Figure 7(a1)-(c1). As shown in the figures, the convex hull can effectively take the annular riverway segments as a whole.  The connection for each riverway segment was achieved by combining each layer of the growth results (S = 3, K = 5), as shown in Figure 7(a2)-(c2). As shown in Figure 7(b2), the construction of a convex hull can avoid the internal connection of annular regions (see green ellipses). For discontinuities caused by dense scattering points near the riverway (see yellow ellipse in Figure 7(c2)), the number of seeds located was limited due to the linearity of the convex hull border, resulting in insufficient connection.
For the qualitative evaluation of the extraction results, the extracted riverway was represented as red and superimposed on the original image, as shown in Figure 7(a3)-(c3). The superposition results show that the riverway was completely extracted, and the boundary of the extraction result agrees well with real riverway boundaries.
For the quantitative evaluation of the object regions, the Sorensen-Dice and Jaccard similarity metrics are used [24,25]. The Sorensen-Dice similarity evaluates the extraction quality based on all pixels inside the riverway regions, which is computed as where Jg is the corresponding ground truth of J, | · | denotes the number of elements in the set, dice(J, Jg) ∈ [0, 1], the higher the value, the better the extraction result. The Jaccard similarity has a relation with the Sorensen-Dice similarity as follows, where jaccard(J, Jg) ∈ [0, 1]; the higher the value, the better the extraction result. The connection for each riverway segment was achieved by combining each layer of the growth results (S = 3, K = 5), as shown in Figure 7a2-c2.
As shown in Figure 7b2, the construction of a convex hull can avoid the internal connection of annular regions (see green ellipses). For discontinuities caused by dense scattering points near the riverway (see yellow ellipse in Figure 7c2), the number of seeds located was limited due to the linearity of the convex hull border, resulting in insufficient connection.
For the qualitative evaluation of the extraction results, the extracted riverway was represented as red and superimposed on the original image, as shown in Figure 7a3-c3. The superposition results show that the riverway was completely extracted, and the boundary of the extraction result agrees well with real riverway boundaries.
For the quantitative evaluation of the object regions, the Sorensen-Dice and Jaccard similarity metrics are used [24,25]. The Sorensen-Dice similarity evaluates the extraction quality based on all pixels inside the riverway regions, which is computed as where J g is the corresponding ground truth of J, | · | denotes the number of elements in the set, dice(J, J g ) ∈ [0, 1], the higher the value, the better the extraction result. The Jaccard similarity has a relation with the Sorensen-Dice similarity as follows, where jaccard(J, J g ) ∈ [0, 1]; the higher the value, the better the extraction result.
For the quantitative evaluation of the object boundaries, the assessment criteria proposed by Modava et al. [26] was used. Hand-drawn riverway boundary lines were used as standard data, and the evaluation region with a 4-pixel radius was established as the center. The cumulative percentage of the extracted riverway boundary falling into the evaluation region with different radius was calculated using the formula where λ is the index of the evaluation region radius, buf λ is the evaluation region with radius λ, bor is the riverway boundary obtained from the experiment when λ = 0, buf is the standard riverway boundary line, and A λ is the cumulative percentage of overlap within λ. The accuracy evaluation results are presented in Table 1. As shown in the table, the evaluation values of each sub-image block for Sorensen-Dice and Jaccard similarity metrics are more than 90%, indicating that the proposed method has good regional accuracy. For the boundary accuracy, about 43% of the boundary lines extracted by each sub-image block are accurately located on the manually drawn boundary lines. When the radius of the evaluation region reached 2 pixels, the locational accuracy of the riverway boundaries increased to more than 90%, indicating high boundary accuracy achieved by the proposed method.
The processing results of each sub-image block are stitched according to the defined stitching line to obtain the riverway extraction result for the entire image, which is represented as red and superimposed on the original image, as shown in Figure 8a. In addition, compared with the 1-D Otsu method of the work in [12], the 2-D Otsu method of the work in [13], and the FCM method of the work in [16], the results are shown in Figure 8c [13], and FCM method in [16].
From the superimposed images of the extraction results of each method, it can be seen that in the 1-D Otsu method of the work in [12], there are many mis-extractions connected to the real riverway. This is because Otsu is a global threshold method, and the optimal threshold determined  [12], 2-D Otsu method in [13], and FCM method in [16].
From the superimposed images of the extraction results of each method, it can be seen that in the 1-D Otsu method of the work in [12], there are many mis-extractions connected to the real riverway. This is because Otsu is a global threshold method, and the optimal threshold determined by Otsu in large-scale images with complex background is not suitable for some local regions. The 2-D Otsu method in [13] considers the texture information of the pixel neighborhood, which improves the extraction accuracy. The FCM method in [16] considers the neighborhood information of pixels and the class fuzziness of pixels, and its extraction is better than 1-D and 2-D Otsu methods. The proposed method can make full use of the local information of the image because of the sub-image blocks processing strategy, so that there is almost no mis-extraction in the result, and because of the discontinuity connection measures, the riverway segments are connected together. This is of great significance for subsequent processing such as hydrological analysis and riverway planning.
The quantitative evaluation of each method is shown in Table 2. It can be seen from the table that the proposed method gives the highest Sorensen dice and Jaccard scores, reaching 93.97% and 88.63%, respectively, indicating that it has good regional extraction accuracy. With the increase of the evaluation radius, the boundary accuracy advantage of the proposed method is more obvious.

Discussion
In order to verify the feasibility and efficiency of the proposed method, the performance of image block processing strategy and sub-image block filtering was analyzed and evaluated.

Block Processing Strategy
The SRAD filter preprocessing was used as an example to determine the usefulness of the image blocking processing strategy proposed in this paper. The necessity and advantage of boundary overlap are illustrated by the stitching results of adjacent sub-image blocks. The extraction processing of all sub-image blocks is used to observe the impact of the operation of classifying sub-image blocks as covered or uncovered riverways through visual interpretation on the final river channel extraction results.
The goal of SRAD is to adjust the intensity value of any pixel using Equation (3) in order to minimize image noise. For the computational complexity of SRAD, the c(q) at (i ± 1, j ± 1) needs to be calculated when smoothing points (i, j). From Equation (4), the process involves square operation, which means that the computational complexity of SRAD filtering (Figure 4) in each iteration would not be less than (2 × M × N) 2 . The image was processed using the proposed blocking strategy, and the study area was partitioned into 16 sub-image blocks (see Figure 5). The computational complexity of the SRAD filter was less than 16(2 × M/4 × N/4) 2 , as the sub-image blocks that do not cover riverways are only used as background and not for filtering. The results were then graphed and compared using the number of pixels as abscissa and the computational complexity as ordinate. As shown in Figure 9, the number of calculations can be significantly reduced using the proposed blocking procedure.
For the filtering effect, the whole image was taken as the input, and the same parameters were adopted. As shown in the filtering result in Figure 10, the speckle noise in the homogeneous region of the filtered image has been smoothed, while some riverway boundaries became blurred due to some strong speckle noises in the image. If these strong noises are filtered out during global processing, the weak noise region must be oversmoothed. When the image is partitioned, the background in each sub-image block becomes relatively simple, and the noise intensity is relatively consistent, effectively alleviating the problem. minimize image noise. For the computational complexity of SRAD, the c(q) at (i ± 1, j ± 1) needs to be calculated when smoothing points (i, j). From Equation (4), the process involves square operation, which means that the computational complexity of SRAD filtering (Figure 4) in each iteration would not be less than (2 × M × N) 2 . The image was processed using the proposed blocking strategy, and the study area was partitioned into 16 sub-image blocks (see Figure 5). The computational complexity of the SRAD filter was less than 16(2 × M/4 × N/4) 2 , as the sub-image blocks that do not cover riverways are only used as background and not for filtering. The results were then graphed and compared using the number of pixels as abscissa and the computational complexity as ordinate. As shown in Figure 9, the number of calculations can be significantly reduced using the proposed blocking procedure. For the filtering effect, the whole image was taken as the input, and the same parameters were adopted. As shown in the filtering result in Figure 10, the speckle noise in the homogeneous region of the filtered image has been smoothed, while some riverway boundaries became blurred due to some strong speckle noises in the image. If these strong noises are filtered out during global processing, the weak noise region must be oversmoothed. When the image is partitioned, the background in each sub-image block becomes relatively simple, and the noise intensity is relatively consistent, effectively alleviating the problem. For the boundary effect of block partition, it is due to image degradation that occurs at the boundary of the block when the sub-image blocks are processed, which makes the objects occupying adjacent sub-image blocks appear dislocation or even discontinuous in the stitching result. Therefore, this paper adopts the sub-image block overlap partition strategy. The adjacent sub-image blocks (B1,2 and B1,3) in the non-overlapping partition and the adjacent sub-image blocks (B'1,2 and B'1,3) in the corresponding overlapping partition are selected as an example, as shown in Figure 11(a1)-(d1), where B'1,2 and B'1,3 share a region of 10 pixels wide at the adjacent boundary relative to B1,2 and B1,3. After SRAD filtering, riverway segments extraction and discontinuity connection, the extraction results are stitched, as shown in Figure 11(a2),(b2), and the stitched For the boundary effect of block partition, it is due to image degradation that occurs at the boundary of the block when the sub-image blocks are processed, which makes the objects occupying adjacent sub-image blocks appear dislocation or even discontinuous in the stitching result. Therefore, this paper adopts the sub-image block overlap partition strategy. The adjacent sub-image blocks (B 1,2 and B 1,3 ) in the non-overlapping partition and the adjacent sub-image blocks (B' 1,2 and B' 1,3 ) in the corresponding overlapping partition are selected as an example, as shown in Figure 11a1-d1,  where B' 1,2 and B' 1,3 share a region of 10 pixels wide at the adjacent boundary relative to B 1,2 and B 1,3 . After SRAD filtering, riverway segments extraction and discontinuity connection, the extraction results are stitched, as shown in Figure 11a2,b2, and the stitched boundary parts are marked by a red rectangle box. To make it clearer, enlarge the red rectangle as shown in Figure 11a3 For the overlap degree v0 and h0, selecting a larger value can ensure smooth connection, but it will increase the amount of calculation. As the sub-image block extraction result is a binary image, it will not involve the color difference of different categories, so we select an overlap of 5 to 10 pixels that can meet the stitching requirements. The overlap degree can also be adjusted adaptively according to the size of the block, such as setting it to 5% of the block size.
The impact of the step of classifying the sub-image blocks by visual interpretation on the final result is verified by processing all sub-image blocks. According to the previous classification, each sub-image block in the set B'l 2 = {B'2,1, B'2,2, B'3,1, B'4,3, B'4,4} is processed by the extraction procedure of this paper, and the results are shown in Figure 12(a1)-(e1). As can be seen from Figure  12(a1)-(e1), there are regions identified as riverways in each sub-image block. This is because the extraction is based on the regions with low overall brightness and continuous distribution, while there are always some regions with relatively low brightness in the sub-image blocks. It can be seen from the comparison with the original image in Figure 4 that these mis-extractions are low brightness object regions in the sub-image blocks, such as asphalt road, shadow, urban artificial lake, sewage pool and so on. All the extraction results are stitched together, as shown in Figure  12(a2). It can be seen from the figure that these mis-extractions are not connected to the real riverway, while those originally interrupted riverway segments are connected together after the treatment of the discontinuity connection measures proposed in this paper. In the last step of this method, the connected component processing of Figure 12(a2) can remove the mis-extractions, and the result is shown in Figure 12(b2).
By processing all the sub-image blocks, it can be found that classifying the sub-image blocks in It can be seen from the comparison that Figure 11a3 has dislocation at the stitching line and the connection is not smooth. Figure 11b3 can discard the boundary degradation part through the shared region between B' 1,2 and B' 1,3 during stitching, so that the object can naturally transition at the stitching line.
For the overlap degree v 0 and h 0 , selecting a larger value can ensure smooth connection, but it will increase the amount of calculation. As the sub-image block extraction result is a binary image, it will not involve the color difference of different categories, so we select an overlap of 5 to 10 pixels that can meet the stitching requirements. The overlap degree can also be adjusted adaptively according to the size of the block, such as setting it to 5% of the block size.
The impact of the step of classifying the sub-image blocks by visual interpretation on the final result is verified by processing all sub-image blocks. According to the previous classification, each sub-image block in the set B' l 2 = {B' 2,1 , B' 2,2 , B' 3,1 , B' 4,3 , B' 4,4 } is processed by the extraction procedure of this paper, and the results are shown in Figure 12a1-e1. As can be seen from Figure 12a1-e1, there are regions identified as riverways in each sub-image block. This is because the extraction is based on the regions with low overall brightness and continuous distribution, while there are always some regions with relatively low brightness in the sub-image blocks. It can be seen from the comparison with the original image in Figure 4 that these mis-extractions are low brightness object regions in the sub-image blocks, such as asphalt road, shadow, urban artificial lake, sewage pool and so on. All the extraction results are stitched together, as shown in Figure 12a2. It can be seen from the figure that these mis-extractions are not connected to the real riverway, while those originally interrupted riverway segments are connected together after the treatment of the discontinuity connection measures proposed in this paper. In the last step of this method, the connected component processing of Figure 12a2 can remove the mis-extractions, and the result is shown in Figure 12b2.
Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 27 characteristic of riverway connectivity and the operation measure of discontinuity connection proposed in this paper, so that the mis-extracted regions can be removed in the final treatment.

Sub-image Block Filtering
To illustrate the advantages of the proposed sub-image block filtering, Lee, Kuan, and Frost filtering methods were used as comparisons, and all methods use the iterative criteria used in this paper. The filtering results are shown in Figure 13. By processing all the sub-image blocks, it can be found that classifying the sub-image blocks in advance can reduce the amount of unnecessary calculations. When human judgment is wrong, these sub-image blocks that do not cover the riverway only increase the amount of calculations, but have no effect on the accuracy of the final extraction results. This is mainly due to the natural characteristic of riverway connectivity and the operation measure of discontinuity connection proposed in this paper, so that the mis-extracted regions can be removed in the final treatment.

Sub-Image Block Filtering
To illustrate the advantages of the proposed sub-image block filtering, Lee, Kuan, and Frost filtering methods were used as comparisons, and all methods use the iterative criteria used in this paper. The filtering results are shown in Figure 13.
Comparing Figure 13 with Figure 6a1-c1, the comparative analysis shows the Lee filtering approach yielded poor results. Lee filtering is a model-based method, where the static noise assumption does not often meet the actual situation. Compared with the Lee method, Kuan and Frost use adaptive local statistical features to consider the nonuniformity of the image. Their results show much better removal of speckle noise, but many stray points still remain. The adaptive SRAD proposed in this paper can fully smoothen the internal noise of the region and keep the edge information of the water-land junction, which yields a more accurate extraction of riverways.
For the quantitative evaluation of the filtering results, two indexes-the Equivalent Numbers of Looks (ENL) and the Contrast Noise Ratio (CNR)-were selected. ENL is an index that measures regional smoothness. A 128 × 128 pixels background region was selected as ROI1 (red rectangular in Figure 13), and four 60 × 60 pixels object regions were used as ROI2-ROI5 (blue rectangular in Figure 13) to calculate the average value of the ENL: where µ i and σ i are the mean and standard deviation of the i-th ROI, respectively. Larger ENL values mean that the image is smoothed well. CNR, which measures image contrast by dividing the intensity difference between the selected object and the background regions by the sum of standard deviations, can be calculated as The larger the value of CNR, the stronger the separability between the object and the background region, which is more conducive for riverway extraction. The summary of the evaluation results is presented in Table 3. Comparing Figure 13 with Figure 6(a1)-(c1), the comparative analysis shows the Lee filtering approach yielded poor results. Lee filtering is a model-based method, where the static noise assumption does not often meet the actual situation. Compared with the Lee method, Kuan and Frost use adaptive local statistical features to consider the nonuniformity of the image. Their results show much better removal of speckle noise, but many stray points still remain. The adaptive SRAD proposed in this paper can fully smoothen the internal noise of the region and keep the edge information of the water-land junction, which yields a more accurate extraction of riverways.
For the quantitative evaluation of the filtering results, two indexes-the Equivalent Numbers of Looks (ENL) and the Contrast Noise Ratio (CNR)-were selected. ENL is an index that measures regional smoothness. A 128 × 128 pixels background region was selected as ROI1 (red rectangular in Figure 13), and four 60 × 60 pixels object regions were used as ROI2-ROI5 (blue rectangular in Figure 13) to calculate the average value of the ENL: where μi and σi are the mean and standard deviation of the i-th ROI, respectively. Larger ENL values mean that the image is smoothed well. CNR, which measures image contrast by dividing the intensity difference between the selected object and the background regions by the sum of standard  As shown in the table, SRAD has larger ENL and CNR values, which verifies the feasibility and advantages of selecting SRAD for filtering sub-image blocks.

Discontinuity Connection
In view of the complex shape and discontinuity of the urban riverway, the minimum convex hull and multi-layer region growth processing mechanism is designed to automatically connect the discontinuities between riverway segments. Although this mechanism can achieve good connection results, it also has limitations.
First, the minimum convex hull processing mechanism requires a postprocessing to perfect the result. For safety and beautification considerations, some bridges, dams, and other human-made infrastructures are provided with guardrails, lamp posts, drainage, etc. on the surface boundaries, causing them to appear as incompletely homogeneous regions on the image. Due to the heterogeneity of these infrastructures at the boundary, the two ends of the extracted riverway segment may not be straight. The minimum convex hull processing mechanism will fill the concave part of the cross section, resulting in these concave parts will be retained after discontinuity connection. Take B' 4,1 as an example to illustrate the limitation of the convex hull processing mechanism, as shown in Figure 14. It can be inferred from the red box in Figure 14a that the small areas with high brightness on both sides of the bridge should be guardrails and lamp posts. This makes the cross section of riverway segment extracted in Figure 14b uneven, while the convex hull in Figure 14c fills these concave parts. Figure 14d is the result region of growth using multi-layer region growth with convex hull as the connecting object. By superimposing Figure 14d over Figure 14b, as shown in Figure 14e, it can be seen that these concave parts are retained. For this reason, this paper again uses the connected component method in Section 3.2 to remove these concave parts, as shown in Figure  14f.
Second, the parameters of multi-layer region growth cannot be adjusted adaptively. Due to the difference of the distance between the river segments in each sub-image block, when the number of layers of the pyramid is small, some wider discontinuities cannot be fully connected. To reduce the debugging of the layer number setting and adapt to all sub-image blocks, a larger layer number is set in this paper. Two sub-image blocks- B'1,1 and B'1,2-with different size discontinuities are selected as an example to illustrate the limitations of parameter adjustment, as shown in Figure 15.
From the comparison in Figure 15, it can be seen that because the discontinuity width of B'1,1 is It can be inferred from the red box in Figure 14a that the small areas with high brightness on both sides of the bridge should be guardrails and lamp posts. This makes the cross section of riverway segment extracted in Figure 14b uneven, while the convex hull in Figure 14c fills these concave parts. Figure 14d is the result region of growth using multi-layer region growth with convex hull as the connecting object. By superimposing Figure 14d over Figure 14b, as shown in Figure 14e, it can be seen that these concave parts are retained. For this reason, this paper again uses the connected component method in Section 3.2 to remove these concave parts, as shown in Figure 14f.
Second, the parameters of multi-layer region growth cannot be adjusted adaptively. Due to the difference of the distance between the river segments in each sub-image block, when the number of layers of the pyramid is small, some wider discontinuities cannot be fully connected. To reduce the debugging of the layer number setting and adapt to all sub-image blocks, a larger layer number is set in this paper. Two sub-image blocks-B' 1,1 and B' 1,2 -with different size discontinuities are selected as an example to illustrate the limitations of parameter adjustment, as shown in Figure 15. In order to reduce the dependence of layer number setting, the morphological dilation and Fast-R-CNN are used to estimate the boundary of discontinuities, so as to adjust the number of layers of each sub-image block adaptively. However, at present, only the detection of similar discontinuity objects in the image is realized, and its accuracy and robustness need to be improved.

Feasibility and Robustness Analysis
To verify the feasibility and robustness of the proposed method, about 50 SAR urban riverway images obtained from the GF-3 satellites are tested in this section. Limited by the layout of the article, four images with different types of urban riverways are selected for display and analysis. The strip imaging modes of these images are UFS, FSII, FSI, and UFS, with sizes of 1100 × 1800, 3000 × 1950, 1400 × 1000, and 1550 × 2950 pixels, respectively, as shown in Figure 16a-d.
From the selected images, it can be seen that the urban riverways have complex geometric shapes such as transverse, longitudinal and crisscross, and so on due to the joint action of nature and human-made, and the span range is large, which makes the background more complex and increases the difficulty of accurate extraction of riverway. To this end, this paper partitions each image into overlapping sub-image blocks. The partitioning strategy can artificially specify the size of or the number of sub-image blocks. For example, in Figure 16a-d, the numbers of sub-image blocks are specified as 6, 12, 8, and 15, respectively; h0 = v0 = 5; and the partition results are displayed as red rectangular boxes in the figures. It can be seen from the partition that some sub-image blocks do not cover riverway and can be directly processed as the background, while the background of the sub-image blocks covering the riverway is relatively simple, which is beneficial to the extraction of the riverway. Each image is processed by the proposed method, and the riverway extraction result is represented as red and superimposed on the original image, as shown in Figure 17. From the comparison in Figure 15, it can be seen that because the discontinuity width of B' 1,1 is small, when the 3-layer search is set, the discontinuity connection can be completed, while the discontinuity of B' 1,2 is only partially connected. B' 1,2 can complete the discontinuity connection when the 5-layer search is set, while B' 1,1 has stopped the connection.
In order to reduce the dependence of layer number setting, the morphological dilation and Fast-R-CNN are used to estimate the boundary of discontinuities, so as to adjust the number of layers of each sub-image block adaptively. However, at present, only the detection of similar discontinuity objects in the image is realized, and its accuracy and robustness need to be improved.

Feasibility and Robustness Analysis
To verify the feasibility and robustness of the proposed method, about 50 SAR urban riverway images obtained from the GF-3 satellites are tested in this section. Limited by the layout of the article, four images with different types of urban riverways are selected for display and analysis. The strip imaging modes of these images are UFS, FSII, FSI, and UFS, with sizes of 1100 × 1800, 3000 × 1950, 1400 × 1000, and 1550 × 2950 pixels, respectively, as shown in Figure 16a-d.
From the selected images, it can be seen that the urban riverways have complex geometric shapes such as transverse, longitudinal and crisscross, and so on due to the joint action of nature and human-made, and the span range is large, which makes the background more complex and increases the difficulty of accurate extraction of riverway. To this end, this paper partitions each image into overlapping sub-image blocks. The partitioning strategy can artificially specify the size of or the number of sub-image blocks. For example, in Figure 16a-d, the numbers of sub-image blocks are specified as 6, 12, 8, and 15, respectively; h 0 = v 0 = 5; and the partition results are displayed as red rectangular boxes in the figures. It can be seen from the partition that some sub-image blocks do not cover riverway and can be directly processed as the background, while the background of the sub-image blocks covering the riverway is relatively simple, which is beneficial to the extraction of the riverway. Each image is processed by the proposed method, and the riverway extraction result is represented as red and superimposed on the original image, as shown in Figure 17. Remote Sens. 2020, 12, x FOR PEER REVIEW 24 of 27 From the superposition of the results, it can be seen that the riverway extracted from each test image has good visual accuracy in terms of completeness and boundary positioning. The quantitative evaluation of each image is shown in Table 4.  From the superposition of the results, it can be seen that the riverway extracted from each test image has good visual accuracy in terms of completeness and boundary positioning. The quantitative evaluation of each image is shown in Table 4.  As can be seen from the table, the average Sorensen-Dice and Jaccard scores of images are 93.06% and 87.05%, respectively, and the boundary accuracy has reached 85% on average when the radius is 2 pixels. Based on the qualitative and quantitative evaluation of the above extraction results, the feasibility and robustness of the proposed method are verified.

f5. Conclusions
In this paper, a method of urban riverway extraction from high-resolution SAR images is proposed. According to the distribution characteristics of urban riverways, this method adopts the sub-image block extraction processing strategy. Sub-image blocks that do not cover riverways are As can be seen from the table, the average Sorensen-Dice and Jaccard scores of images are 93.06% and 87.05%, respectively, and the boundary accuracy has reached 85% on average when the radius is 2 pixels. Based on the qualitative and quantitative evaluation of the above extraction results, the feasibility and robustness of the proposed method are verified.

Conclusions
In this paper, a method of urban riverway extraction from high-resolution SAR images is proposed. According to the distribution characteristics of urban riverways, this method adopts the sub-image block extraction processing strategy. Sub-image blocks that do not cover riverways are treated as background. For those that do cover riverways, the designed extraction method first uses an iterative adaptive SRAD and optimized Sauvola local threshold algorithm to quickly obtain the riverway segments in sub-image blocks. The minimum convex hull of each riverway segment is then taken as the connection object, and the discontinuity connection is realized by multi-layer region growth.
The qualitative and quantitative evaluations of the riverway extraction results show the effectiveness of the proposed method. However, for different sub-image blocks, some operations in the extraction method would have to be adjusted and selected manually, reducing the automated nature of the extraction method. For future research, the adaptive adjustment of the extraction parameters for each sub-image block need to be studied further to achieve automatic extraction of urban riverways from large-scale high-resolution SAR images.