Harbor Aquaculture Area Extraction Aided with an Integration-Enhanced Gradient Descent Algorithm

Abstract

In this study, the harbor aquaculture area tested is Zhanjiang coast, and for the remote sensing data, we use images from the GaoFen-1 satellite. In order to achieve a superior extraction performance, we propose the use of an integration-enhanced gradient descent (IEGD) algorithm. The key idea of this algorithm is to add an integration gradient term on the basis of the gradient descent (GD) algorithm to obtain high-precision extraction of the harbor aquaculture area. To evaluate the extraction performance of the proposed IEGD algorithm, comparative experiments were performed using three supervised classification methods: the neural network method, the support vector machine method, and the maximum likelihood method. From the results extracted, we found that the overall accuracy and F-score of the proposed IEGD algorithm for the overall performance were $0.9538$ and $0.9541$, meaning that the IEGD algorithm outperformed the three comparison algorithms. Both the visualized and quantitative results demonstrate the high precision of the proposed IEGD algorithm aided with the CEM scheme for the harbor aquaculture area extraction. These results confirm the effectiveness and practicality of the proposed IEGD algorithm in harbor aquaculture area extraction from GF-1 satellite data. Added to that, the proposed IEGD algorithm can improve the extraction accuracy of large-scale images and be employed for the extraction of various aquaculture areas. Given that the IEGD algorithm is a type of supervised classification algorithm, it relies heavily on the spectral feature information of the aquaculture object. For this reason, if the spectral feature information of the region of interest is not selected properly, the extraction performance of the overall aquaculture area will be extremely reduced.

1. Introduction

Given that the scientific planning of harbor aquaculture areas is an effective and sustainable way to develop fishery resources, the scientific determination of the area’s spatial distribution and the real-time monitoring of it are necessary and of critical importance. Moreover, the use of remote sensing techniques allows one to observe objects on the surface of the water through satellite images and overcome the shortcomings of traditional monitoring technology. These techniques have become an important means of dynamic monitoring in harbor aquaculture areas [1]. In recent years, some scholars have achieved excellent extraction results and application prospects in harbor aquaculture area extraction experiments, greatly promoting the development of remote sensing technology [2,3]. The state-of-the-art extraction methods used for harbor aquaculture areas can be divided into three main categories: visual interpretation methods, object-oriented methods, pixel-oriented methods, and neural network methods.

The visual interpretation method is the most commonly used extraction method [4,5]. For example, Zeng et al. used medium-resolution multispectral images for aquaculture pond extraction [6], with the results indicating that the visual interpretation method makes extraction of the target site more difficult. Furthermore, Luo et al. presented a dynamic monitoring method [7] that has important implications for future upgrades to improve aquaculture and other issues. By using Landsat images and through visual interpretation, Duan et al. identified the dynamic nature of the growth of aquaculture pond areas in coastal areas of China [8].

In general, the area provided to aquaculture ponds in most provinces initially increases rapidly and then stabilizes or begins to decline. Therefore, monitoring the extent of aquaculture areas is a matter of great urgency. Zhang et al. used GF-2 satellite images to achieve a superior extraction of aquaculture areas in turbid waters, providing technical support to the aquaculture industry [9]. However, the existence of visual interpretation relies on the visual interpreter’s own experience, which is not conducive to the monitoring needs of aquaculture areas because of the large and time-consuming workload involved.

Next, the object-oriented extraction method comprehensively considers the space, spectrum, texture, and shape features of classified objects in remote sensing images [10,11,12]. To demonstrate the advantages of this method, Peng et al. used the CBR (Case-Based Reasoning) method [13] to successfully accomplish the segmentation of multi-scale images to complete classification and extraction. The process of object-oriented pond culture information extraction consists of image preprocessing, edge extraction, image segmentation, feature analysis, and extraction, etc. [14]. Liu et al. [14] and Xu et al. [15] further optimized the object-oriented extraction technique to improve the classification accuracy [15] while quantitatively describing the change in sea area usage in the study area.

In addition, Wei integrated visual interpretation with the object-oriented approach to optimize the classification and extraction accuracy, thus, optimizing the object-oriented extraction [16]. However, the object-oriented approach also possesses non-negligible drawbacks, such as different types of features presenting the same spectral characteristics in a certain spectral interval, leading to their incorrect classification.

Additionally, the pixel-oriented extraction method can make good use of the spectral reflection characteristics of the aquaculture area and use the threshold to extract the aquaculture area automatically [17]. This was confirmed in the research of Hussain et al. [18] and Zheng el al. [19]. Affected by different water quality factors in the same culture area, the reflection characteristics of aquaculture areas differ, making the pixel-oriented extraction method more difficult to perform. In addition, the setting of the extraction threshold in the pixel-oriented extraction method also needs to be specified manually. When the aquaculture extraction range becomes large, it is difficult to select the threshold that will produce the best extraction results. All these factors make it difficult for the pixel-oriented extraction method to extract aquaculture areas independently and accurately.

With the development of deep learning, more neural networks are being applied for the monitoring of aquaculture areas. The reason for this is that neural networks have shown excellent capabilities in the areas of automated processing, self-learning, and parsing in complex environments. Therefore, in order to achieve the monitoring and extraction of aquaculture areas using neural networks, Cui et al. identified the shortcomings of visual interpretation and proposed a method based on the automatic extraction of floating raft culture areas using a fully convolutional neural network (CNN) model [20], which effectively improved the accuracy of floating raft culture area identification.

Moreover, Jiang et al. upgraded the CNN model to a 3D-CNN model [21], which is an extraction model with a high extraction accuracy and strong spatial migration ability in complex water backgrounds. Therefore, the presented model is suitable for extracting large-scale, multi-temporal offshore rafting areas from remote sensing images. Further, Liu et al. used GF-2 satellite images to construct a deep learning richer convolutional feature network model for the water and soil separation of areas [3], thus achieving effective extraction in areas with more sediment and waves.

On top of that, Cheng et al. provided a hybrid dilated convolution U-Net model [22], indicating that their method possesses excellent distinguishability. In addition, Fu et al. constructed a hierarchical cascaded homogeneous neural network for marine aquaculture extraction. This network showed a superior classification performance to that of other existing methods [23]. However, in the above neural network algorithms, it can be seen that they all require a certain number of training samples and take a long time, making them very unfavorable for large-scale image observation.

After the above explanations and analyses, this paper aims to address the matter of aquaculture extraction with a generalizable performance, simple operation, explicit implementation, and high accuracy. For this reason, an implementable and feasible constrained energy minimization (CEM) scheme for aquaculture extraction is considered. Moreover, combined with the actual application effect of remote sensing based on satellite images and the actual demands of the Zhanjiang aquaculture area, this research selected GF-1 satellite PMS images as the data source. It is important to bear in mind that gradient-type algorithms have been proven to possess an accessible implementation paradigm for use in various optimization problems.

Therefore, in order to further improve the extraction accuracy of the CEM scheme in aquaculture areas, this paper proposes the use of an integrated-enhanced gradient descent (IEGD) algorithm derived from the gradient descent (GD) algorithm. The innovations made by this work are as follows. The proposed IEGD algorithm is derived from a major innovation in control theory, which adds an integration error summation term to improve the computational precision. In other words, the integration error summation term can effectively eliminate the deviation between the computed and theoretical solutions. In general, the proposed IEGD algorithm demonstrates highly accurate extraction results in aquaculture areas. The main contributions and highlights are summarized as follows:

1.

Based on control theory, the proposed IEGD algorithm represents a major breakthrough of the traditional gradient descent algorithm for solving extraction problems and can be regarded as a novel algorithmic paradigm.

2.

The proposed IEGD algorithm possesses relatively excellent robustness.

3.

The proposed IEGD algorithm improves the shortcomings of the traditional CEM algorithm with an insufficient accuracy.

2. Research Principles and Methods

2.1. Technical Route

In view of the complex and diverse spectra of remote sensing images of objects on the ground in the research area, the method proposed in this paper is mainly divided into three stages—namely, spectral feature selection, an IEGD solving algorithm, and ground object extraction. First of all, this paper uses the feature index method and the Bhattacharyya Distance (BD) method to screen the sensitive wavebands to form the spectral feature set. Secondly, the CEM scheme aided by the proposed IEGD algorithm is used to enhance the target ground object and weaken the background spectral information. Finally, Otsu’s method is used to calculate the division threshold to achieve ground object extraction. Figure 1 shows the technical route map.

2.2. Spectral Feature Selection

Due to the different reflection characteristics of ground objects in different wave bands, the characteristic index method can be used to amplify the small brightness differences and highlight the spectral characteristics of ground objects. For situations where there is a high chlorophyll concentration and high suspended sediment content in Zhanjiang nearshore waters, the NDWI (Normalized Difference Water Index) feature index,

$$NDWI=\frac{{\mathrm{band}}_{green}-{\mathrm{band}}_{nired}}{{\mathrm{band}}_{green}+{\mathrm{band}}_{nired}},$$

is applied in this study [24] to reduce the influence of chlorophyll and suspended sediment on aquaculture area information.

After the above band expansion method has been applied, the new bands (NDWI) are obtained as well as the blue, green, red, and near-infrared bands of the GF-1 image for a total of 17 wavebands. These wavebands contain not only the characteristic information of the surface features in the culture area, but also the information of suspended sediment and chlorophyll concentration in non-culture waters. However, not all bands are sensitive to the target ground objects in the culture area.

In order to eliminate the band information that has a weak relation to the target ground object, we analyze the spectral characteristics of ground objects and use the Bhattacharyya Distance to screen the bands. Bhattacharyya Distance measurement determines the similarity of two discrete or continuous probability distributions by comparing the overlap between species in order to analyze the similarity between two statistical samples or populations [25].

$$BD=\frac{1}{4}\frac{{({\mu }_1-{\mu }_2)}^2}{{\sigma }_1^2+{\sigma }_2^2}+\frac{1}{2}\lg \left(\frac{{\sigma }_1^2+{\sigma }_2^2}{2{\sigma }_1{\sigma }_2}\right),$$

where $\mu$ represents the mean value of the gray-level eigenvalues of two different ground objects and $\sigma$ is the standard deviation between ground objects. The larger the BD value is, the greater the spectral difference between the ground objects is, and the easier it is to distinguish. In this study, all the above seven bands were screened using the Bhattacharyya Distance method, three bands with great differences from the ground objects were retained, and a band feature set for the ground object was formed using band synthesis technology.

Table 1. The BD ratio of each waveband.

BD	Value ($\mathsf{\mu }\mathbf{m}$)	BD	Value ($\mathsf{\mu }\mathbf{m}$)
${\mathrm{BD}}_{blue}$	1.133	${\mathrm{BD}}_{{\mathrm{band}}_{green}/{\mathrm{band}}_{red}}$	0.179
${\mathrm{BD}}_{green}$	0.632	${\mathrm{BD}}_{{\mathrm{band}}_{green}/{\mathrm{band}}_{nired}}$	5.850
${\mathrm{BD}}_{red}$	0.556	${\mathrm{BD}}_{{\mathrm{band}}_{red}/{\mathrm{band}}_{blue}}$	0.039
${\mathrm{BD}}_{nired}$	2.766	${\mathrm{BD}}_{{\mathrm{band}}_{red}/{\mathrm{band}}_{green}}$	0.173
${\mathrm{BD}}_{NDWI}$	0.289	${\mathrm{BD}}_{{\mathrm{band}}_{red}/{\mathrm{band}}_{nired}}$	2.696
${\mathrm{BD}}_{{\mathrm{band}}_{blue}/{\mathrm{band}}_{green}}$	0.800	${\mathrm{BD}}_{{\mathrm{band}}_{nired}/{\mathrm{band}}_{blue}}$	0.512
${\mathrm{BD}}_{{\mathrm{band}}_{blue}/{\mathrm{band}}_{red}}$	0.044	${\mathrm{BD}}_{{\mathrm{band}}_{nired}/{\mathrm{band}}_{green}}$	1.815
${\mathrm{BD}}_{{\mathrm{band}}_{blue}/{\mathrm{band}}_{nired}}$	1.928	${\mathrm{BD}}_{{\mathrm{band}}_{nired}/{\mathrm{band}}_{red}}$	1.143
${\mathrm{BD}}_{{\mathrm{band}}_{green}/{\mathrm{band}}_{blue}}$	0.419

shows the BD ratios of each waveband, where ${\mathrm{band}}_{blue},{\mathrm{band}}_{green},{\mathrm{band}}_{red}$, and ${\mathrm{band}}_{nired}$ represent the gray-level values of GF-1 satellites from blue, green, red, and near infra-red bands after atmospheric correction.

Combined with the spectral characteristics of the ground objects in the culture area, we can see that the spectral difference between the cage and the sea water in the raft culture area is the smallest. Thus, the BD value between the cage and the sea water is selected as the selection standard of the spectral feature set in the raft culture area. Therefore, the spectral feature set of the raft culture zone is mainly composed of ${\mathrm{band}}_{green}/{\mathrm{band}}_{nired}$, ${\mathrm{band}}_{red}/{\mathrm{band}}_{nired}$, and the near-infrared band.

2.3. IEGD Solving Algorithm

The Constrained Energy Minimization (CEM) scheme [26] is a classical hyperspectral image target detection algorithm that is widely used in the field of remote sensing image ground object extraction [27,28,29]. The traditional CEM scheme regards the input remote sensing image as a limited observation signal set and supposes that a hyperspectral image can be arranged as a matrix $S=[r_1,r_2,\cdots ,r_N]\in {\mathbb{R}}^{l\times N}$, where each column of S is a spectral vector and $r_i\in {\mathbb{R}}^{l\times 1}$ is a spectral vector representing the spectrum of targets of interest. In this, N is the number of pixels and l is the number of wavebands. The CEM scheme designs a finite impulse response (FIR) linear filter under the assumption that the prior information d of the target pixel spectrum is known and minimizes the average output “energy” of the image under the constraint of the following formula:

$${\omega }^{\mathrm{T}}d=\sum _{i=1}^l{\omega }_i{d}_i=1,$$

(1)

where $\omega =[{\omega }_1,{\omega }_2,\cdots ,{\omega }_N]\in {\mathbb{R}}^l$ is the filter coefficient, which is a l-dimensional vector. Suppose the output of the filter to input $r_i$ is $y_i$:

$$y_i={\omega }^{\mathrm{T}}{r}_i.$$

(2)

Then, the average output “energy” corresponding to $S=[r_1,r_2,\cdots ,r_N]\in {\mathbb{R}}^{l\times N}$ can be written as:

$$\begin{array}{cc}\hfill y& =\frac{1}{N}\sum _{i=1}^N{y}_i^2\hfill \\ \hfill & =\frac{1}{N}\sum _{i=1}^N{\omega }^{\mathrm{T}}{r}_i{r}_i^{\mathrm{T}}\omega \hfill \\ \hfill & ={\omega }^{\mathrm{T}}R\omega ,\hfill \end{array}$$

(3)

where $R=(1/N){\sum }_{i=1}^N{r}_i{r}_i^{\mathrm{T}}\in {\mathbb{R}}^{l\times l}$ is the autocorrelation matrix of the remote sensing image and $y=[y_1,y_2,\cdots ,y_N]\in \mathbb{R}$ represents the output “energy” of each pixel in the image after passing through the filer. Therefore, the CEM scheme can be described as a linearly constrained optimization mathematical model as follows:

$$\begin{array}{c}\hfill \underset{\omega \in {\mathbb{R}}^l}{min}\left\{{\omega }^{\mathrm{T}}R\omega \right\}\\ \hfill s.t.{\omega }^{\mathrm{T}}d=1.\end{array}$$

(4)

By using the Lagrange multiplier method, the linearly constrained optimization mathematical model (4) of the filter output can be transformed into an unconstrained optimization mathematical model as follows:

$$F\left(\omega \right)={\omega }^{\mathrm{T}}R\omega +\lambda (d^{\mathrm{T}}\omega -1),$$

(5)

where $\lambda$ is the Lagrangian multiplier. The mathematical model of unconstrained optimization (5) is transformed into a mathematical model of linear equality equation as follows:

$$Gs\left(t\right)=b,$$

(6)

where the autocorrelation coefficient matrix $G=[2R,d^{\mathrm{T}};d,0]\in {\mathbb{R}}^{(l+1)\times (l+1)}$, the coefficient vector $b={[0,1]}^{\mathrm{T}}\in {\mathbb{R}}^{l+1}$, $s\left(t\right)={[\omega \left(t\right),\lambda \left(t\right)]}^{\mathrm{T}}\in {\mathbb{R}}^{l+1}$ is the unsolved vector, $\omega \left(t\right)={[{\omega }_1\left(t\right),{\omega }_2\left(t\right),\cdots ,{\omega }_l\left(t\right)]}^{\mathrm{T}}$ is the l-dimensional vector composed of filter coefficients, and $\lambda \left(t\right)\in \mathbb{R}$ is the Lagrangian function multiplier. Then, the error function that defines Equation (6) is as follows:

$$e\left(t\right)=Gs\left(t\right)-b.$$

(7)

In order to better solve $s\left(t\right)$ for the CEM scheme (3), we must define a scalar-valued norm-based error function,

$$\epsilon =\frac{{\parallel e\left(t\right)\parallel }_F^2}{2},$$

(8)

where ${\parallel G\parallel }_F=\sqrt{tr\left(G^{\mathrm{T}}G\right)}$ represents the Frobenius norm of matrix G. Moreover, $tr(·)$ denotes the trace operator of a matrix. Additionally, a lemma that we must know [30] is:

$$\frac{\partial \epsilon }{\partial s\left(t\right)}=G^{\mathrm{T}}e\left(t\right).$$

Finally, the formula follows the form shown above:

$$\dot{s}\left(t\right)=-\alpha \frac{\partial \epsilon }{\partial s\left(t\right)}.$$

Combined with the algorithm mentioned above, the gradient descent algorithm can be obtained as:

$$\dot{s}\left(t\right)=-\alpha G^{\mathrm{T}}e\left(t\right).$$

(9)

According to automatic control theory, adding an integration error summation term to Formula (9) can increase the robustness of this formula:

$$\dot{s}\left(t\right)=-G^{\mathrm{T}}\left(\alpha e\left(t\right)+\beta {\int }_0^{t}e\left(\tau \right)\mathrm{d}\tau \right).$$

(10)

Finally, the discrete form of the integration-enhanced gradient descent (IEGD) algorithm is constructed as follows:

$$s_{k+1}=s_k-\alpha G^{\mathrm{T}}{e}_k-\beta G^{\mathrm{T}}\sum _0^k{e}_k,$$

(11)

where k denotes the iteration number. Through recursive calculation, the filter output coefficient obtained is inverted, and thus the aquaculture area target enhancement process is completed. Furthermore, the pseudo-algorithm of the proposed IEGD algorithm (10) is provided as Algorithm 1:

Algorithm 1: IEGD algorithm solving procedure.

1. Initial set $maxstep$, $\alpha$, $\beta$, $\omega$, $\lambda$, b, R, and d

2. Circular iteration

for ($k=$; k $<=$ $maxstep$; $k++$) do

calculate $G^{\mathrm{T}}=[2R,d;d^{\mathrm{T}},0]$   calculate $s=[\omega ;\lambda ]$   calculate $e=Gs-b$   calculate ${\sum }_0^k{e}_k=e_1+e_2+\cdots +e_k$   calculate $s_{k+1}=s_k-G^{\mathrm{T}}(\alpha e_k+\beta {\sum }_0^k{e}_k)$  end for

3. Output $s_{k+1}$

2.4. Theoretical Analyses

In this section, the convergence performance of the proposed model will be discussed at zero noise.

2.4.1. Conversion

In order to prove the convergence performance more conveniently, Equation (11) should be transformed. First, multiply matrix G on both sides:

$$G\dot{s}\left(t\right)=-\alpha G{G}^{\mathrm{T}}e\left(t\right)-\beta G{G}^{\mathrm{T}}{\int }_0^{t}e\left(\tau \right)\mathrm{d}\tau .$$

(12)

Based on the matrix theory, the real symmetric matrix $G{G}^{\mathrm{T}}$ can then be generated:

$$G{G}^{\mathrm{T}}\simeq \left[\begin{array}{cccc}{\gamma }_1& 0& \cdots & 0\\ 0& {\gamma }_2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\gamma }_{lN}\end{array}\right]\in {\mathbb{R}}^{l\times N},$$

where ${\gamma }_i(i=1,2,\cdots ,lN)$ is the eigenvalue of the positive definite matrix $G{G}^{\mathrm{T}}$ and the symbol ≃ denotes the similarity of the two matrices. In line with above derivation, Equation (11) can be transformed:

$$G\dot{s}\left(t\right)=-\alpha {\mathsf{\Lambda }}_{i}e\left(t\right)-\beta {\mathsf{\Lambda }}_i{\int }_0^{t}e\left(\tau \right)\mathrm{d}\tau ,$$

(13)

where $\mathsf{\Lambda }=G{G}^{\mathrm{T}}\in {\mathbb{R}}^{l\times N}$. Based on Equation (7), Equation (13) can be simplified to:

$$\dot{e}\left(t\right)=-{\epsilon }_{i}e\left(t\right)-{\zeta }_i{\int }_0^{t}e\left(\tau \right)\mathrm{d}\tau ,$$

(14)

where ${(·)}_i$ denotes the $i^{th}$ member of a vector, ${\epsilon }_i=\alpha {\mathsf{\Lambda }}_i$, and ${\zeta }_i=\beta {\mathsf{\Lambda }}_i$.

2.4.2. Stability Analyses

Theorem 1.

For a randomly generated initial state of the proposed IEGD algorithm (10), the residual error $e\left(t\right)$ converges to zero globally.

Proof of Theorem 1. 

Take the $i^{th}$ element arbitrarily from Equation (14) as a subsystem where:

$${\dot{e}}_i\left(t\right)=-{\epsilon }_i{e}_i\left(t\right)-{\zeta }_i{\int }_0^t{e}_i\left(\tau \right)\mathrm{d}\tau ,\forall i\in 1,\cdots ,lN.$$

(15)

Define a Lyapunov for Equation (15). This can be:

$${\dot{v}}_i\left(t\right)=-{\epsilon }_i{e}_i^2\left(t\right)-{\zeta }_i{\left({\int }_0^t{e}_i\left(\tau \right)\mathrm{d}\tau \right)}^2.$$

(16)

Then, begin the derivation of the constructed Lyapunov function (16):

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\dot{v}}_i\left(t\right)}{\mathrm{d}t}& =2e_i\left(t\right){\dot{e}}_i\left(t\right)+2{\zeta }_i{e}_i\left(t\right){\int }_0^t{e}_i\left(\tau \right)\mathrm{d}\tau \hfill \\ \hfill & =-2{\epsilon }_i{e}^2\left(t\right)\hfill \\ \hfill & ⩽0.\hfill \end{array}$$

(17)

Therefore, based on the Lyapunov theory, it can be concluded that $e\left(t\right)$ converges to zero globally. In summary, the theorem has been proven. □

2.4.3. Convergence Analyses

Theorem 2.

For the solution of the CEM model (3), the residual error of the proposed IEGD algorithm (10) converges exponentially.

Proof. 

Let $\theta \left(t\right)={\int }_0^{t}e\left(\tau \right)\mathrm{d}\tau$. Combining this with Equation (15), Equation (15) can be rewritten as:

$${\ddot{\theta }}_i\left(t\right)=-{\epsilon }_i{\dot{\theta }}_i\left(t\right)-{\zeta }_i{\theta }_i\left(t\right).$$

We arrange the above equation as:

$${\zeta }_i{\theta }_i\left(t\right)+{\epsilon }_i\frac{\mathrm{d}{\theta }_i\left(t\right)}{\mathrm{d}t}+\frac{{\mathrm{d}}^2{\theta }_i\left(t\right)}{\mathrm{d}{t}^2}=0.$$

(18)

According to the solution of the differential equation, two feature roots can be obtained as:

$$\begin{array}{cc}\hfill {\theta }_1& =\frac{-{\epsilon }_i+\sqrt{{\epsilon }_i^2-4\zeta }}{2},\hfill \\ \hfill {\theta }_2& =\frac{-{\epsilon }_i-\sqrt{{\epsilon }_i^2-4\zeta }}{2}.\hfill \end{array}$$

After that, in order to solve this differential Equation (18), initial values of ${\theta }_i\left(0\right)=0$ and ${\dot{e}}_i\left(0\right)=e_i\left(0\right)$ are given. There are three situations in Equation (18), which are discussed as follows:

• ${\theta }_1\ne {\theta }_2$, where ${\theta }_1,{\theta }_2\in \mathbb{R}$, i.e., ${\epsilon }_i^2>4{\zeta }_i$.
  For the solution of ${\epsilon }_i\left(t\right)$, $e_i\left(t\right)$ can be written as:

  $$\begin{array}{cc}\hfill {\epsilon }_i\left(t\right)& =\frac{e_i\left(0\right)(e^{{\theta }_{1}t}-e^{{\theta }_{2}t})}{\sqrt{{\epsilon }_i^2-4{\zeta }_i}},\hfill \\ \hfill e_i\left(t\right)& =\frac{e_i\left(0\right)({\theta }_1{e}^{{\theta }_{1}t}-{\theta }_2{e}^{{\theta }_{2}t})}{\sqrt{{\epsilon }_i^2-4{\zeta }_i}}.\hfill \end{array}$$
  Combining the above formulas, $e\left(t\right)$ can be obtained:

  $$e\left(t\right)=\frac{e\left(0\right)({\theta }_1{e}^{{\theta }_{1}t}-{\theta }_2{e}^{{\theta }_{2}t})}{\sqrt{{\epsilon }_i^2-4{\zeta }_i}}.$$
• ${\theta }_1\ne {\theta }_2$, where ${\theta }_1=\varphi +i\psi \in \mathbb{C},{\theta }_2=\varphi -i\psi \in \mathbb{C}$—i.e., ${\epsilon }_i^2<4{\zeta }_i$.
  For the solution of ${\epsilon }_i\left(t\right)$, $e_i\left(t\right)$ can be written as:

  $$\begin{array}{cc}\hfill {\epsilon }_i\left(t\right)& =e_i\left(0\right)sin\left(\psi t\right)e^{\frac{\varphi t}{\psi }},\hfill \\ \hfill e_i\left(t\right)& =e_i\left(0\right)e^{\varphi t\left(\frac{\varphi sin\left(\psi t\right)+cos\left(\psi t\right)}{\psi }\right)}.\hfill \end{array}$$
  Combining the above formulas, $e\left(t\right)$ can be obtained:

  $$e\left(t\right)=e\left(0\right)e^{\varphi t\left(\frac{\varphi sin\left(\psi t\right)+cos\left(\psi t\right)}{\psi }\right)}.$$
• ${\theta }_1={\theta }_2$, where ${\theta }_1,{\theta }_2\in \mathbb{R}$, i.e., ${\epsilon }_i^2=4{\zeta }_i$.
  The solution of the ${\epsilon }_i\left(t\right)$, $e_i\left(t\right)$ can be written as:

  $$\begin{array}{cc}\hfill {\epsilon }_i\left(t\right)& =e_i\left(0\right)e^{{\theta }_{1}t},\hfill \\ \hfill e_i\left(t\right)& =e_i\left(0\right)e^{{\theta }_{1}t}(1+{\theta }_1).\hfill \end{array}$$
  Combining the above formulas, $e\left(t\right)$ can be obtained

  $$e\left(t\right)=e\left(0\right)e^{{\theta }_{1}t}(1+{\theta }_1).$$

Based on the above three situations and according to the conclusion of Theorem 1, it can be concluded that the convergence speed converges exponentially to the conclusion from Theorem 1. In summary, the proof of this theorem is completed. □

2.5. Ground Object Extraction

Otsu’s method mainly divides the image into a target and background according to the gray-level difference of the image and determines the optimal threshold by obtaining the inter-class variance ${\sigma }_B^2$ between the ground objects. The algorithm is efficient and fast, and the execution efficiency is high [31,32]. The calculation formula is as shown in (19), where $P_1$ and $P_2$ are the probability that the image pixels are divided into the target pixels and the background pixels and $m_1$ and $m_2$ are the average gray-level values of the image pixels divided into target pixels and the background pixels.

$${\sigma }_B^2=P_1{P}_2{(m_1-m_2)}^2.$$

(19)

The enhanced remote sensing image is used to calculate the division threshold using Otsu’s method. In this process, because the overall gray level of some sea water pixels is similar to that of aquaculture objects, it is mistakenly categorized as the area of aquaculture features to improve the false detection rate of aquaculture areas. Therefore, aiming at the sea water retention phenomenon, this study used the characteristic that the gray-level value of sea water in a single band is smaller than that of cultured objects and carried out single-band threshold segmentation.

In view of the fact that the gray value of the sea water pixels in the raft culture area in the blue waveband is smaller than that in the cage, after many experiments, the blue waveband gray threshold $T_1>2000$ and the near-red waveband gray threshold $T_2>1300$ are used as the single-band threshold values in the raft culture area to segment the sea water and the raft. The preliminary extraction effect of the raft culture area is better, as the raft culture area separates most of the sea water from the raft, leaving only a small part of the high-turbidity water at the edge of the coast.

It is difficult to remove this water because the suspended sediment is too high. Finally, we use the object-oriented area attributes to remove the sea water. In this process, sea water with an area smaller than 43 pixels or larger than 700 pixels is removed to obtain the final extraction result.

3. Materials

3.1. Research Area

As shown in Figure 2, Zhanjiang is located on Leizhou Peninsula—the third largest peninsula in China, at the southernmost tip of the Chinese mainland. It is also an important aquaculture base for Guangdong Province and even for China [33]. The research object of this paper is mainly the raft culture area, which is divided into aquaculture rafts for oyster culture and cages for captive fish culture. Aquaculture rafts are mainly composed of bamboo or wood, and there is often a rope hanging down from the aquaculture raft to cultivate shellfish, such as oysters.

In the process of culturing, the rafts are closely arranged, the oysters cling onto the rope and are immersed in the water, and only the aquaculture raft floats on the water surface, showing a bright white rectangular object on the remote sensing image, as shown in Figure 2. The cage is usually composed of wooden boards and plastic buckets, and there are small fish ponds inside, which are shown as a rectangular grid similar to the color of the water in remote sensing images, as shown in Figure 2.

3.2. Data Sources

The GF-1 satellite is the first high-resolution and wide-width Earth observation satellite developed by China Aerospace Science and Technology Corp. The highest spatial resolution of the satellite is 2 meters (m), and the maximum imaging width is more than 800 kilometers (km), meaning that it can meet the needs of remote sensing data with various spatial and spectral resolutions [34,35]. The GF-1 PMS level 1 data are made up of five wavebands—namely, one panchromatic band and four multispectral wavebands [35]. The specific parameters are shown in Appendix A

Table A1. Parameters of GF-1 PMS multi-spectral satellite images.

Band Order	Value ($\mathsf{\mu }$m)	Spatial Resolution (m)
Pan 1-Panchromatic	0.450–0.900	2
Band 1-Blue	0.450–0.520	8
Band 2-Green	0.520–0.590	8
Band 3-Red	0.630–0.690	8
Band 4-NIR	0.770–0.890	8

.

The ENVI (Environment for Visualizing Images) 5.3 software is used to preprocess the remote sensing images, performing radiometric calibration, atmospheric correction, orthophoto correction, and so on. At the same time, in order to identify the culture area in the image more clearly, the GS (Gram–Schmidt) method is used to fuse the multi-spectral and panchromatic images, thus improving the spatial resolution of the multi-spectral image to 2 meters (m).

3.3. Accuracy Evaluation

In order to evaluate the accuracy of the detection results, the overall accuracy (OA), precision, recall, and F-score are used. The OA, precision, recall, and F-score are computed as follows [36]:

$$\left\{\begin{array}{c}\mathrm{OA}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}};\hfill \\ \mathrm{Precision}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}};\hfill \\ \mathrm{Recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}};\hfill \\ \mathrm{F}-\mathrm{score}=(1+{\beta }^2)·\frac{\mathrm{Precision}·\mathrm{Recall}}{{\beta }^2·\mathrm{Precision}+\mathrm{Recall}},\hfill \end{array}\right.$$

(20)

where TP denotes, if the true category of the sample is a positive case and the model predicts it to be a positive case; TN denotes, if the true category of the sample is a negative case and the model predicts it to be a negative case; FP denotes, if the true category of the sample is a negative case but the model predicts it to be a positive case; FN denotes, if the true category of the sample is a positive case but the model predicts it to be a negative case [37], and $\beta$ represents the weight. When $\beta =1$, this indicates that both precision and recall are important.

4. Result

In this section, we apply five extraction methods in the study region to compare their performance and judge the accuracy of their extraction rafts against the accuracy mentioned above. The five methods are the integration-enhanced gradient descent (IEGD) algorithm (10), the constrained energy minimization (CEM) scheme (3), neural network (NN) [38], support vector machine (SVM) [39], and maximum likelihood estimation (MLE) [40], which are presented in this paper. In addition, Figure 3 represents the original image of the study area.

4.1. Results of the Extraction Process

Not counting the time spent on image preprocessing in the ENVI software, the time taken to run the Algorithm 1 of the IEGD algorithm (10) in Matlab 2017A was $1.4863$ seconds. This time was less than two seconds for a large-scale image with $7575\times 8740$ pixels and 2-meter spatial resolution. It is worth noting that all the experiments were performed using MATLAB 2017A and ENVI 5.3 on a computer with Windows 10, the AMD Ryzen 5 3600 6-Core Processor @3.60 GHz and 32 GB RAM. In addition, certain results of the extraction process are shown below:

$$\left\{\begin{array}{c}d=[0.662;0.709;0.804];\hfill \\ \omega =[-9.1553;9.9557;0.0027];\hfill \\ \kappa =0.2863;\hfill \\ R=\left[\begin{array}{ccc}0.1717& 0.1701& 0.0419\\ 0.1701& 0.1703& 0.0419\\ 0.0419& 0.0419& 0.0224\end{array}\right],\hfill \end{array}\right.$$

where d is the prior information, $\omega$ is the filter coefficient, $\kappa$ is the optimal threshold, and R is the autocorrelation matrix of the remote sensing image.

4.2. Overall Performance

First, Figure 3 shows the original image of the whole raft breeding area, where the white color represents the aquaculture rafts. Figure 4 shows the graph of the results after extraction has been carried out using the IEGD algorithm (10).

Table 2. Comparison of the accuracy evaluation index of various methods in raft farming areas, where OA denotes the ratio of the predicted correct result samples to the total samples, precision denotes the ratio of the actual positive samples to the predicted positive samples, recall denotes the ratio of the predicted positive samples to the actual positive samples, and the F-score is a combined consideration of the precision and recall.

Accuracy Evaluation	Method
	IEGD	CEM	NN	SVM	MLE
OA	0.9538	0.9180	0.9070	0.9046	0.8290
Precision	0.9817	0.9839	0.8521	0.8461	0.7486
Recall	0.9279	0.8556	0.9847	0.9890	0.9940
F-score	0.9541	0.9153	0.9136	0.9120	0.8540

represents the comparison of the accuracy evaluation of various methods in the whole raft culture. From the overall view, the extraction effect is very satisfactory, basically extracting all the aquaculture rafts. This is demonstrated by the F-scores given in Table 2, which are above 0.96. Furthermore, it is necessary to explore the extraction accuracy of various methods locally.

4.3. Local Performance

In this subsection, we will apply five extraction methods in three different regions, of which this region is the region mentioned in the previous subsection. We will compare the performance of these five extraction methods and judge the accuracy of their extraction of rafts against the accuracy mentioned above.

Table 3. The accuracy evaluation index of the use of various methods in different regions, where OA denotes the ratio of the predicted correct result samples to the total samples, precision denotes the ratio of actual positive samples to predicted positive samples, recall denotes the ratio of the predicted positive samples to the actual positive samples, and the F-score is a combined consideration of the precision and recall.

Method	OA	Precision	Recall	F-Score
Region a
IEGD	0.9908	0.9907	0.9767	0.9837
CEM	0.9317	0.9919	0.7647	0.8636
NN	0.9169	0.7742	0.9966	0.8714
SVM	0.9146	0.7693	0.9966	0.8683
MLE	0.8466	0.6484	0.9990	0.7864
Region b
IEGD	0.9982	0.9951	0.9980	0.9965
CEM	0.9681	0.9978	0.8821	0.9364
NN	0.9445	0.8275	1.0000	0.9056
SVM	0.9402	0.8167	1.0000	0.8991
MLE	0.8000	0.5711	1.0000	0.7270
Region c
IEGD	0.9565	0.8908	0.9516	0.9202
CEM	0.9565	0.8908	0.9516	0.9202
NN	0.8201	0.5944	0.9987	0.7453
SVM	0.8183	0.9995	0.7537	0.8594
MLE	0.6394	0.4222	0.9993	0.5937

indicates the different accuracies of these five methods in different regions. In the following pictures of the comparison of different regional methods, there will be three different colors, where green stands for misclassified aquaculture areas, red stands for omitted aquaculture areas, and black represents the background.

4.3.1. Region a

In Region a, we can clearly observe that the proposed IEGD algorithm (10) possesses superior performance in the extraction of Region a. It can be clearly seen in Figure 5 that only the proposed IEGD algorithm (10) neither misses aquaculture rafts nor misidentifies other objects that are not aquaculture rafts as aquaculture rafts. In the accuracy for each method in each region in Table 3, it is clear that the accuracy of the proposed IEGD algorithm (10) is above 0.97, while the accuracies of the rest of the methods are worse to different degrees.

In terms of details, although the overall extraction accuracy of the IEGD algorithm (10) is the best, in terms of recall, MLE [40] performs better and exhibits no missed extraction. However, because of the defects of this method, it results in the lowest precision rate among the five methods, at 0.6484. In turn, it can be demonstrated in Table 3 that the precision and recall rates of the proposed IEGD algorithm (10) are 0.9907 and 0.9767, respectively, while the F-score value of 0.9837 demonstrates the excellence of this algorithm. Therefore, in Region a, the proposed IEGD algorithm (10) performs the best out of the five algorithms.

4.3.2. Region b

Figure 6 shows the image information of Region b. Due to there being less interference in the image information, the correct extraction accuracy of these five methods rebounded compared to that obtained in region a. The extraction performance of the proposed IEGD algorithm (10) is still excellent; it can perfectly extract the rafts in this region, and its accuracy is above 0.99. Compared with Region a, its general accuracy OA is improved except for MLE.

Compared with Region a, the range missed by the CEM scheme (3) was drastically reduced, while, for NN [38], the SVM [39] misjudgment was also significantly reduced. This is because the CEM scheme (3) still maintains a very high precision of 0.9978, while that of the proposed IEGD algorithm (10) differs by only 0.0059. However, in terms of recall, the CEM scheme (3) differs from the other methods by more than 0.1, which is a large error. Nevertheless, the proposed IEGD algorithm (10) achieved a recall of 0.9980, representing a difference of only 0.0020 from the remaining three methods. According to the overall F-score, the proposed IEGD algorithm (10) still performs well.

4.3.3. Region c

The accuracy of Region c decreased compared to that of Regions a and b. This is most likely due to the fact that the rafts in this region are too dense with cargo as well as with interference from some objects. As shown in Table 3, the accuracy of all methods decreased. Surprisingly, there is no difference between the accuracy of the proposed IEGD algorithm (10) and that of the CEM scheme (3).

However, in the raft extraction of region c, there is a missed extraction and a false extraction. The accuracy of SVM [39] surpasses that of the CEM scheme (3) at 0.9995, but its recall is the lowest among these five methods, at only 0.7537. In terms of recall, MLE [40] is the best without exception, achieving a recall of 0.9993. In terms of combined precision and recall, the proposed IEGD algorithm (10) and CEM scheme (3) perform the best in region c extraction.

5. Discussion

5.1. Robustness Discussion

1.

From the viewpoint of optimization theory, the proposed IEGD algorithm (10) and the other four compared methods are supervised learning classification algorithms that are considered as optimization models. Differing from the others, the proposed IEGD algorithm (10) improves the model computing accuracy and enhances the model robustness by adding an integration error summation term, which can be supported by the stability analysis and convergence analysis described in Section 2.4.2 and Section 2.4.3.

2.

Our proposed method should also be discussed from the viewpoint of experimental results. The three representative Regions a, b, and c all show the phenomena of “same object with different spectra” and “same spectrum with different objects”. These are undoubtedly the best testing examples for evaluating the robustness performance of the actual extraction effectiveness. According to the corresponding visual and quantitative results shown in Figure 4, Figure 5, Figure 6 and Figure 7 and Table 2 and Table 3, the robustness of the proposed IEGD algorithm (10) is the best.

5.2. Applicability Discussion

1.

The proposed IEGD algorithm (10) makes full use of the spectral, textural, and spatial geometric feature information of the Zhanjiang offshore aquaculture area and increases the feature dimension of aquatic image elements in the aquaculture area by constructing various feature indexes. Moreover, the proposed IEGD algorithm (10) enhances the target feature information in the aquaculture area, making the difference between the aquaculture objects and non-aquaculture objects more distinct. The proposed IEGD algorithm (10) can better expand the feature information between the target object and background object to improve the performance of supervised learning classification, allowing it to effectively overcome the phenomena of “same object with different spectra” and “same spectrum with different objects”.

This algorithm is relatively reliable in cases where there are rich spectral features and a high extraction accuracy. It can be used as a self-selected algorithm for extracting aquaculture areas from high-resolution remote sensing images. In cases where there are fewer spectral features, it can also be combined with other existing methods to achieve a satisfying extraction performance.

2.

Differing from most existing extraction methods that only focus on local small areas, the proposed IEGD algorithm (10) can be directly employed for large-scale remote sensing images and can achieve an overall extraction of full-frame images.

3.

The proposed IEGD algorithm (10) is an important breakthrough for supervised learning classification algorithms, which is attributed to the integrated processing of the Matlab and ENVI software. When processing large-scale remote sensing images, the use of this algorithm can guarantee the accuracy of local feature extraction and provide a fast extraction speed. This can be seen in the extraction results of the experimental process. Not counting the time spent on image preprocessing in the ENVI software, the time generally taken to run the Algorithm 1 of the proposed IEGD algorithm (10) in Matlab 2017A is $1.4863$ s.

This time is less than two seconds for a large-scale image with $7575\times 8740$ pixels and a 2-meter spatial resolution. In addition, the overall accuracy and F-score of the proposed IEGD algorithm (10) in terms of the overall performance are $0.9538$ and $0.9541$, meaning that it outperforms the other four comparison algorithms. This demonstrates the excellent extraction performance of the proposed algorithm in aquaculture areas.

5.3. Expansibility Discussion

1.

The proposed IEGD algorithm (10) belongs to the category of optimization methods. Specifically, it improves the traditional gradient descent method by introducing an integration error summation term to help its optimal solution process and obtain a higher-precision computational solution. In this way, the proposed algorithm can be extended to other remote sensing-like algorithms that are applicable to the gradient descent method to help improve their solution accuracy. It is well known that the gradient descent method is a widely used algorithm in application scenarios; thus, the proposed IEGD algorithm (10) possesses high expansibility, good implementability, and acceptable feasibility.

2.

The proposed IEGD algorithm (10) is currently used for single aquaculture objects (i.e., rafts), but it can be extended to multiple aquaculture object extraction tasks by adding the prior spectral information of multiple targets.

3.

The proposed IEGD algorithm (10) can not only be employed for GF-1 remote sensing images but also for other multi-source remote sensing images, especially hyperspectral remote sensing images.

4.

The proposed IEGD algorithm (10) can be used not only for extraction from aquaculture objects but also for other areas, such as mining areas, surface water on land, crops, etc.

5.4. Limitations

1.

The proposed IEGD algorithm (10) relies heavily on the spectral features of the aquaculture features as prior information and still lacks the ability to fully exploit and utilize the local image element dependencies. As a result, some preprocessing of the remote sensing images of the aquaculture features is needed in order to render the spectra features sufficiently reliable.

2.

Supervised learning classification requires the manual selection of regions of interest (ROI), which can be easily influenced by manual subjectivity.

3.

For supervised learning classification, our algorithm can only determine the ROI in defined regions that have been manually selected. It relies on human subjective selectivity and can easily miss some tiny regions, leading to a significant reduction in the extraction performance for the overall aquaculture area.

6. Conclusions

In this study, taking the offshore raft culture area of Zhanjiang City, Guangdong province, as the research object, we used domestic high-resolution GF-1 PMS images as the data source; studied the extraction rules for raft farming areas in high-resolution remote sensing images; and proposed an IEGD algorithm (10) that combines spectral features of remote sensing image features, the threshold method, etc. The neural network method, support vector machine method, and maximum likelihood estimation of traditional supervised classification methods were compared with the proposed IEGD algorithm (10) in this paper.

The results show that, even in culture areas where the phenomenon of “foreign bodies with the same spectrum” is obvious, this method can effectively overcome the interference of background ground objects and obtain high-precision raft culture area extraction results. Compared with supervised classification, the overall accuracy of this method is more than 0.98. At the same time, for the extraction of other ground objects (including different crops and other types of culture areas), it is expected that this method will also provide a better application effect.

In addition, through the research in this paper, the IEGD algorithm (10) can better help others to observe the distribution of rafts in the farming area, which makes the layout of the farming area reasonable. The limitation of this paper is that the IEGD algorithm (10) relies on the spectral information of the farming features as a priori information to be extracted, which happens to be a great drawback. Since the a priori information needs to be selected manually, subjective human judgment is extremely important. Applying this method to more areas and improving it will be subjects of our future work.