Offshore Bridge Detection in Polarimetric SAR Images Based on Water Network Construction Using Markov Tree

Abstract

It is difficult to detect bridges in synthetic aperture radar (SAR) images due to the inherent speckle noise of SAR images, the interference generated by strong coastal scatterers, and the diversity of bridge and coastal terrain morphologies. In this paper, we present a two-step bridge detection method for polarimetric SAR imagery, in which the probability graph model of a Markov tree is used to build the water network, and bridges are detected by traversing the graph of the water network to determine all adjacent water branch pairs. In the step of the water network construction, candidate water branches are first extracted by using a region-based level set segmentation method. The water network is then built globally as a tree by connecting the extracted water branches based on the probabilistic graph model of a Markov tree, in which a node denotes a single branch and an edge denotes the connection of two adjacent branches. In the step of the bridge detection, all adjacent water branch pairs related to bridges are searched by traversing the constructed tree. Each bridge is finally detected by merging the two contours of the corresponding branch pair. Three polarimetric SAR data acquired by RADARSAT-2 covering Singapore and Lingshui, China, and by TerraSAR-X covering Singapore, are used for testing. The experimental results show that the detection rate, the false alarm rate, and the intersection over union (IoU) between the recognized bridge body and the ground truth are all improved by using the proposed method, compared to the method that constructs a water network based on water branches merging by contour distance.

1. Introduction

Bridge detection in synthetic aperture radar (SAR) images is essential in applications such as autonomous navigation, urban monitoring and planning, flood monitoring, and disaster prevention, among others. With respect to geometric structures, a bridge is a long strip-shaped area that spans a water tributary connecting two land regions. Usually, the two sides of a bridge are straight and parallel to each other. In term of physical scattering characteristics, the double-bounce scattering generated from various components of man-made structures on bridges make large bridges show as bright regions in SAR images. As shown in Figure 1, some bridges marked in red circles are located on a single tributary. Because the scattering intensity of water is much lower than that of land in SAR images, bridges seem easy to detect by extracting water tributaries and finding the connecting land regions on the tributaries. However, it is not a simple procedure, as the imaging geometry may play a factor in significantly reducing the radar returns from the bridge so that it is difficult to discriminate between the bridge and speckle noise in the water (region 1 in Figure 1). Additionally, the interference caused by strong scatterers along the coast and the diversity of bridge and coastal terrain morphologies [1] (region 2 in Figure 1) are also influencing factors.

In almost all existing bridge detection methods, water branches are extracted by sea-land segmentation or classification first before bridges are detected by geometric and scattering feature extraction and recognition. Luo et al. [2] extracted water regions using Gauss Markov random fields (MRFs) in combination with support vector machine methods and detected bridges by extract line across the water area based on the prior knowledge of the bridge width. Chaudhuri et al. [3] extracted river branches by classifying the image into one of three types, including water, building, and background, and recursively scanning the extracted water. Bridges were then detected based on the spatial relationship of the bridge and river. Wang et al. [4] extracted water regions using an edge detection method and detected bridges by the spatial relationship of the bridge and water area. Chen et al. [5] tracked water contours using a particle filter and detected bridges by scanning the binary segmentation result of water regions. Chen et al. [6] extracted water regions using a global thresholding segmentation method and detected bridges by a combination of morphological processing and the geometric model of the bridge. Because the geometric feature of the bridge is salient, the methods based on spatial relationships and geometric structures can easily detect bridges. However, the performance of the methods depends on the accuracies of the water–land segmentation and the geometric feature extraction. It is error-prone in SAR images due to strong coherent speckle and scattering interference. Except the methods based on the geometric structure and spatial distribution of bridges, Wang et al. [7] extracted and identified the region of interest (ROI) of bridges by a combination of the context information, the scattering intensity and the spatial structure of the bridge in the SAR image. Song et al. [8] detect bridges using a constant false alarm rate (CFAR) detector based on the high intensity of the bridge in the SAR image. The problem with the intensity-based method is that only some bridges are strong scatterers. In polarimetric SAR imagery, Yu et al. [9] first extracted water using a polarimetric parameter to enhance the contrast between bridges and water. Bridges were then detected by examining the side lines of the bridge extracted by a ratio edge detector.

Accurate extraction and connection of tributaries of water are crucial steps for bridge detection. When all the branches of water area are accurately extracted, we only need to determine the small land areas that connect different water branches to detect bridges. In the existing methods, water and land are segmented only by using thresholding segmentation or MRF segmentation methods. It is not applied when the water branches are so complex that they form a network. For the extraction of a complex water network, the existing methods are mainly based on digital elevation model (DEM) images and optical images. Fern et al. [10] first proposed extracting a water network by building a local network structure model using the elevation of water. Bai et al. [11] extracted and connected the water network by building a binary search tree model using DEM data. In the method based on optical remote sensing images, the water area is first extracted using the global thresholding segmentation method with a normalized difference water index (NDWI) [12] parameter. The water network is then constructed by extracting the centerline of the water branches using morphological processing and merging the centerlines by distance. Isikdogan et al. [13] extracted water regions with a modified NWDI index, then determined water branches using a proposed curve singularity index, and finally built the water network by tracing the centerline of the water branches by a non-maximum suppression. Yang et al. [14] also extracted water regions using the modified NWDI index, but connected the water component regions using a multi-point marching method. Chen et al. [15] proposed a supervised water network extraction method based on multi-temporal data, which counts the frequency of the water regions that appear in different time phases. The main and central lines of the water network were obtained from water in a drought period, whereas the small branches were extracted from water in a flood period. Isikdogan et al. [16] proposed a water network extraction method based on a supervised deep neural network. However, the methods based on DEM and optical images are not applied to SAR images due to the strong speckle and scattering interference. In SAR and polarimetric SAR images, Klemenjak et al. [17] proposed a method based on the morphological processing and a supervised classification method. Liu et al. [18,19] proposed two water and land segmentation methods using multi-level and multi-scale level set segmentation. Liu et al. [20] segmented the water and land using a proposed volume scattering polarimetric parameter to perform the offshore port detection. Liu et al. [21] proposed an improved level set segmentation method by limiting the initial segmentation to segment the water and offshore oil platforms. Obida et al. [22] proposed a water extraction method based on a combination of DEM and SAR data. However, the existing SAR-based methods are not applied to the extraction of a complex water network. The problem of mixing DEM and SAR data is that the DEM data are often inaccurate or incomplete [22], particularly when observing waterways, making it difficult to fuse the two data types for an effective water detection technique.

If we represent a water branch as a set of centerlines with constant width, the connection of water branches is similar to a network of roads. Tupin et al. [23] connected the candidate sets of road lines extracted by edge detection using a Markov graph model. Krishnamachari et al. [24] defined energy terms related to distance and direction for the extracted edge lines and extracted the boundaries of buildings using an MRF model. Tu et al. [25,26] first defined the energy terms of a single curve branch, parallel curve branch and curve group with tree structure, respectively. Images were then parsed to region, curve, and curve group using the Markov Monte Carlo method. However, when using the existing line segments connection model, it is difficult to achieve high accuracy in SAR images due to the strong speckle and scattering interference.

There are two main challenges for offshore bridge detection in SAR images. One is some branches are too narrow and too small to be correctly extracted. The other is the scatterings of some branches are too ambiguous due to the scattering interference caused by coastal building to be correctly segmented. If some branches are wrongly classified as land, the tributary that contains these branches will be broken for a long distance and cannot be correctly connected. However, the whole topology of those water branches is usually very similar to a tree. Thus, to improve the performance of the extraction and connection of water branches, this paper proposes an offshore bridge detection method in SAR and polarimetric SAR imagery by connecting water networks based on a Markov tree. Unlike traditional methods directly connecting water branches by the local spatial relationship of branches, the proposed method globally constructs the water network by a tree using the probabilistic graph model of Markov, in which branches are taken as nodes and the adjacent branches are taken as edges. Based on the tree structure of branches, the ROIs of bridges are easy to extract by traversing the tree to determine all the edges. Bridges are finally detected by merging the feature points of the two contours of each adjacent branch pair.

The rest of this paper is organized as follows. In Section 2, the general segmentation model of the water network is introduced. In Section 3, the proposed method is introduced in detail. Experimental results are shown and discussed in Section 4. The discussion is given in Section 5. The conclusion is given in Section 6. The nomenclature of all symbols in the article is shown in Appendix A. The structure diagram of all symbols used in the water network construction is shown in Appendix B.

2. General Segmentation Model of the Water Network

Bridges cross over water, and the minimum distance between two water branches adjacent to a bridge is near equal to the width of the bridge. To detect bridges, we only need to extract water branches and connect the branches to determine the break part of the tributary of the branches. For a coastal image, the land is composed of several homogeneous regions, and the water network is composed of inland water regions, a main water trunk, and a number of tributaries. The shape of the water network is so similar to a tree that the topological structure can be represented by a tree structure, in which a node denotes a single branch, the main trunk denotes the root of the tree, and an edge denotes two branches are adjacent. Thus, the image can be modeled as a set including some homogeneous regions, some single branches, some water tributaries with a tree structure, and the topological structure of different regions. By building the probability model of each component of the set, we can get the maximum a posteriori probability (MAP) model for the image parse based on Bayesian criterion.

2.1. General Model of the Global Segmentation

If $\mathit{R}$ denotes the image plane, and the given polarimetric SAR image is represented as $\mathbf{T}$, the objective of image segmentation is to parse the image into a set of terms, including some homogeneous regions $R^r$, free curve branches $R^c$, tree regions $R^{\psi }$, and the topological relationship $TP$. If the parameter set of the segmentation is represented as $\Theta =\left(R^r,R^c,R^{\psi },TP\right)$ and $\Omega$ denotes the parameter space, the MAP model can be represented as

$${\Theta }^{*}={}_{\Theta \in \Omega }^{argmax}-\log \left(p\left(\mathbf{T}|\Theta \right)p\left(\Theta \right)\right),$$

(1)

where $p\left(\Theta \right)$ denotes the prior probability of the model, and it is related to the prior model of different terms; $p\left(\mathbf{T}|\Theta \right)$ denotes the conditional likelihood probability distribution of the coherent matrices, and it is related to the probability distribution of regions.

2.2. Probability Distribution of a Region

If both regions $R^r$ and all curve branches are homogeneous, the coherent matrix of a homogeneous region obeys the complex Wishart distribution in polarimetric SAR imagery (for SAR images, the intensity of a homogeneous branch obeys a Gamma distribution). If the average coherent matrix is $\Sigma$ and the number of looks is L, then the coherent matrix $\mathit{T}\sim W\left(\Sigma ,L,p\right)$ [27],

$$f\left(\mathit{T}|\Sigma ,L,p\right)=\frac{L^{pL}{\left|\mathit{T}\right|}^{L-p}\exp \left\{-L\mathrm{tr}\left({\Sigma }^{-1}\mathit{T}\right)\right\}}{K\left(L,p\right){\left(\left|\Sigma \right|\right)}^L},$$

(2)

where p denotes the number of polarimetric channels; $\mathrm{tr}\left(.\right)$ denotes the trace of matrix; $K\left(L,p\right)={\pi }^{p(p-1)/2}\Gamma \left(L\right)\dots \Gamma \left(L-p+1\right)$; and $\Gamma \left(.\right)$ is the Gamma function.

The segmentation regions are limited by the length of the contour and the area of the region in prior. Thus, the prior probability of regions satisfies

$$p\left(R\right)\propto \exp \left\{-{\gamma }_r{\left|D\left(R\right)\right|}^{\rho }-{\lambda }_r\left|\partial R\right|\right\},$$

(3)

where ${\gamma }_r,\rho ,{\lambda }_r$ are constant, $\left|D\left(R\right)\right|$ denotes the area of region R, and $\left|\partial R\right|$ denotes the contour length of region R.

2.3. Model of a Single Water Branch

Similar to reference [25], a single water branch is represented as a long strip area with a curve as its centerline, and its width is nearly constant. Based on the probability distribution of the coherent matrix in the homogeneous region (2), if the average coherent matrix is ${\Sigma }_j$, and the coherent matrix of a single water area obeys $\mathit{T}\sim W\left({\Sigma }_j,L,p\right)$. Given curve branch C, if its whole length is S and the arc coordinates are $s=\left[0,S\right]$, it can be represented as a region $D\left(C\right)$ where the center curve is $\mathit{c}\left(s\right)=\left(x_m\left(s\right),y_m\left(s\right)\right)$ and the width is $w\left(s\right)$ for region $D\left(C\right)$. The value of $w\left(s\right)$ is almost invariant in the region and $w\left(s\right)\ll S$. Then, the prior probability of the curve branch C is

$$p\left(C\right)\propto p\left(D\left(C\right)\right)p\left(\mathit{c}\left(s\right)\right)p\left(w\left(s\right)\right)\propto \exp \left\{-E\left(C\right)\right\}.$$

(4)

As $E\left(C\right)$ depends on the area, length, and width of the region, we define it as

$$E\left(C\right)={\gamma }_c{\left|D\left(C\right)\right|}^{\rho }-{\lambda }_{c}S+{{\mu }_{c}E}_o\left(w\right),$$

(5)

where ${\gamma }_c,\rho ,{\lambda }_c,{\mu }_c$ are constant; $\left|D\left(C\right)\right|$ denote; the area of region $D\left(C\right)$; $E_o\left(w\right)$ is a function of the consistency of $w\left(s\right)$. If $w\left(s\right)$ is high, then $E_o\left(w\right)$ is small. The expression ${\gamma }_c{\left|D\left(C\right)\right|}^{\rho }$ denotes the energy related to the area, ${\lambda }_{c}S$ denotes the energy related to the length, and ${{\mu }_{c}E}_o\left(w\right)$ denotes the energy related to the consistency. The energy $E\left(C\right)$ limits the water branch to be a region with small area, large length, and high consistency.

2.4. Model of Water Tributaries with a Tree Structure

Because the region of the water tributaries with a tree structure is homogeneous, its coherent matrix satisfies $\mathit{T}\sim W\left(\Sigma ,L,p\right)$. Given a tree tributary $\psi$, if it consists of n single branches $\left\{C_1,C_2,\dots ,C_n\right\}$ that denote the node of the tree, the structure relationship $TP$ formed by the branches is $〈{\alpha }_1,{\beta }_k〉,〈{\alpha }_k,{\beta }_l〉,\dots ,〈{\alpha }_l,{\beta }_n〉$. Both ${\alpha }_i$ and ${\beta }_i$ are the index of branch $C_i$; $〈{\alpha }_k,{\beta }_l〉$ denotes the edge of the tree such that ${\alpha }_k$ is the father node and ${\beta }_l$ is the child node. Then we have

$$\psi =\left(n,\left\{C_1,C_2,\dots ,C_n\right\},TP\right).$$

(6)

Supposing that $\psi$ satisfies the hypothesis of a Markov field that the label of a single node is only determined by the labels of the adjacent parent and children nodes, then $p\left(\psi \right)$ obeys the Gibbs distribution according to the Hammersley Clifford theorem [28]:

$$p\left(\psi \right)=\frac{1}{Z}\prod _{c\in C_{\psi }}\phi \left(X_c\right)\propto \exp \left\{-\sum _{c\in C_{\psi }}E\left(X_c\right)\right\},$$

(7)

where $C_{\psi }$ is the set of a series of cliques in graph $\psi$, and Z is the normalization constant.

From (7), we get the prior probability of $\psi$

$$p\left(\psi \right)\propto \exp \left\{-{\lambda }_{\psi }n-\sum _{i=1}^{n}E\left(C_i\right)-\sum _{{\alpha }_i\ne \varnothing }E\left(C_i,C_{{\beta }_i}\right)\right\},$$

(8)

where $E\left(C_i\right)$ denotes the energy of branch $C_i$, $E\left(C_i,C_{{\beta }_i}\right)$ denotes the energy of the branch pair $C_i$ and $C_{{\beta }_i}$, and ${\sum }_{{\alpha }_i\ne \varnothing }E\left(C_i,C_{{\beta }_i}\right)$ denotes the total energy of all branches adjacent to $C_i$; ${\lambda }_{\psi }n$ is the regularization term related to the number of branches.

According to (3), (4) and (8), if the number of regions $R^r$ is $N^r$, the number of single branches is $N^c$, and the number of trees is $N^{\psi }$, the prior probability of image $\mathbf{R}$ is

$$p\left(\Theta \right)=\prod _{i=1}^{N^r}p\left(R_i\right)\prod _{j=1}^{N^c}p\left(C_j\right)\prod _{k=1}^{N^{\psi }}p\left({\psi }_k\right).$$

(9)

Because a free single branch can be considered as a special tree with one node, if the total number of trees including single branches and branches in tree structures is M, the likelihood probability is

$$\begin{array}{cc}\hfill p\left(\mathbf{T}|\Theta \right)=& \prod _{i=1}^{N^r}\prod _{\left(x,y\right)\in R^{r_i}}f\left({\mathit{T}}_i\left(x,y\right)|{\Sigma }_i,L,p\right)\hfill \\ \hfill & \prod _{j=1}^M\prod _{\left(x,y\right)\in R^{c_j}}f\left({\mathit{T}}_j\left(x,y\right)|{\Sigma }_j,L,p\right),\hfill \end{array}$$

(10)

where ${\mathit{T}}_i\left(x,y\right)$ denotes the coherent matrix of point $\left(x,y\right)$.

Based on the prior and likelihood probability Formulas (9) and (10) as well as the prior model Formulas (3), (4) and (8), we get the final general model of the segmentation from (1).

3. The Proposed Method

3.1. Overview of the Proposed Method

The optimum parameter $\Theta$ can be searched using the Markov Monte Carlo method based on (1). However, the complexity of the algorithm is very high because the parameter space is huge. To reduce the complexity of the general model, a hierarchical model is proposed to divide the model into two layers: one is water branch extraction by segmentation, and the other is a water network connection, as shown in Figure 2b. The general model of the segmentation presented in Section 2 can be represented in Figure 2a, where the latent random variables $R,C,\psi ,TP$ denote the four parameters of the model $\left(R^r,R^c,R^{\psi },TP\right)$, and variable $\mathit{T}$ denotes the observed coherent matrices. The hierarchical model is shown in Figure 2b, where the latent variable K denotes whether the water branches belong to the water network, and the latent variable Y denotes whether the segmented regions belong to the land or the water. The hierarchical model can be greedily solved layer by layer based on each graph.

The flowchart of the proposed method is shown in Figure 3. The first part of the algorithm performs the water network construction. The first step of this part performs the water branch extraction, which corresponds to the first layer of the hierarchical model. The second step performs the water branch connection, which corresponds to the second layer of the hierarchical model. The second part of the algorithm performs the bridge detection and bridge body extraction.

The first part of the algorithm is described in Section 3.2 and Section 3.3. The steps of the water branch extraction are shown in box a, where the water and land are segmented using a level set method, and all water branches are extracted by a post-process of the segmentation result. The steps of the water branch connection are shown in box b, where the branches are scanned, tracked, and connected recursively to construct the water network. First, the centerline, endpoints, and directions of each branch are extracted by a morphological processing in a combination of the spline interpolation. Next, the main water trunk and its bifurcations are extracted from the branches using a binary image scanning method. Then, the initial graph structure of a local tributary takes a bifurcation as a root is constructed by connecting branches based on the distance and direction angle of different centerlines. The initial whole graph structure of the water is constructed by connecting different local tributaries. Finally, by defining the node and edge energy of the graph, the optimal Markov tree can be dynamically solved using the simulated annealing algorithm.

The second part of the algorithm is described in Section 3.4. In this part, the adjacent branch pairs for the bridge detection are extracted by traversing the edges of the tree first. The feature points of each branch contour are then extracted using a curve polygon approximation algorithm. Finally, merging the feature points on the contours of the branch pair by distance, the bridge body is extracted to finish the bridge detection.

3.2. Water Branch Extraction

If curve $\Gamma$ denotes the boundary of the water and land segmentation, then the objective of the two-region segmentation is as follows according to (1):

$${\Gamma }^{*}={}_{\Gamma \in \Omega }^{argmax}-\log \left(p\left(\mathbf{T}|\Gamma \right)p\left(\Gamma \right)\right),$$

(11)

where $p\left(\mathbf{T}|\Gamma \right)$ denotes the probability of each region given the segmentation $\Gamma$, and $p\left(\Gamma \right)$ denotes the prior probability of $\Gamma$.

According to (3), (10) and (11), the objective energy function of the segmentation can be represented as

$$\begin{array}{cc}\hfill E\left(\Gamma \right)=& {\gamma }_r{\left|D\left(R_1\right)\right|}^{\rho }+{\lambda }_r\left|\Gamma \right|-{\int }_{R_1}\log f\left({\mathit{T}}_1\left(x,y\right)|R_1\right)dxdy\hfill \\ \hfill & -{\int }_{R_2}\log f\left({\mathit{T}}_2\left(x,y\right)|R_2\right)dxdy,\hfill \end{array}$$

(12)

where $R_1$ is the region inside the curve $\Gamma$, $R_2$ is the region outside the curve $\Gamma$, and $f\left({\mathit{T}}_i\left(x,y\right)|R_i\right)$ is the probability distribution function of region $R_i\left(i\in \left\{1,2\right\}\right)$.

If the curve $\Gamma$ is implicitly embedded into a level set function $\Phi$, the optimal segmentation can be solved by iteratively evolving $\Phi$ as the derivative of the energy function to $\Phi$ [29]. When the segmentation is finished, the water and land regions can be identified by comparing the average scattering total power of the two segmented regions. After the final binary result of the land and water segmentation $\mathit{B}$ is obtained, all water branches $\mathit{D}=\left\{d_j\right\}\left(j\in \left\{1,\dots ,N^d\right\}\right)$ ($N^d$ denotes the number of branches) can be extracted by determining all connection regions from the water segmentation result.

3.3. Water Branch Connection

If $\mathit{W}$ denotes the water network, the class labels of the branches $\mathit{D}$ are $\mathit{K}=\left\{k_j\right\}\left(j\in \left\{0,1\right\}\right)$, where $k_j=1$ means that $d_j$ belongs to $\mathit{W}$, and the value space of $\mathit{W}$ is $\mathbf{W}$, then the general connection model is

$${\mathit{K}}^{*}={}_{\mathit{K}\in \mathbf{W}}^{argmax}-\log \left(p\left(\mathit{D}|\mathit{K}\right)p\left(\mathit{K}\right)\right).$$

(13)

According to (9), the prior probability of labels $\mathit{K}$ is

$$p\left(\mathit{K}\right)=\prod _{i=1}^{N^c}p\left(C_i\right)\prod _{j=1}^{N^{\psi }}p\left({\psi }_j\right),$$

(14)

where $\mathit{D}$ is composed of several single water branches $C_i\left(i\in \left\{1,\dots ,N^c\right\}\right)$ and water tributaries with tree structure ${\psi }_j\left(j\in \left\{1,\dots ,N^{\psi }\right\}\right)$.

Because a single water branch can be considered as a water tributary with a tree structure only containing a root node, formula (14) is equivalent to

$$p\left(\mathit{K}\right)=\prod _{i=1}^{M}p\left({\psi }_i\right)\propto \exp \left\{-\sum _{i=1}^{M}E\left({\mathit{K}}_i\right)\right\},$$

(15)

where M is the sum of $N^c$ and $N^{\psi }$; ${\mathit{K}}_i$ is the label vector of tributary ${\psi }_i$; and $E\left({\mathit{K}}_i\right)$ is the prior energy of tributary ${\psi }_i$, which can be computed from Formula (8).

According to (10), the likelihood probability model of labels $\mathit{K}$ is

$$p\left(\mathit{D}|\mathit{K}\right)=\prod _{i=1}^{M}f\left({\mathit{T}}_i\left({\mathit{D}}_i\right)|{\psi }_i\right)\propto \exp \left\{-\sum _{i=1}^{M}E\left({\mathit{D}}_i|{\mathit{K}}_i\right)\right\},$$

(16)

where $f\left({\mathit{T}}_i\left({\mathit{D}}_i\right)|{\psi }_i\right)$ denotes the likelihood probability of the average coherent matrix ${\mathit{T}}_i\left({\mathit{D}}_i\right)$ of region ${\mathit{D}}_i$ in case of tributary ${\psi }_i$, and $E\left({\mathit{D}}_i|{\mathit{K}}_i\right)$ is the likelihood energy of the tributary ${\psi }_i$.

Similarly, joining (15) and (16), the objective energy function of water connection can be represented as

$${\mathit{K}}^{*}={}_{\mathit{K}\in \mathbf{W}}^{argmin}\sum _{i=1}^M\left(E\left({\mathit{K}}_i\right)+E\left({\mathit{D}}_i|{\mathit{K}}_i\right)\right).$$

(17)

To construct the whole water network of the branch set $\mathit{D}$ based on (17), we need to get the geometric representation of the branches first. We then construct the initial graph structure of the branches based on the spatial adjacent relationship of different branches. Finally, we can search the optimal graph structure by defining the prior and connection energy of the branch and solving the connection probabilistic graph model.

3.3.1. Geometric Representation of a Branch

From Section 2.3, we know that to represent a branch, the centerline and its endpoints and directions need to be obtained first. It is difficult to extract the centerline because the shape of branch is varied and some branches are winding. To get a reasonable centerline, we first extract the raw skeleton of the branch using morphological processing and then fit the skeleton into a curve using the spline interpolation. If the extracted centerline of a branch is ${\mathit{c}}_m\left(s\right)=\left(x_m\left(s\right),y_m\left(s\right)\right)$, and $s\in \left[0,S\right]$, the two endpoints of the branch are ${\mathit{c}}_m\left(0\right)$ and ${\mathit{c}}_m\left(S\right)$, and the two directions of the branch are the tangent directions of points ${\mathit{c}}_m\left(0\right)$ and ${\mathit{c}}_m\left(S\right)$.

3.3.2. Main Trunk and Bifurcations Extraction

The whole water feature is composed of a trunk, several tributaries, and some isolated branches. If the extracted branches set is $\mathit{D}=\left\{d_j\right\}\left(j\in \left\{1,\dots ,N^d\right\}\right)$, where $N^d$ denotes the number of branches, the branch with the maximum area $d_m$ is selected as the main trunk of the water. There are several bifurcations in the main truck, each of which forms the starting point of a tributary. To determine the adjacent relationship of branches on a tributary, the bifurcations need to be extracted first. The bifurcations are areas protruding from the trunk, flanked by land. Taking a pixel on a bifurcation as the center, the pixels with a certain distance along the two sides of the horizontal or vertical direction are on the land. Thus, we can get the bifurcation by a scanning method. An example segmentation result is shown in Figure 4a, where the region in white is water and the region in black is land. One bifurcation of the main trunk is shown in Figure 4b. We find that the bifurcation is a long strip with land on both sides. Because the direction of the bifurcation is varied, the scanning is carried out in multiple directions at the same time. In this paper, we scanned in two group orthogonal directions, where one group is horizontal and vertical, and the other group is 45 degrees and 135 degrees counterclockwise. As shown in Figure 4b, for pixel P on a bifurcation, the scanning is carried out in four directions. For scanning along the direction $\overrightarrow{\tau }$, if there are land pixels in the range from positive ${\delta }_{rad}$ to negative ${\delta }_{rad}$ with P as the center; P is considered as a point in the bifurcation branch. If P is scanned as a bifurcation pixel in the two orthogonal directions, P is classified as a bifurcation pixel in the group directions. If P is classified as a bifurcation pixel in one of the two groups, P is finally classified as a bifurcation pixel.

Because the coastal terrain is uneven, some convex regions are wrongly scanned as bifurcations. To reduce the false regions, the scanned bifurcations are identified by the ratio of length to width, area, and the consistency of width. If the perimeter of a region is $l_P$ and the area is A, the ratio of length to width ${l_L/l}_W$ can be represented as follows:

$${\rho }_{LW}=\frac{l_L}{l_W}=\frac{\left(l_P+\sqrt{l_P^2-16A}\right)}{\left(l_P-\sqrt{l_P^2-16A}\right)}.$$

(18)

There are some bifurcations with large width and area. As shown in Figure 4a, the tributary marked in the red box is broken into two parts by a bridge. The tributary cannot be extracted as a bifurcation using the scanning method. For this kind of tributary, if a branch in the tributary satisfies that the distance between the contour of the branch and the main water is less than ${\delta }_{ref}$ and its area is greater than ${\Delta }_{ref}$, then it is classified as a bifurcation.

3.3.3. Graph Structure Construction

The N-tree is used to represent the graphical structure of a water tributary. As shown in Figure 5, each node $X_i$ in the N-tree denotes a water branch $d_i$. A single node contains one parent node and $j\left(j\le N\right)$ children nodes. Taking a bifurcation as the root node, a tree is constructed by searching for children nodes according to the distance and direction angle between the centerlines of two branches.

For the branch set $\mathit{D}$, if $d_m$ denotes the extracted main trunk, and $\mathit{H}=\left\{h_x\right\}\left(x\in \left\{1,\dots ,N^h\right\}\right)$ denotes the extracted bifurcations on $d_m$, the graph structure of the water network can be recursively built based on the construction algorithm of the tree. First, the angle and distance between each two branches in $\mathit{D}$ are calculated. We take $d_m$ as the root node. Using a flag set to mark whether the branches have been built, all branches are initially marked as 0 except branch $d_m$. Then, starting with bifurcation $h_x$, the successor nodes of N-tree are searched in distance limitation ${\delta }_{ref}$ and direction limitation ${\theta }_{ref}$. The children nodes are set to root nodes and the successor nodes are iteratively searched until the leaf node is reached. We mark all the nodes in the N-tree that have been built. Finally, by following the same steps, each sub-tree of the root node $d_m$ is successively constructed by starting with a bifurcation that is not searched.

3.3.4. Node and Edge Energy Definition

To construct the graph of branches using the probability model, the prior energy of the branch and the connection energy of two adjacent branches need to be defined according to (8) and (15). As a branch is represented as a node and the connection of two branches is represented as an edge, the prior energy of a branch is called node energy $E_{node}$, and the connection energy is called edge energy $E_{edge}$.

The expression of the node energy of a branch $E_{node}$ is defined in (5). The values of the area $\left|D\left(C\right)\right|$ and length S are easy to calculate based on the geometric representation of the branch. However, the consistency function of a branch $E_o\left(w\right)$ needs to be defined. Given branch $d_j$, if ${\mathit{c}}_i\left(s\right)$ denotes the centerline curve of $d_j$, the left boundary is ${\mathit{c}}_l\left(s\right)$ and right boundary is ${\mathit{c}}_r\left(s\right)$, then the width function $w\left(s\right)$ is equal to $\left|{\mathit{c}}_r\left(s\right)-{\mathit{c}}_l\left(s\right)\right|$. Because $w\left(s\right)$ may gradually increase or decrease as s for a water branch, $E_o\left(w\right)$ is defined as the average of the local width consistency of each point on the centerline. The local width consistency of a point is defined as the difference between its width and the average width of the point in a window region. The average width of each point is calculated using a sliding window. If the average width function is $w_o\left(s\right)$, the width consistency of the branch is defined as follows:

$$E_o\left(w\right)={\left(\frac{1}{\left|{\mathit{c}}_i\left(s\right)\right|}\sum _s{\left(w_o\left(s\right)-w\left(s\right)\right)}^2\right)}^{1/2},$$

(19)

where $\left|{\mathit{c}}_i\left(s\right)\right|$ denotes the length of ${\mathit{c}}_i\left(s\right)$.

The edge energy of two adjacent branches is determined by the distance and direction angle of the branches. Supposing that the distance between node Q and node P is ${\delta }_{PQ}$ and the angle of direction is ${\theta }_{PQ}$, the energy of the edge $PQ$ can be defined as

$$E_{edge}\left(PQ\right)=\epsilon sin{\theta }_{PQ}+\eta {\delta }_{PQ}-E_{ref},$$

(20)

where the term $\epsilon sin{\theta }_{PQ}$ restricts the angle of two connected branches; $\eta {\delta }_{PQ}$ restricts the distance of two connected branches; $E_{ref}$ is the reference edge energy; and $E_{ref}$ is subtracted in $E_{edge}$ is a mean normalization process so as to keep a balance between the energy of $E_{edge}$ and the energy of $E_{node}$. If the reference angle is ${\theta }_{ref}$ and the reference distance is ${\delta }_{ref}$, $E_{ref}$ is $\epsilon sin{\theta }_{ref}+\eta {\delta }_{ref}$.

3.3.5. Markov Tree Construction

To construct the tree using (17), the prior and likelihood energies of the whole tree need to be defined.

(1)

Prior energy of Markov tree

From the N-tree structure, we see that the clique of a N-tree $\psi$ is composed of a node, a parent node, and several children nodes adjacent to the node.

As shown in Figure 5, given node Q, if the father node is P and there are $j\left(j\le N\right)$ children nodes, then the prior energy of Q is the sum of node energy and all edge energies related to Q, and the expression is as follows:

$$E_{prior}\left(Q\right)=E_{node}\left(Q\right)+\frac{1}{2}\left(E_{edge}\left(PQ\right)+\sum _{i\in \mathit{O}}{E}_{edge}\left(Qi\right)\right)$$

(21)

where $\mathit{O}$ is the children nodes set of Q. Because edge $PQ$ is the common edge of two nodes Q and P, the edge energy is only set to half for Q.

Calculating the energy of each node respectively and according to (8), the total energy of N-tree $\psi$ is

$$E_{prior}\left(\psi \right)={\lambda }_{\psi }n+\sum _{Q\in \psi }{E}_{prior}\left(Q\right),$$

(22)

where n is the number of nodes.

(2)

Likelihood energy of Markov tree

Given node Q, if it belongs to tree $\psi$, the likelihood energy is determined by the average coherent matrix ${\Sigma }_{\psi }$ of region $\psi$:

$$E_{like}\left(Q\right)=E\left(d_Q|k_Q\right)=-\log f\left(\overline{\mathit{T}}\left(d_Q\right)|{\Sigma }_{\psi }\right),$$

(23)

where $k_Q$ is the class label of Q; $d_Q$ denotes the region of branch Q; and $\overline{\mathit{T}}\left(d_Q\right)$ is the average coherent matrix of $d_Q$.

Similarly, the total likelihood energy of N-tree is the sum of the likelihood energy of each node:

$$E_{like}\left(\psi \right)=\sum _{Q\in \psi }{E}_{like}\left(Q\right).$$

(24)

The structure diagram of the symbols defined in the model of the water network construction from Section 2, Section 3.2 and Section 3.3 is shown in Appendix B.

(3)

Network construction based on Markov tree in global searching

According to (22) and (24), the energy of the whole water network (17) can be obtained. The optimum solution of labels $\mathit{K}$ can be solved by searching the minimum of the objective energy function in space $\mathbf{W}$ using the Monte Carlo method or the simulated annealing method [30].

If we add the rate parameter $\xi$ in a Gibbs distribution, i.e., $p\left(\mathit{K}\right)\propto \exp \left(-\xi E\left(\mathit{K}\right)\right)$, the steps of the simulated annealing algorithm (Algorithm 1) are as follows.   

Algorithm 1: Water network construction based on Markov tree in global searching.

Step 1: Initialization. Initialize the label vector of water network branches to ${\mathit{K}}^{\left(0\right)}$ randomly. Set the limited angle of branches to ${\theta }_{ref}$ and the limited distance to ${\delta }_{ref}$. Set the maximum number of the children nodes of tree to N, the initial parameter $\xi$, the learning rate $lr$, the step size of parameter update $N_{up}$, and the maximum number of iterations $N_{max}$. Step 2: Build the network graph model. For branches with label 1 in ${\mathit{K}}^{\left(i\right)}$, the angle and distance between each two branches are calculated first. Take the branch with the maximum area as the root node. The successor nodes of the N-tree are then searched in distance limitation ${\delta }_{ref}$ and direction limitation ${\theta }_{ref}$. The children nodes are set to root nodes, and the successor nodes are iteratively searched until the leaf node is reached. Mark all the nodes in the N-tree that have been searched. Following the same steps, the unmarked node with the maximum area is set as the new root node to search a new tree until all the branches with label 1 are marked. After all the sub-trees ${\psi }^{\left(i\right)}$ have been searched, the total energy of the network $E^{\left(i\right)}\left(\mathit{K}\right)$ is calculated using ${\psi }^{\left(i\right)}$. Step 3: Update the state variables. Flip the label $k_j$ of branch $d_j$ randomly to update ${\mathit{K}}^{\left(i\right)}$ to ${\mathit{K}}^{\left(i+1\right)}$. Build the network graph again according to step 2 and update the energy of network to $E^{\left(i+1\right)}\left({\mathit{K}}^{\prime }\right)$. Accept the change of label $k_j$ according the following probability: $$a_{\xi }=min\left(1,\exp \left(-\xi \left(E^{\left(i+1\right)}\left({\mathit{K}}^{\prime }\right)-E^{\left(i\right)}\left(\mathit{K}\right)\right)\right)\right).$$ (25) Step 4: Update parameter $\xi$. Repeat step 3 until the maximum number of iterations is reached, where $\xi$ is updated to $\xi \ast lr$ every $N_{up}$ iterations.

(4)

Network construction based on Markov tree in local searching

Because the water network is complex, as the number of branches is large and the sizes of different branches are varied, global searching is difficult to converge to the optimal value. For SAR images, the area of some tributaries is very small, and the contours are uneven due to speckle noise. It is easy to ignore the branches with weak energy if the global searching algorithm is used. However, an important feature of the N-tree is that the global optimum must be the local optimum. According to (21) and (24), the energy of an N-tree can be decomposed into the energy of the root node, adding all the energies of the sub-tree, taking each child node as root. Thus, the global optimum can be dynamically searched locally. If $\mathit{H}=\left\{h_x\right\}$ denotes the extracted bifurcations of the main trunk, the water network construction algorithm based on local iterative searching (Algorithm 2) can be described as follows.

Algorithm 2: Water network construction based on Markov tree in local searching.

Step 1: Extract the main trunk and bifurcations. The branch with the maximum area $d_m$ from the branch set $\mathit{D}=\left\{d_j\right\}\left(j\in \left\{1,\dots ,N^d\right\}\right)$ is selected as the main water. Extract the main trunk and bifurcations $\mathit{H}=\left\{h_x\right\}$ using the scanning method in a scanning radius ${\delta }_{rad}$. Step 2: Build the network graph of a tributary. Take $d_m$ as the root node and start with bifurcation $h_x$ in $\mathit{H}$ to build a tree. Calculate the angle and distance between $h_x$ and others branches $d_j\left(j\ne m\right)$. Search the children nodes of $h_x$ in the tree under the limitation of ${\theta }_{ref}$ and ${\delta }_{ref}$. Then, take each child node as the root node to search its children nodes until the leaf node is reached to get the current sub-tree ${\psi }_x$. Step 3: Search the Markov tree using simulated annealing. Take ${\psi }_x$ as a water network and search the optimal tree using simulated annealing similar to Algorithm 1. Step 4: Iteration. Repeat steps 2 and 3 until all elements in $\mathit{H}$ complete the tree construction.

3.4. Bridge Detection and Bridge Body Extraction

When the final water network is obtained, all the adjacent branch pairs can be easily extracted by determining all edges of the network using the depth traversal algorithm of the tree because each edge denotes an adjacent branch pair. For an extracted branch pair, if the ROI of the bridge is directly extracted by merging the points on the two contours by distance, some boundary points that are not on the bridge will be wrongly classified to the ROI of bridges when the branch is narrow enough. As shown in Figure 6a, two narrow branches are adjacent to each other, where the region in gray is the land. The region marked by the black rectangle box is the ROI extracted by distance merging. We find that part of the boundary on the land is wrongly extracted.

Because the angle between a bridge and the surrounding land is nearly 90 degrees, the junction between the bridge and the water contour is distinct. To avoid the wrong extraction, the ROI of a bridge is extracted by merging the distinct points extracted from the water contour.

The feature points of the water contour are extracted using the Douglas–Peucker (DP) splitting and merging algorithm [31]. The DP algorithm recursively splits a curve to two segments at the middle point of the curve until the distance between each point on the curve and the straight line linking two adjacent feature points is less than a distance tolerance ${\delta }_{tol}$. Given two adjacent branches shown in Figure 6b, if the set of feature points of branch 1 is ${\mathbf{P}}_1$ and branch 2 is ${\mathbf{P}}_2$, the feature point pairs in close distance can be extracted under the threshold of the reference distance ${\delta }_{ref}$. The feature points on the boundary can then be obtained by the spatial coordinate relationship. If the feature points of branch 1 are $\left\{o_{11},o_{12}\right\}$ and branch 2 are $\left\{o_{21},o_{22}\right\}$, the pixels on the land inside the quadrilateral surrounded by points $\left\{o_{11},o_{12},o_{21},o_{22}\right\}$ are recognized as the bridge body. Sometimes only one feature point can be detected in branch 1 or 2, so the bridge body can be recognized as follows.

In the case of $o_{11}=o_{12}$ and $o_{21}=o_{22}$, the segment linking $o_{11}$ to $o_{21}$ is recognized as the bridge body.

In the case of $o_{11}\ne o_{12}$ and $o_{21}=o_{22}$, given a land pixel s in the quadrilateral, if $\left(o_{11}{o}_{12}\times o_{11}s\right)·\left(o_{11}s\times o_{11}{o}_{21}\right)>0$, $\left(o_{12}{o}_{21}\times o_{12}s\right)·\left(o_{12}s\times o_{12}{o}_{11}\right)>0$ and $\left(o_{21}{o}_{11}\times o_{21}s\right)·\left(o_{21}s\times o_{21}{o}_{12}\right)>0$, then s is inside the triangle $o_{11}{o}_{12}{o}_{21}$ and is recognized as a pixel on the bridge body, where “×” denotes the vector cross product, “·” denotes the vector dot product, and if $\mathbf{a}=\left(a_x,a_y\right)$, $\mathbf{b}=\left(b_x,b_y\right)$ then $\mathbf{a}\times \mathbf{b}=a_x{b}_y-a_y{b}_x$.

In the case of $o_{11}=o_{12}$ and $o_{21}\ne o_{22}$, the rule is the same as the case of $o_{11}\ne o_{12}$ and $o_{21}=o_{22}$.

In the case of $o_{11}\ne o_{12}$ and $o_{21}\ne o_{22}$, given a land pixel s in the quadrilateral, if $\left(o_{11}{o}_{12}\times o_{11}s\right)·\left(o_{11}s\times o_{11}{o}_{22}\right)>0$, $\left(o_{12}{o}_{21}\times o_{12}s\right)·\left(o_{12}s\times o_{12}{o}_{11}\right)>0$, $\left(o_{21}{o}_{22}\times o_{21}s\right)·\left(o_{21}s\times o_{21}{o}_{12}\right)>0$, and $\left(o_{22}{o}_{11}\times o_{22}s\right)·\left(o_{22}s\times o_{22}{o}_{21}\right)>0$, then s is inside the quadrilateral $o_{11}{o}_{12}{o}_{21}{o}_{22}$ and is recognized as a pixel on the bridge body.

4. Experimental Results and Analysis

Three quad-polarization SAR data acquired by RADARSAT-2 and TerraSAR-X sensors are used to test the proposed algorithm. The ground truth (GT) was drawn with reference to the Google earth map. The algorithm parameters were obtained by making statistics of true bridges and water branches first. Five experiments were then carried out to evaluate the performance of the proposed algorithm. In the first experiment, data over part of Singapore region was tested to show the detailed process flow of the proposed method. In the second experiment, the proposed method was tested for different sites. In the third experiment, the proposed method was compared with the method that constructed the water network by merging water branches by contour distance [32]. In the fourth experiment, the performance of the proposed method using quad-polarization SAR images was compared with that of three different single-polarization channel data. Finally, the bridge body recognition performance of the proposed method was compared with that of the traditional method based on spatial relationships [3]. In all the tests, the performance of the bridge detection was evaluated using detection rate and false rate. The performance of the bridge body recognition was evaluated using the intersection over union (IoU) index.

4.1. Data Description

Single look quad-polarization SAR data acquired by RADARSAT-2 sensor over the Singapore region in 2013 and the Lingshui, Hainan province, China region in 2014 and Multi-look quad-polarization SAR data acquired by TerraSAR-X sensor over the Singapore region in 2014 are used to test the proposed algorithm. When the proposed method is carried out, the coherent matrix (directly transformed by the Sinclair matrix for single look data) is calculated as the input data. The details of the data are listed in

Table 1. Experimental dataset details; UTC stands for the international standard time; AOI stands for the angle of incident; $m\times m$ stands for meter times meter.

-	Scene	Sensor	Size	Resolution (m × m)	UTC	AOI (°)
1	Singapore	RADARSAT-2	$6161\times 4256$	$4.73\times 4.80$	19 January 2013 11:31:08	47.3
2	Lingshui	RADARSAT-2	$5421\times 3500$	$4.73\times 5.47$	12 June 2014 10:49:49	37.2
3	Singapore	TerraSAR-X	$5500\times 2500$	$2.06\times 6.59$	10 March 2014 11:07:06	34.7

. The size of image is $6161\times 4256$ for the Singapore image from RADARSAT-2, $5937\times 3920$ for the Lingshui image, and $5500\times 2500$ for the Singapore image from TerraSAR-X. The corresponding resolutions are 4.73 m $\times 4.80$ m, 4.73 m $\times 5.47$ m, and 2.06 m $\times 6.59$ m, respectively. The detailed parameters of data are listed in Table 1. The Pauli pseudocolor images and geographic legends of the three data are shown in Figure 7, where “R”, “G”, and “B” denote the three components of the Pauli vector. To keep the balance of color contrast, each component is divided by its mean when the image is shown. The data of Singapore from RADARSAT-2 is shown in Figure 7a, the data of Lingshui is shown in Figure 7b, and the data of Singapore from TerraSAR-X is shown in Figure 7c. The color code of the image is shown in Figure 7d. The geographic legends are shown in Figure 7e–g. The topology of the water in Singapore is so complex that it forms a network from Figure 7a,c. The topology of the water in Lingshui is simple, but the tributary is so narrow that the bridges are difficult to detect from Figure 7b. If the main water branch is taken as the root of a tree, the structure of the whole water network is very similar to a tree.

4.2. Parameters Setting

According to Section 3, the parameters of the proposed method include four parts: level set segmentation, the main water scanning, the water network connection, and the bridge bodies recognition. With respect to the water and land segmentation using the level set method, the curve regularization parameter ${\lambda }_r$ is set to 0.2 and ${\gamma }_r$ is set to 0.1, as in reference [18]. With respect to the main water scanning, supposing that the input image resolution is $R_x\times R_y$, the equivalent resolution in x and y directions is $R_{=}=\sqrt{\left(R_x^2+R_y^2\right)/2}$. The maximum length of the bridge is $L_b$. The maximum width of the bridge is $W_b$. Then, the equivalent value of $L_b$ is $L_b^{\prime }=L_b/R_{=}$ The equivalent value of $W_b$ is $W_b^{\prime }=W_b/R_{=}$. With respect to the main water scanning, the scanning radius ${\delta }_{rad}$ is initially set to $W_b^{\prime }/2$. As a scanned bifurcation is also a long strip area, its length is usually far larger than the width, so that the ratio of length to width ${\rho }_{LW}$ is usually far greater than 2. To avoid missing the scanned bifurcations of some water tributaries, the ratio of length to width ${\rho }_{LW}$ is set to a small value, 3, according to experimental statistics. Thus, when selecting the broken adjacent branches of the main water, the reference distance ${\delta }_{ref}$ is set to $W_b^{\prime }$ and the reference area ${\Delta }_{ref}$ to $\left(L_b^{\prime }\times 4L_b^{\prime }\right)$.

With respect to the water network connection, when building the graph of the tree structure, the maximum number of children nodes N is set to 2 or 3 and the reference angle to 45${}^{\circ }$. To set the optimum parameters for the energy calculation of the tree structure, we first extracted 270 water branches by segmenting data 1 in Table 1, as shown in Figure 8. The histogram distribution of three geometric parameters, including length $l_L$, area, and the defined width consistency $E_o\left(w\right)$, of the 270 water branches were then computed through statistics, as shown in Figure 9, where the area is processed by quadratic square root. Statistics results show that the median of the length is 28.0, the area is 19.6, and the width consistency is 0.68. We observe that the histograms of length and area are concentrated, but the histogram of consistency is scattered from Figure 9. When calculating the energy of the tree structure, to keep a balance among the three parts of energies, the parameter of node energy ${\mu }_c$ is set to 10, ${\gamma }_c$ to 1, $\rho$ to 0.5, and ${\lambda }_c$ to 1 because the counted median of the length and the square root area is about 10 times the width consistency. To keep a balance between the edge energy related to distance and the edge energy related to direction, parameter $\epsilon$ is set to 50 and $\eta$ to 1 because the average distance between two adjacent branches is almost 50 times to the average sine of the angle. When calculating the total energy, ${\lambda }_{\psi }$ is set to 1. When the simulated annealing algorithm is carried out, the initial parameter $\xi$ is set to 0.05, the learning rate $lr$ is set to 1.1, the updating step parameter $N_{up}$ is set to 50, and the maximum number of iterations $N_{max}$ is set to 1000, according to experimental statistics.

With respect to bridge bodies extraction, the distance tolerance for the extraction of feature points ${\delta }_{tol}$ is set to $0.1\sqrt{L_b^{\prime 2}+W_b^{\prime 2}}$ to ensure the endpoint junction of the bridge is correctly extracted. The distance tolerance for the merging of feature points is also set to ${\delta }_{ref}$.

4.3. Example Results

An image covering part of the Singapore region obtained by a $4\times 4$ multi-look processing of data 1 is tested to show the algorithm flow of the proposed method. The size of the image is $1000\times 1000$. The Pauli pseudocolor image is shown in Figure 10a. From Figure 10a, we find that the coastal topography of the data is complex. Bridges are distributed in various tributaries with different shapes. There is a large sea-crossing bridge and some small bridges over narrow tributaries. It is difficult to segment and connect the branches on the narrow tributaries. The result of water and land segmentation using level set is shown in Figure 10b. The observation result shows that water and land are correctly separated, but some low scattering regions in the land are wrongly classified as water regions. The extracted main water is shown in Figure 10c. The branches other than the main water are shown in Figure 10d. We observe that there are some tributaries in the main water from Figure 10c. The shapes of tributaries are varied, and some of them are so narrow that the branches on the tributary cannot be connected from Figure 10d. The bifurcation scanning results of the main water is shown in Figure 10e, where the regions in gray marked in the red box are the scanned bifurcations. The observation result shows that the bifurcations of the main water are well scanned by the orthogonal direction scanning method. The final connection result of the water network is shown in Figure 10f. We find that all tributaries are correctly extracted and some small branches on the tributaries are also correctly connected. The results of bridge detection and bridge body extraction are shown in Figure 10g, where the detected bridge bodies are marked in red. The corresponding regions of the detected bridges in the original image are shown in Figure 10h, where the detected bridges are marked by yellow boxes and the bridge bodies are drawn in green. We find that most of the bridges are correctly detected and the corresponding bridge bodies are also accurately extracted.

4.4. Performance Evaluation of the Proposed Method under Different Sites

To validate the robustness of the proposed method, the method was tested for different sites using the three testing datasets listed in Table 1. The results of data 1, data 2, and data 3 are shown in Figure 11, Figure 12 and Figure 13, respectively. The ground truth of the bridges is shown in Figure 11a, Figure 12a and Figure 13a, where the bridges are marked by green boxes. Different bridges are located on different branches in Figure 11a, Figure 12a and Figure 13a. An enlarged view of part of the image near a narrow tributary is presented in Figure 11b, Figure 12b and Figure 13b. Some bridges are too small to be clearly observed from Figure 11b, Figure 12b and Figure 13b. The water branch extraction result is shown in Figure 11c, Figure 12c and Figure 13c. There are so many branches of different sizes that they are difficult to connect. The water branch connection result is shown in Figure 11d, Figure 12d and Figure 13d. We observe that the water network is well constructed, and even some small branches on the narrow tributaries are correctly connected in Figure 11d, Figure 12d and Figure 13d. The results of bridge detection and bridge body extraction are shown in Figure 11e, Figure 12e and Figure 13e, where the detected bridges are marked by yellow boxes and the bridge bodies are drawn in green. The detection results of the enlarged part of Figure 11e, Figure 12e and Figure 13e corresponding to Figure 11b, Figure 12b and Figure 13b are shown in Figure 11f, Figure 12f and Figure 13f. Comparing Figure 11b, Figure 12b and Figure 13b with Figure 11f, Figure 12f and Figure 13f, we find that almost all the small bridges over the narrow tributary are correctly detected.

According to the ground truth, the detection performance of the three datasets are listed in

Table 2. Detection performance of the proposed method under different sites, where Num denotes the number of targets, Pd denotes the detection rate, and Pf denotes the false alarm rate.

Data	Targets	Correct	Pd (%)	False Alarm	Pf (%)
1	16	13	81.3	1	7.14
2	5	5	100	0	0
3	11	10	90.9	1	9.09

. For data 1, 13 of the 16 bridges are correctly detected and only 1 false alarm target is detected. The detection rate is 81.3%, and the false alarm rate is 7.14%. For data 2, all 5 bridges are correctly detected, and the number of false alarm targets is 0. The detection rate is 100%, and the false alarm rate is 0. For data 3, 10 of the 11 bridges are correctly detected, and only 1 false alarm target is detected. The detection rate is 90.9%, and the false alarm rate is 9.09%.

4.5. Comparison of the Proposed Method with Branch Merging Method

In the third experiment, the performance of the proposed method was compared with the one that constructs a water network by merging water branches by contour distance [32] (branch merging method). The comparison results are shown in Figure 14. The comparison results of the three testing datasets are shown in the three lines of Figure 14, respectively. Figure 14a1–a3,b1–b3 are results of the comparison method. Figure 14c1–c3,d1–d3 are results of the proposed method. Figure 14a1–a3,c1–c3 are results of the water network construction, where the detected bridge bodies are marked in red. Figure 14b1–b3,d1–d3 are results of the bridge detection, where the bridges are marked by yellow boxes and the bridge bodies are drawn in green. Compared with Figure 14a1–a3,c1–c3 as well as Figure 14b1–b3,d1–d3, the result of the proposed method is significantly better. The detection performance of the comparison method is listed in

Table 3. Detection performance of the branch merging method.

Data	Targets	Correct	Pd (%)	False Alarm	Pf (%)
1	16	12	75	68	85
2	5	5	100	16	76.2
3	11	10	90.9	57	85

. For data 1, 12 bridges are correctly detected but with a higher number of false alarm targets, at 68. The detection rate is 75%, but the false alarm rate is 85%. For data 2, all bridges are also correctly detected by the comparison method but with a higher number of false alarm targets, at 16. The detection rate is 100%, but the false alarm rate is 76.2%. For data 3, 10 bridges are correctly detected but with a higher number of false alarm targets, at 57. The detection rate is 90.9%, but the false alarm rate is 85%. Compared with the proposed method listed in Table 2, the number of correctly detected targets of the proposed method is larger, while the number of false alarms is greatly reduced.

4.6. Comparison of Quad-Polarization with Single-Polarization Data

To evaluate the gains of polarization data, the performance of quad-polarization data is compared with that of three different single-polarization channels. Compared with quad-polarization data, all the algorithm steps and parameters of the single-polarization data are the same except replacing the coherent matrix with intensity for the computation of likelihood energy. The results of data 1 are shown in Figure 15, where the detected bridges marked in red are shown on the water network result. The result of quad-polarization data is shown in Figure 15a. The results of HH, HV, and VV channels are shown in Figure 15b–d, respectively. We observe that there exist some false alarm targets for the HH and VV channels and some missing targets for the HV channel. As the procedures of water network construction and bridge body recognition are almost the same for the two kinds of data, the cause for that is the difference of the water–land segmentation results. For HH and VV channels, some low intensity regions along the ocean are wrongly recognized as bridge. That is because the low intensity regions connecting different small water regions inside the land are wrongly segmented, and then the area and direction of the water branch are changed. For the HV channel, the cause for error is that some pixels on the bridge are wrongly segmented as water, causing the bridge regions to be broken, and the two adjacent branches are then directly connected. The detection performance is listed in

Table 4. Performance comparison of the proposed method using quad-polarization and three single-polarization channels data in Singapore, where Polar denotes quad-polarization data, and HH, HV, VV denote HH, HV, and VV channels, respectively.

Data	Targets	Correct	Pd (%)	False Alarm	Pf (%)
Polar	16	13	81.3	1	7.14
HH	16	10	62.5	9	47.3
HV	16	8	50.0	5	38.5
VV	16	4	25.0	18	81.8

. The detection rate of the quad-polarization data is $81.3\%$, which is higher than that of the three single-polarization data. The highest detection rate is $62.5\%$ in the HH channel, and the lowest detection rate is $25.0\%$ in the VV channel among the three single-polarization channels. The false alarm rate of the quad-polarization data is $7.14\%$, which is less than that of the single-polarization data. The lowest false alarm rate is $38.5\%$ in the HV channel among the three single-polarization channels.

4.7. Comparison of the Bridge Body Recognition between the Proposed Method and the Spatial Method

As the summary described in Section 1, most of the existing bridge detection methods extract bridge bodies by the spatial relationship of bridge with water based on the water segmentation result. However, some boundary points that are not on the bridge will be wrongly classified to the bridge when the water branch is narrow if those spatial methods are used, as shown in Figure 6. In the last experiment, we compared the bridge body recognition performance of the proposed method with that of the traditional method based on spatial relationship (the spatial method). Both methods obtain the water result from the proposed water network construction method. The spatial distance threshold is set to ${\delta }_{ref}$, as in Section 3.2, for the spatial method. The results are shown in Figure 16, where the recognized bridge bodies marked in red are all mapped on the water network result. The result of the spatial method is shown in Figure 16a, where the recognized bridge bodies are marked in a blue box, and the extracted contour pixels by spatial relationship are marked in green. The local region marked in cyan in Figure 16a is shown in Figure 16b. We observe that some pixels of the recognized bridge bodies are not on the bridge. The result of the proposed method is shown in Figure 16c, where the extracted feature points are marked in green circles. The result of the same local region as Figure 16b is shown in Figure 16d. We find that all the recognized bridge bodies are accurately located on the bridge. We compute the IoU between the detected box and the ground truth shown in Figure 11a. The mean IoU index of the proposed method is $99.5\%$, but the index of the spatial method is only $68.9\%$. The performance of the proposed method is much better than that of the spatial method.

5. Discussion

With respect to the offshore bridge detection problem in SAR images, the proposed method introduces a novel solution by constructing the water network using the probability graph model of a Markov tree. The experimental results of two RADARSAT-2 data and one TerraSAR-X datum show that the false alarm rates of the proposed method are less than 10%, and the detection accuracies are larger than 85%.

The difficulties of offshore bridge detection include the extraction and connection of water branches. By defining the internal and external energies of the segmentation curve and embedding the energy into a level set function, all connected branches are accurately extracted by the level set segmentation method. The segmentation results of the three experimental data showed that the extracted water branches were accurate and smooth. Because the number of water branches are large, the structure of the branches is complex, and there are many small and narrow branches, it is difficult to connect the branches merging directly by distance. By building the graph structure of the connection model using the tree structure, the branches are effectively organized because the tree structure is ordered from root to leaf so that each parent branch only has two or three child branches. The observations of the three experimental data showed that the structure of the water of each datum is similar to a tree. By defining the node energy and the edge energy of the tree and searching the extremum of the energy using the simulated annealing algorithm, the global optimum solution of the water branch connection is solved. Experimental results showed that the false alarm rates of the proposed method were greatly reduced compared to the method of merging branches contour by distance. Because the sizes of bridges are varied, there will some boundary points of the land wrongly recognized as the ROI of bridges. By extracting the feature points of the branches and merging them by distance, the bridge bodies are more accurately recognized. Experimental results showed that the mean IoU index of the proposed method was 99.5%, which is much larger than 68.9% of the spatial method.

Although the water connection algorithm is complex, the time-consuming step is mainly in the water and land segmentation. For the segmentation, if the size of the image is $n\times n$, the time complexity is $O\left(n^2\right)$. For the bifurcation scanning, if the scanning windows are designed before scanning, the time complexity of the scanning of one pixel is $O\left(1\right)$. Because the size of the image is $n\times n$, the time complexity of the whole scanning is $O\left(n^2\right)$. However, it is easy to scan in parallel, so the time complexity of the entire scanning is $O\left(1\right)$. For the graph construction, the time complexity of the tree construction algorithm is $O\left(n\right)$. For the Markov tree construction, if the number of nodes is n and the number of iterations is m, the complexity of the simulated annealing algorithm is $O\left(mn\right)$ because the time complexity of the tree construction algorithm is $O\left(n\right)$. Thus, the most time-consuming step is segmentation, and the total time complexity of the proposed algorithm is $O\left(n^2\right)$. In the test, the time of each step of the proposed algorithm was calculated. The testing platform was MATLAB v9.5, and the CPU was an Intel Xeon at 3.6 GHz with 16 GB RAM. For data of Singapore in the example test, the time of segmentation is 135.7 s, bifurcation scanning is 0.85 s, graph construction is 2.46 s, and Markov tree construction is 2.78 s.

6. Conclusions

We have proposed a bridge detection method for quad-polarization SAR images, which is based on a water network connection using the Markov tree probability graph model. By establishing the general segmentation model of the whole water network using Bayesian criterion, the bridge detection problem is converted to a two-layer hierarchical model, including the segmentation of land and water and the connection of different water branches. According to the similarity between the structure of water network and a tree, the connection of water branches is converted to the construction of tree by taking each water branch as a node. By defining the energy of the node and edge of the graph, the optimum solution of the tree is recursively searched using a tree-structure simulated annealing algorithm. Because a level set segmentation model and the Markov tree probability graph connection model are used, the extracted water network is accurate and robust to speckle noise and scattering interference. Finally, to avoid wrong results caused by merging contour by distance, bridge bodies are extracted by merging the feature points of the contours of each branch pair of interest. Evaluations of the proposed method as well as comparisons with the bridge detection method of merging water branches by contour distance using three quad-polarization data acquired by RADARSAT-2 and TerraSAR-X show the potential of the proposed method. The detection rate of the proposed method is higher than 85%, and the false alarm rate is lower than 10%. Compared with the branch merging method, the false alarm rate of the proposed method is far less, whereas the detection rate is slightly higher. Compared with the spatial method, the IoU index between the recognized bridge body and ground truth of the proposed method is also much improved. Compared to the three different single-polarization channels of the quad-polarization data, the performance of quad-polarization data is much better because of accurate segmentation.

From the analysis and verification of the experimental results, we find that the proposed method still cannot detect bridges over the tributaries interfered with by strong scatterers. In future work, polarimetric scattering parameters that can reduce the strong scattering interference components in the water area will be extracted to segment water and land. The proportion of error segmented water regions will be effectively decreased by those polarimetric parameters. At the same time, before the water network is connected by the proposed method, we will take a path tracking algorithm, such as particle filter, to track the tributary broken by large distances.