Improving Satellite Retrieval of Coastal Aquaculture Pond by Adding Water Quality Parameters

: Coastal aquaculture is an important supply of animal proteins for human consumption, which is expanding globally. Meanwhile, extensive aquaculture may increase nutrient loadings and environmental concerns along the coast. Accurate information on aquaculture pond location is essential for coastal management. Traditional methods use morphological parameters to characterize the geometry of surface waters to differentiate artiﬁcially constructed conventional aquaculture ponds from other water bodies. However, there are other water bodies with similar morphology (e.g., saltworks, rice ﬁelds, and small reservoirs) that are difﬁcult to distinguish from aquaculture ponds, causing a lot of omission/commissioning errors in areas with complex land-use types. Here, we develop an extraction method with shape and water quality parameters to map aquaculture ponds, including three steps: (1) Sharpen normalized difference water index to detect and binarize water pixels by the Otsu method; (2) Connect independent water pixels into water objects through the four-neighbor connectivity algorithm; and (3) Calculate the shape features and water quality features of water objects and input them into the classiﬁer for supervised classiﬁcation. We selected eight sites along the coast of China to evaluate the accuracy and generalization of our method in an environment with heterogeneous pond morphology and landscape. The results showed that six transfer characteristics including water quality characteristics improved the accuracy of distinguishing aquaculture ponds from salt pans, rice ﬁelds, and wetland parks, which typically had F1 scores > 85%. Our method signiﬁcantly improved extraction efﬁciency on average, especially when aquaculture ponds are mixed with other morphological similar water bodies. Our identiﬁed area was in agreement with statistics data of 12 coastal provinces in China. In addition, our approach can effectively improve water objects when high-resolution remote sensing images are unavailable. This work was applied to open-source remote sensing imagery and has the potential to extract long-term series and large-scale aquaculture ponds globally.


Introduction
Fish products, rich in nutrients, are indispensable foods for human beings. According to the Food and Agriculture Organization of the United Nations, food fish consumption per capita grew from 9.0 kg to 20.5 kg in the past 50 years [1,2]. Aquaculture is an important way to meet people's demand for fish products and has undoubtedly become the fastestgrowing food production sector in the world [3,4]. The coastal aquaculture industry has been one of the most dynamic growth sectors over the past 20 years for reasons such as addressing employment, developing the economy, and eradicating poverty [5,6]. Expansion of aquaculture ponds has been found along the coastline worldwide [7][8][9][10][11][12][13][14][15][16].
(i) test the usefulness of adding water quality characteristics in aquaculture pond extraction, (ii) develop a new method in aquaculture pond extraction with large area applicability, and (iii) test our method along the coast of China to reveal the expansion of aquaculture. These results have important implications for the extraction of aquaculture ponds and the management of aquatic ecosystems in coastal areas globally. All acronyms used in the paper are given in Table 1.

Study Area
Our study area was the 50 km buffer zone from the coast of China (Figure 1). The landscape in coastal areas of China shows a spatial variability from the north to the south [33]. In the north, coastal areas are low elevation plains (≤200 m) featured with regular and large aquaculture ponds (Figure 2f). In the south, the drainage system is well-developed in tropical and subtropical climates. The aquaculture ponds are irregularly distributed along the numerous waterways, and the average area of each pond varies dramatically according to local terrain ( Figure 2b). Therefore, we considered the shape, location, area of the aquaculture ponds, and the complexity of the water body types. To test if adding water quality parameters can efficiently improve the identification of aquaculture ponds, we chose 8 sites with diverse types of water bodies covering North and South China. For the convenience of representation, we use regions of interest (ROI) 1-8 in Table 2 to correspond to each site in order.
The test sites in Figure 2a-d are located in South China.
• Donghai Island, Zhanjiang City, Guangdong Province ( Figure 2d): extensive aquaculture pond system built on the tidal flats of small straits: irregular shape, large boundary curvatures, and the area of each pond vary dramatically. • Gulao Water Town, Heshan City, Guangdong Province ( Figure 2c): a typical intensive aquaculture system with a river tourism industry across the main river. These ponds naturally developed according to the meandering rivers. Furthermore, the water supply system and pathway are intertwined with aquaculture ponds for tourism. These characteristics increased difficulties in teasing out aquaculture ponds only. • Wukan Port, Lufeng City, Guangdong Province ( Figure 2b): a semi-intensive aquaculture system based on reclaimed bays: ponds are approximately rectangular, and have wide ranges in area variability. Reeds field around the pond is a factor that confuses the aquaculture ponds extraction other than ditches (Field survey photos: Figure S1). • Pubagang Town, Taizhou City, Zhejiang Province ( Figure 2a): a natural aquaculture system on the coastline. Pathways and channels are built to introduce seawater for aquaculture, and the water quality is similar to that of seawater. Some ponds have an irregular shape, which was constrained by the natural coastline. The test site in are managed uniformly. A modern water purification system is used to balance the nutrient content in the pond water (Field survey photos: Figure S2). • Donggang, Rizhao City, Shandong Province ( Figure 2f): an intensive aquaculture system that used to be saltworks decades ago. In this system, there are many narrow paths between aquaculture ponds, which are difficult to tease out. (Field survey photos: Figure S3). • Chengkou Saltworks, Binzhou City (Figure 2g), Shandong Province: an intensive aquaculture system with both fish ponds and saltworks. The geometry and area of the evaporation ponds and crystallization ponds are similar to those of aquaculture ponds. In addition, the ponds were next to river banks, with irregular shapes (Field survey photos: Figure S4). • Tianjin Hangu Saltworks ( Figure 2h): a coastal aquaculture system next to saltworks. Large area ponds are restricted by coastlines and saltworks. The shape is not regular rectangles, and the curvature of the boundary is high. The salt farm and the aquaculture pond are jointly operated here (Field survey photos: Figure S5).
Each site has its own difficulties, which confuse aquaculture ponds and other water bodies (  The water objects produced in each area are displayed in (a1-h1), and each color represents an object. The water object is marked as T when accurately depicting the shape of the pond, otherwise, it is marked as F. The main error factors are marked in red letters below for the traditional shape parameter classification method.   Morphological features and surrounding landscape of aquaculture ponds in eight test areas (a-h). The water objects produced in each area are displayed in (a1-h1), and each color represents an object. The water object is marked as T when accurately depicting the shape of the pond, otherwise, it is marked as F. The main error factors are marked in red letters below for the traditional shape parameter classification method.  The Sentinel-2 Multispectral Instrument (MSI) samples 13 spectral bands: visible and NIR at 10 m, red edge and SWIR at 20 m, and atmospheric bands at 60 m spatial resolution [34]. Its high-resolution multispectral products have been widely used in mapping water bodies and water quality inversion [35]. The Google Earth Engine platform [36] provided all-weather radar images and high-resolution optical images from Sentinel 2. The mapping of aquaculture ponds is generally performed once a year, as aquaculture ponds are generally filled with water throughout the year. We calculated the median per pixel for all Sentinel-2 MSI images over the March to November period each year. Since aquaculture ponds are mostly full of water during this period, the median can characterize the spectral reflectance of the water surface over a long-term series. The band and spatial resolution information of the Sentinel-2 MSI image are listed in Table 3. Moreover, the number of images used every year and the availability of the image is shown in Figures S6 and S7.

Validation Datasets
We used GEE's random point function to generate 1000 random points in the eight test sites and took the water objects corresponding to the random points as training samples. Aquaculture ponds are usually closed reservoirs artificially constructed or transformed from natural water bodies. We referred to the optical high-resolution images to determine the attributes of the water bodies corresponding to random points (aquaculture or nonaquaculture). To balance the number between the two types of samples (aquaculture or non-aquaculture), we used an oversampling method: the samples of the minority class were repeatedly sampled multiple times. This balances the numbers between the two classes of samples (aquaculture or non-aquaculture), since random forests are sensitive to large numbers of differences between sample classes. The categories of water objects in each dataset were verified by communicating with the local residents. In addition, some water objects were confirmed during field surveys in Wukan Port, Rudong county, Donggang, Chengkou, and Tianjian Hangu Saltworks during the winter and summer of 2021.

Water-Object Extraction
We used the object-based method that comprehensively considers the adaptive imagery binarization threshold so that the connectivity algorithm could combine water pixels into a homogeneous water object [21]. First, the Normalized Difference Water Index (NDWI; Equation (1)) [37] was used to distinguish water from non-water: where Band3 and Band8 represent the Green and Near-infrared bands of Sentinel-2 MSI imagery. Second, we applied a 3 × 3 high-pass filter (Equation (2)) to sharpen the NDWI image (Figure 3a), which could enhance the contrast between pathway pixels and pond pixels ( Figure 3b) [38]. Figure 3g,h suggested that the sharpening filter enabled a more prominent  [39]. The sharpening filter did not significantly affect the adaptive binarization threshold of the NDWI image. Under the same binarization threshold, the pathway pixels marked as "non-water" in the sharpened image were more distinct (Figure 3c,d). In the end, we applied the "connected components" four-neighbor connectivity function in GEE to the binarized image to obtain the water object (Figure 3e,f).
Remote Sens. 2022, 14, x FOR PEER REVIEW 8 of 21 more distinct (Figure 3c,d). In the end, we applied the "connected components" fourneighbor connectivity function in GEE to the binarized image to obtain the water object (Figure 3e,f). (e,f) are water objects, each color represents an object. The water object produced by NDWI imagery is partially missing. The water object produced by Sharpened NDWI imagery has suitable separation and coverage.
The water object produced by the sharpened NDWI images had a higher degree of separation and coverage when the pond did not mix with surrounding rivers (Figure 3f). The maximum connection range of the GEE communication function is 1024 pixels. If the mixed pixels due to pathways were incorrectly classified as water, a large portion of aquaculture ponds would be missing (Figure 3e). The high-quality water object increased the ability to project real water bodies, including morphological features and water quality parameters. Moreover, the MNDWI (Based on Sentinel-2 MSI) and Sentinel-1 VH bands were also tested to create water objects. Although these two indicators were sensitive to water and non-water land types as well, we chose the NDWI method according to its performance (see Support text 1 and Figure S8).

Classification Parameter Selection
First, we selected five shape parameters for the water object: Area, Perimeter, Compactness ratio C, Aspect Ratio, and Landscape Shape Index (LSI). The Area and Perimeter of the water object are the most basic morphological parameters. Compactness ratio (e,f) are water objects, each color represents an object. The water object produced by NDWI imagery is partially missing. The water object produced by Sharpened NDWI imagery has suitable separation and coverage.
The water object produced by the sharpened NDWI images had a higher degree of separation and coverage when the pond did not mix with surrounding rivers (Figure 3f). The maximum connection range of the GEE communication function is 1024 pixels. If the mixed pixels due to pathways were incorrectly classified as water, a large portion of aquaculture ponds would be missing (Figure 3e). The high-quality water object increased the ability to project real water bodies, including morphological features and water quality parameters. Moreover, the MNDWI (Based on Sentinel-2 MSI) and Sentinel-1 VH bands were also tested to create water objects. Although these two indicators were sensitive to water and non-water land types as well, we chose the NDWI method according to its performance (see Support text 1 and Figure S8).

Classification Parameter Selection
First, we selected five shape parameters for the water object: Area, Perimeter, Compactness ratio C, Aspect Ratio, and Landscape Shape Index (LSI). The Area and Perimeter of the water object are the most basic morphological parameters. Compactness ratio C is a dimensionless constant in the interval [0, 1], which can effectively distinguish artificial water bodies from natural water bodies. The larger value of Compactness ratios C indicates a more complicated shape [40]; the Aspect Ratio is calculated based on the circumscribed rectangle of the object [41], and the value of aquaculture ponds is generally greater than 1. Likewise, Landscape Shape Index (LSI) is a quantitative measure of landscape complexity or heterogeneity. In general, more complex landscapes have larger perimeters or edges for a given area and therefore a higher LSI [42].
We selected 5 common water quality parameters: Normalized Difference Chlorophyll Index (NDCI), Phycocyanin (2BDA), Maximum Chlorophyll Index (MCI), Colored Dissolved Organic Matter (CDOM), and Surface Algal Bloom Index (SABI). Chl and PC reflect the trophic status of a water body. PC contains a distinct spectral absorption feature centered on 620 nm, and Chl corresponds to 685 nm. To make the PC compatible with Sentinel-2 MSI, Johansen et al. (2019) [43] simplified the empirical equation for the inversion of phycocyanin proposed by Wynne (2008) [44] with the ratio of the red band and the VRE5 band. We chose two indicators that reflect the chlorophyll concentration. The first one was the normalized difference chlorophyll index proposed by Mishra et al. (2012) [45], which is often used to predict the chlorophyll concentration in turbid coastal waters. The other was the Maximum Chlorophyll Index (MCI) [46] based on the fluorescence height baseline algorithm. Here, we used 705 nm as the fluorescent line. As the concentration of colored dissolved organic matter increases, the color of water also changed from dark green to brownish brownish-yellower [47], which was an effective indicator to distinguish natural water bodies from aquaculture ponds. Kutser et al. (2005) [48] mapped the CDOM value through the power function of the ratio of ALI band 2 (565 nm) and band 3 (660 nm). Accordingly, we applied the green (560 nm) and red (665 nm) bands of Sentinel-2 MSI to the power function to predict the CDOM concentration. Similarly, we also included SABI, which can directly detect the spatial distribution of surface and near-surface algae [49]. All shape and water quality parameters are listed in Table 4.

Classification Model Settings and Validation
We used the classification and regression tree model (CART) in a random forest (RF) classifier to extract the aquaculture pond and return the optimized parameter combination [51]. As an ensemble classifier, random forest has the advantages of better accuracy, strong ability to resist the noise of training data, and strong resistance to overfitting. It has been widely used to map land-use types in remote sensing data [52,53]. For building an RF classifier, the Number of Trees (Ntree) and the number of variable divisions (Mtry) are the two most sensitive variables. An increase in the Ntree value does not carry the risk of overfitting, and 500 is usually an acceptable value [54]. Considering that GEE has a computational memory limit, we choose 400, which can achieve a similar effect to 500 and save computational memory. The Mtry parameter is usually set to the square root of the number of input variables [28].
We randomly divided the basic sample dataset of each site into a training dataset and a test dataset, of which the training dataset accounted for 70% and the test dataset accounted for 30%. Such random splits were repeated 10,000 times in the base dataset at each site and fed into a random forest classifier for cross-validation. The purpose of this step is to reveal the effect of adding water quality parameters on the performance and stability of the aquaculture pond classifier.
Firstly, we selected 5 morphological parameters as the control feature combination 1shape feature combination (SF), because shape parameters are the most commonly used classification features for classifying aquaculture ponds. Secondly, set three experimental feature combinations: experimental feature combination 1-water quality parameter feature combination (WQF); experimental feature combination 2-mixed feature recursive feature elimination (MFRFE); experimental feature combination 3-mixed feature non-recursive feature elimination (MFNRFE).
The experimental feature combination has the same number of parameters as the control feature combination. In addition, in each cross-validation, the random forest classifier can sort the importance of all morphological parameters (5 in total) and water quality parameters (5 in total) according to the GINI index, eliminating unimportant features to save computing memory and improve the classification accuracy. There are two main methods of feature elimination, feature recursive feature elimination and feature nonrecursive feature elimination. The details of these two methods are mentioned in Guyon et al. (2003) [55]. Experimental feature combinations 2 and 3 are obtained by different elimination strategies.
The core of accuracy evaluation is the confusion matrix [56,57]. Some measures of the consistency between evaluation model predictions and actual surveys are proposed by the confusion matrix [58], user's accuracy (UA), producer's accuracy (PA): Producer s accuracy = N ij N : N i: represents the amount of sample data of class i. N :j represents the amount of class j in the prediction result of the model. PA represents the probability that a pixel is correctly classified, and reflects the reliability of a certain type of prediction. UA represents the probability that the classified pixels represent the class, reflecting the correct attribution of a certain class of sample data. These two measures evaluate whether the model is biased to correctly classify the majority class, resulting in a low accuracy rate for the minority class [56]. Helldén (1980) [59] developed a harmonic average of producer and user accuracy, called the average accuracy index. We assume that PA is equally important to UA, so it is defined as the F1-score as Equation (5): The F1-score is a metric that comprehensively considers the contribution of each category to the accuracy rate [60,61]. It is more convincing when applied to the feature elimination strategy than the overall accuracy rate. We chose the F1-score as the accuracy assessment metric.

Quantitative Assessment of Accuracy
The error matrix provided an effective basis for quantitative interpretation of the agreement between the classified map and the ground truth. Our cross-validation returned an F1-score for the four test feature combinations (SF, WQF, MFNRFE, and MFRFE; MFRFE, and MFNRFE are collectively referred to as MF). The Number of best performance (NBP), Median F1-score (F1M), and Standard Deviation (SD) of the test feature combinations in all cross-validation are used to illustrate performance and stability. The results are visualized in Figure 4 and Table 5.
The F1-score is a metric that comprehensively considers the contribution of each category to the accuracy rate [60,61]. It is more convincing when applied to the feature elimination strategy than the overall accuracy rate. We chose the F1-score as the accuracy assessment metric.

Quantitative Assessment of Accuracy
The error matrix provided an effective basis for quantitative interpretation of the agreement between the classified map and the ground truth. Our cross-validation returned an F1-score for the four test feature combinations (SF, WQF, MFNRFE, and MFRFE; MFRFE, and MFNRFE are collectively referred to as MF). The Number of best performance (NBP), Median F1-score (F1M), and Standard Deviation (SD) of the test feature combinations in all cross-validation are used to illustrate performance and stability. The results are visualized in Figure 4 and Table 5.   In all the test sites, MFRFE always received the highest F1M (average 0.94 is 0.1 higher than WQF and 0.3 higher than SF), the smallest SD (average 2.31 × 10 −2 is 0.12 × 10 −2 lower than WQF and 0.48 × 10 −2 lower than SF), and the most NBP (average 5.28 × 10 3 is 1.08 × 10 3 more than WQF and 3.95 × 10 3 more than SF). When there was correlation between optimized parameters, RFE could usually return more stable feature set [62]. This result suggested that the mixed feature combination (MFNRFE or MFRFE) was more reliable than SF and WQF.
Excluding the ROI8 area, SF has the lowest F1M (average 0.91). The lowest NBP (average 1.33 × 10 3 ) and the largest SD (average 2.79 × 10 −2 ) indicated that there were limitations and instability in extracting ponds based only on shape parameters. The more chaotic the morphological features of the water object at the site (ROI2, ROI3; Figure 4), the worse the SF classification effect.
In ROI1, the model trained by WQF obtained the smallest SD, and the largest NBP; in ROI2-6, the performance of WQF was close to that of mixed feature combination (MFNRFE or MFRFE), sometimes better than MFNRFE (ROI3). At ROI7-8 with saltworks, WQF's performance was the worst among the four test feature combinations. In the overall quantitative evaluation, the performance of WQF was unstable.

Accuracy Assessment by Visual Interpretation
Comparing the classification results under the same spatial scale with the physical reference was an effective approach to evaluate the classification accuracy. To evaluate the classification results, we visually interpreted the high-resolution images of Google Earth, as well as field survey photos.

Water Quality Feature Combination
It might problematic for WQF to handle mixed object classification. For example, a complex interaction between ponds and other water bodies could cause misclassification of the WQF approach. In Figure 5(a2), the reclamation ponds were identified as bays by mistake. In Figure 6(b2), diversion ditches were identified as aquaculture ponds by mistake. Although WQF could distinguish crystallization ponds, evaporation ponds with aquaculture ponds caused mistakes (ROI7, Figure 5

Shape Feature Combination
SF's ability to classify mixed water objects was limited due to the diverse shape aquaculture ponds (Figure 6(b1,e1,f1)). In the saltworks, the water bodies shared a si lar shape between the aquaculture ponds, crystallization ponds, and evaporation pon It was inevitable to make identification mistakes (Figures 5(b1) and 6(d1)). In South C na, the shapes of aquaculture ponds were irregular, which could be easily confused w other water bodies. Reed field (or rice field) was another factor that limited the classif tion accuracy (Figure 5(a1)). In some cases, the water object could depict the dive shapes of the ponds, in the ROI1. However, SF neither accurately identified the exp sion ponds nor the intensive ponds.

Shape Feature Combination
SF's ability to classify mixed water objects was limited due to the diverse shapes of aquaculture ponds (Figure 6(b1,e1,f1)). In the saltworks, the water bodies shared a similar shape between the aquaculture ponds, crystallization ponds, and evaporation ponds. It was inevitable to make identification mistakes (Figures 5(b1) and 6(d1)). In South China, the shapes of aquaculture ponds were irregular, which could be easily confused with other water bodies. Reed field (or rice field) was another factor that limited the classification accuracy ( Figure 5(a1)). In some cases, the water object could depict the diverse shapes of the ponds, in the ROI1. However, SF neither accurately identified the expansion ponds nor the intensive ponds.
In most cases, MF could keep the advantages of both SF and WQF. In ROI1, MFNRFE did not have a large-scale classification error. MFRFE misclassified a section of a river and a pond (Figure 6(a3,a4)). Compared with SF and WQF, the area of water objects with incorrect MF classification was reduced. In ROI2, MF had the same ability as SF to correctly classify freshwater rivers and reservoirs (Figure 6(b3,b4)). In ROI3, MF kept the ability to correctly classify marine aquaculture ponds and reed fields near the sea ( Figure 5(a3,a4)). WQF could not correctly determine the class of the former, and SF could not correctly determine the class of the latter. In ROI4, MF had the ability to identify mariculture bays ( Figure 6(c3,c4)), while WQF could not. WQF had a superior ability on the classification of mixed water bodies. More importantly, MF inherited the feasibility of WQF for the classification of mixed water objects. In ROI5 with a large amount of sample data, both MF and WQF could correctly classify mixed water bodies (Figure 6(e3,e4)). However, due to the small sample size of ROI6, the classification effect of MF was poor (Figure 6(f3,f4)), which was similar to the results of SF. MF could avoid the confusion between aquaculture ponds, evaporation ponds, and crystallization ponds (Figure 6(d3,d4)), especially when the sample size was large ( Figure 5(b3,b4)).

Generalization and Transferability of MF
We used probability maps of aquaculture ponds in each coastal province to evaluate MF's generalization ability, i.e., the adaptability of the classification approach in fresh samples [63]. We classify and train the south model and the north model. The training samples of the south model are all from the stations in the south, and the training samples of the north model are all from the stations in the north. Cross-validation was used to quantitatively explore the generalization performance of MF. The F1M of MF was 0.947 in the southern model and 0.908 in the northern model ( Figure 4). Then, we optimized 43 different MFs in the four southern location models and the southern model and applied them to the southern total sample dataset to retrain the aquaculture pond classification model (Table S1). The aquaculture pond classification model returns 43 classification results for each water object. We drew a probability map of aquaculture ponds in seven southern provinces based on the probability that each water object is classified as an "aquaculture pond". The same 79 different MFs were used to map the probability of aquaculture ponds in the five northern provinces (Table S2).

Validation with Yearbook of Statistics
We used provincial-level data in the "China Fishery Statistical Yearbook" to validate the area of aquaculture ponds extracted by MF. Our study area was focused on a 50 km buffer zone along the coastal line of each province (Hong Kong and Macau were not listed). Since the inland aquaculture ponds (distance to the coast greater than 50 km) in Guangxi and Jiangsu have wide distribution, it was reasonable that the yearbook statistic data were much higher than the area extracted by our algorithm (Figure 7). There are fewer inland aquaculture ponds in Tianjin and Hebei. The data in the yearbook were within the cumulative probability area of 50-100%. For the remaining provinces, the yearbook data were within the cumulative probability area of 25-100% (Figure 7). The area data extracted by our algorithm and the yearbook data are shown in Table S3.

Advantages in Adding and Water Quality Parameters
Our extraction approach suggested that the combination of shape and water quality parameters is efficient to identify large-scale aquaculture ponds. Taking 2017 as an ex-

Advantages in Adding and Water Quality Parameters
Our extraction approach suggested that the combination of shape and water quality parameters is efficient to identify large-scale aquaculture ponds. Taking 2017 as an example, we compare our results with other studies [27]. The provincial-level results showed that our extraction areas were roughly the same as these two studies (Figure 8). Due to the spatial resolution of the Sentinel-2 MSI, our results were improved to a higher resolution (i.e., 10 m). Different evaluation and verification procedures (e.g., field survey photos) have confirmed that the method of mapping national aquaculture pond probabilities and extracting aquaculture pond area through MF was effective. Because the salt drying farm and aquaculture ponds are similar in shape, the traditional shape parameter extraction of aquaculture ponds encountered troubles from the joint operation of the saltworks and aquaculture. However, different from the water used for aquaculture, the raw material of the salt farm is seawater without chemical agents. The water quality parameters can be used as a classification basis to monitor the difference between the two, regardless of the quality of the water object. Secondly, because the water quality characteristics of aquaculture ponds are similar, MF has the potential to correctly classify mixed water objects under the effect of water quality characteristics. Finally, the shadow noise of urban buildings affects the effect of NDWI in extracting water bodies [64], MF is not affected by urban shadow noise. Then MF still shows instability in some aspects. (1) Although the classification of mixed objects has been improved, the probability of correct classification is still within a large uncertainty range (25-100%). Some natural lakes are also classified as aquaculture ponds because of their probability of being classified as aquaculture ponds, making the overall pond area too large; (2) The sample data set is also an important factor affecting the MF's classification of aquaculture ponds. When the surrounding landscape is complicated, the error is particularly obvious. Therefore, it may be a better choice to allocate the number of different types of sample points based on the area ratio [65].
Our results also demonstrated advantages in three major deltas of China. In the Yellow River Delta, where saltworks are distributed, MF can more accurately realize the classification of aquaculture ponds and saltworks. The probability that the ponds in the saltworks are classified as aquaculture ponds is less than 25% (Figure 9(b1,b2)). MF's Because the salt drying farm and aquaculture ponds are similar in shape, the traditional shape parameter extraction of aquaculture ponds encountered troubles from the joint operation of the saltworks and aquaculture. However, different from the water used for aquaculture, the raw material of the salt farm is seawater without chemical agents. The water quality parameters can be used as a classification basis to monitor the difference between the two, regardless of the quality of the water object. Secondly, because the water quality characteristics of aquaculture ponds are similar, MF has the potential to correctly classify mixed water objects under the effect of water quality characteristics. Finally, the shadow noise of urban buildings affects the effect of NDWI in extracting water bodies [64], MF is not affected by urban shadow noise. Then MF still shows instability in some aspects.
(1) Although the classification of mixed objects has been improved, the probability of correct classification is still within a large uncertainty range (25-100%). Some natural lakes are also classified as aquaculture ponds because of their probability of being classified as aquaculture ponds, making the overall pond area too large; (2) The sample data set is also an important factor affecting the MF's classification of aquaculture ponds. When the surrounding landscape is complicated, the error is particularly obvious. Therefore, it may be a better choice to allocate the number of different types of sample points based on the area ratio [65].
Our results also demonstrated advantages in three major deltas of China. In the Yellow River Delta, where saltworks are distributed, MF can more accurately realize the classification of aquaculture ponds and saltworks. The probability that the ponds in the saltworks are classified as aquaculture ponds is less than 25% (Figure 9(b1,b2)). MF's unstable classification of mixed objects has been relieved here, but it still exists. The probability that some mixed objects are classified as aquaculture ponds is less than 50% (Figure 8(b1)), but most are in the 50-75% range. (Figure 9(b1)). The Yangtze River Delta has the highest river network density in China, with an average river network length of 4.8-6.7 km per square kilometer. There are more than 200 lakes, and aquaculture ponds are operated around the inland lakes. The classification results of MF under the landscape surrounding inland lakes were still effective, and the probability of classifying small lakes and diversion ditches as aquaculture ponds is less than 50%, and the probability of aquaculture pond objects being correctly classified is mostly above 75% (Figure 9(c2,c3,c4)). Due to the different morphological characteristics of small aquaculture ponds in the Pearl River Delta, the sample dataset could not represent all aquaculture ponds. The probability interval for unrepresented objects to be classified as aquaculture ponds was usually between 50 and 75%, and some are between 25% and 50% (Figure 9(a1)). This was also an important source of uncertainty in model classification. However, MF had a better classification effect on rivers (Figure 9(a1)), ditches (Figure 9(a1)), and wetland parks (Figure 9(a2)). The probability of these non-aquatic objects being classified as aquaculture ponds is usually less than 25%. small lakes and diversion ditches as aquaculture ponds is less than 50%, and the probability of aquaculture pond objects being correctly classified is mostly above 75% ( Figure  9(c2,c3,c4)). Due to the different morphological characteristics of small aquaculture ponds in the Pearl River Delta, the sample dataset could not represent all aquaculture ponds. The probability interval for unrepresented objects to be classified as aquaculture ponds was usually between 50 and 75%, and some are between 25% and 50% ( Figure  9(a1)). This was also an important source of uncertainty in model classification. However, MF had a better classification effect on rivers (Figure 9(a1)), ditches (Figure 9(a1)), and wetland parks (Figure 9(a2)). The probability of these non-aquatic objects being classified as aquaculture ponds is usually less than 25%.

Implications for Further Development of the Extraction Approach
The MF classification framework developed in our study presents a baseline for aquaculture detection across various spatiotemporal scales. Currently, MF has been applied to annual composite images, yet the accuracy of such implementation can be affected by the impoundment of new ponds within a year and by the seasonal variation of combining rice and shrimp in the ponds [66]. Thus, aquaculture pond detection at a higher temporal resolution (i.e., seasonal) is needed to alleviate the omission and commission errors. Further warranted development includes near real-time monitoring of aquaculture ponds by leveraging high-resolution images from multiple satellite missions

Implications for Further Development of the Extraction Approach
The MF classification framework developed in our study presents a baseline for aquaculture detection across various spatiotemporal scales. Currently, MF has been applied to annual composite images, yet the accuracy of such implementation can be affected by the impoundment of new ponds within a year and by the seasonal variation of combining rice and shrimp in the ponds [66]. Thus, aquaculture pond detection at a higher temporal resolution (i.e., seasonal) is needed to alleviate the omission and commission errors. Further warranted development includes near real-time monitoring of aquaculture ponds by leveraging high-resolution images from multiple satellite missions such as Sentinel-2, Landsat, and CubeSats (e.g., Planet Dove) [67]. The associated products can accommodate agricultural and aquatic ecosystem management across spatial scales. It is also worth noting that our framework can serve as a benchmark for deep learning algorithm development such as the convolutional neural network, which, similar to MF, also detects targets based on their colors and geometries.
The developed data product from our work can provide an essential base map for various downstream applications. The water quality condition of aquaculture ponds is of great importance for fish production. Using the pond extent from our product, several water quality indicators can be derived at the field scale using multi-spectral satellite images including water surface temperature [68] and other optically active constituents [48]. This additional information can assist precision aquaculture as well as water quality management in connected water bodies that may suffer nutrient enrichment and other adverse impacts from aquaculture [69].

Uncertainties and Limitations
The aim of this work is to investigate the role of water quality parameters in aquaculture pond extraction. In order to test such a general hypothesis, we have focused on eight regions with typical characteristics. This approach allows the illustration of the effects of adding water quality parameters, but has limitations in the spatial validation along the coast of China. Despite these limitations, the advantages of adding water quality parameters to aquaculture pond extraction are clear. In future, the statists of annual aquaculture ponds area should be validated by county and species (i.e., fish, shrimp, and mussels) to reduce uncertainties, spatially and temporally.

Conclusions
Our study tested the hypothesis that water quality features serve as an important complement to shape features, which ensures that water bodies that are confused with the shape of aquaculture ponds (e.g., saltworks, small reservoirs, and rice fields). Results of testing at eight sites along the coast of China showed that the combination of water quality characteristics and six transfer characteristics improved the accuracy of distinguishing aquaculture ponds from salt pans, rice fields, and wetland parks, which typically had F1 scores > 85%. Several key conclusions were made: 1.
The sharpened NDWI image improves the performance of water-object masked in water body extractions.

2.
The combination of water quality parameters and water body geometry enables aquaculture ponds from other regular water objects to be efficiently distinguished.

3.
Coastal aquaculture in three major river deltas (i.e., Yellow River, Yangtze River, and Pearl River Deltas) shows distinct morphological characteristics and competing land use.
Distinguishing different types of surface water bodies is of great significance for clarifying the water cycle mechanism, explaining the geo-biological-chemical circulation, and managing the coastal water ecological environment. Due to the generality of water quality inversion for different sensors, our method can be applied to other medium-resolution remote sensing datasets. Since aquaculture ponds are seasonally filled by biological and chemical active water, more work is necessary to refine the extraction approaches of different types of aquaculture ponds by the integration of different sources of datasets in the specific environment.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/rs14143306/s1, Figure S1. Wukan Port field survey photos, Figure S2. Rudong field survey photos, Figure S3. Donggang field survey photos, Figure S4. Chengkou Saltworks field survey photos, Figure S5. Hangu Saltworks field survey photos, Figure S6. Number of Sentinel-2 remote sensing images used annually from 2015 to 2020, Figure S7. Availability of Sentinel-2 remote sensing imagery used annually from 2015 to 2020, Figure S8. The water mask made by S1 VH, S2 NDWI, and S2MNDWI at test site 3, and the comparison of the sharpening effect, Text S1. Contain Equation (S1), Table S1. MF optimized from the four models located in the southern site and the southern model, Table S2. MF optimized from the four models located in the northern site and the northern model, Table S3. The area of each year proposed by this research algorithm is compared with the pond area in the China Fishery Statistical Yearbook.

Data Availability Statement:
The datasets generated during the study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.