The topographic variability in a mountain ecosystem forms a micro-habitat with diverse microclimates influencing the variations in local and regional vegetation patterns thereby supporting high levels of biodiversity [5]. In recent years, international research on mountain vegetation has broadened from earlier studies that focused mostly on temperate zones of Europe (e.g., the Alps and the Scandes) [8] to other mountainous regions in the world including Alaska and the State of California, and the Rocky Mountains in North America [9,10], tropical mountainous regions of Africa (Kenya, Mount Kilimanjaro, Mount Wilhelm in New Guinea, etc.) [11,12], the subtropical Andean regions [13,14], and the southern and eastern Himalayas [15]. Meanwhile, the research field has also expanded from the physiological ecology of alpine plants to issues of elevational zonality and the impact of global climate change on alpine vegetation patterns [16,17]. In particular, studies of alpine vegetation in the eastern Himalayas and across the QTP have made great progress in recent years through the efforts of scholars in China, which not only studied the vegetation pattern in the Three River Headwater Region, Hengduan Mountains, but also included high latitudes, such as the Tianshan Mountains, Changbai Mountains, and Qinling Mountains, and mainly deliberated the response mechanisms for conservation of mountain vegetation in the face of global climate change and human activities, in particular those major ecological projects [18–21].

Several environmental factors have been widely identified as drivers for vegetation cover changes. Natural (e.g., climatic factors, topography and etc.) and social factors (e.g., socioeconomics, ecological engineering policies and settlements) can have a strong influence on vegetation growth processes, the extent of which varies from one place to another [5,22]. Most studies believe that precipitation and temperature are the main climatic factors affecting regional vegetation change [23], in addition to topography [24]. In mountainous areas with complex topography, topographic features (elevation, slope, aspect, etc.) have noticeable effects on vegetation patterns by themselves but also by controlling other environmental factors, such as solar radiation, wind, precipitation, snow cover, and edaphology distribution, a combination of these environmental factors jointly determines the spatial heterogeneity of vegetation pattern [17,19]. Besides, human activities (such as roads construction, settlements, and hydraulic engineering) also have an impact on mountain vegetation cover [22], which may have a negative impact on vegetation in low elevation regions, while less interference on the vegetation in the middle elevation [7,24]. At present, the analysis of drivers affecting vegetation patterns in the Tibetan Plateau region is mostly concentrated in the analysis of climatic factors, and relatively little research has been conducted on topographic control mechanisms [25–27]. In the context of global climate change, topography is the only relatively constant environmental factor, and deeper insights on how it controls the vegetation change pattern is particularly necessary to understand the vegetation dynamics and strengthen the response analysis of vegetation growth and distribution to topographic factors.

Given their inaccessibility in most mountainous regions, both remote sensing (RS) and geographic information system (GIS) technology stand out jointly as powerful tools for monitoring vegetation cover changes in mountainous areas by providing continuous, spatially detailed satellite data on mountainous vegetation cover [28]. With the development of RS technology, the types of remote sensing sensors have become more and more diversified, and the large-scale long time series vegetation dynamics has gradually

developed into a hot spot for global change research, among which MODIS vegetation cover product data is regarded as one of the most effective data products for vegetation productivity analysis [29]. NDVI has been proved to be a comprehensive index to describe the ecological functional characteristics of vegetation growth, net primary productivity and phenology at regional, continental, and global scales [22,30]. Thus far, it has been widely used to detect the response of vegetation dynamics to climate change, human activities and other driving factors at multiple spatio-temporal scales [31,32]. However, many drivers are known to exhibit spatial non-stationarity with distinct spatio-temporal characteristics [33], while geographical weighted regression (GWR) model has improved the traditional model by combining spatial correlation with linear regression to better reflect the spatial non-stationarity of the relationship between NDVI and drivers by calculating the local parameters of the regression model [34].

The TPRR, located in the heart of the Hengduan Mountains in the southeastern Tibetan Plateau, has been paid widespread attention for being a world-renowned scenic spot and World Natural Heritage Site. The region is characterized by paralleling alternately high mountains and canyons, and the influence of topographic elevation differences has resulted in significant vertical temperature variation and contributed to the rich biodiversity of the region. Mountain ecosystems are signal "amplifiers" of biological responses to global change [4,35], and exploring the unique geographical environment and vertical zone spectrum of mountain bioclimates favors our understanding of regional biodiversity, ecosystem functions and services, vegetation restoration and reconstruction, as well as other issues. Existing studies from the TPRR mainly focus on geology, dynamic changes of alpine vegetation, land use/cover changes, and biodiversity conservation [36–39], and tourism resources [40–42]. Although a few scholars have studied vegetation changes and its driving factors analysis in the whole basin of the TPRR and various sub-basins (e.g., Lancang river and Jinsha river basins and etc.) [43], both of the vegetation dynamics in the Jinsha River Basin and the Lancang River Basin are jointly affected by temperature and precipitation, and the influence of precipitation is more significant in the Jinsha River Basin, while in the Lancang River Basin is more sensitive to temperature [44–46]. Furthermore, anthropological factors such as the establishment of large dams and reservoirs, hydropower stations, and the implementation of ecological projects have great and obvious specific effects on the vegetation pattern of the basin [2]. Only a few previous studies highlight the relationship between vegetation change and meteorological factors in the source areas of the Three Parallel Rivers which confirmed a significant correlation of vegetation cover with air temperature in this area [43–46]. The TPRR is one of the important ecological barriers and water resources conservation areas in China with unique arid valleys and distinct topography differences. Therefore, it is of vital importance to understand the influence of topographic factors on vegetation cover in this typical geographical unit.

Based on the MODIS NDVI image data from 2000 to 2019, combined with DEM and vegetation type data for trend method and geographically weighted model, our study discusses the spatio-temporal dynamic changes of vegetation and its response to topographic factors (i.e., elevation, slope, aspect, and other factors) in the TPRR. In terms of the spatial heterogeneity, the difference of vegetation restoration effect and its causes, and multiple drivers were analyzed by using the geographically weighted model. The specific research questions addressed in this study include: (1) What are the dynamic characteristics of NDVI at spatial and temporal scale? (2) How do different topographic factors effect on vegetation pattern? (3) What is the overall performance of all topographic factors by generalized analysis in terms of Geographically Weighted Regression modelling?

#### **2. Materials and Methods**

#### *2.1. Study Area*

The TPRR (24◦00'~36◦00' N, 90◦20'~102◦20' E) is located in the Hengduan Mountains in southwest China, spanning four provinces/Autonomous Region, Qinghai, Sichuan, Yunnan, and Tibet, bordering Qinghai-Tibet provincial boundary in the north, SichuanYunnan provincial boundary in the east, and China-Myanmar national boundary in the west, with a total watershed area of about 50.07 million km2 (Figure 1). In this paper, the division of the study area was obtained based on the DEM data and watershed data provided by the Resource and Environment Science and Data Center of the Chinese Academy of Sciences (available at https://www.resdc.cn accessed on 12 August 2020). Coupled with the existing reference, a series of spatial data analyses were conducted using DEM data, including fill, flow direction, and flow accumulation, which ultimately determined the watershed extent of the TPRR [47–49]. Three major Asian rivers, the Yangtze (Jinsha River), the Mekong River (Lancang River) and the Salween River (Nujiang River), run parallel through the longitudinal valleys from north to south, with a vertical elevation difference of nearly 6000 m [47]. The unique environmental gradients make the region rich in biodiversity and one of the most unique areas of species diversity in China [50]. The climate type in the basin is diverse, including cold, temperate, subtropical, and tropical zones from north to south, with temperature and precipitation decreasing from south to north with increasing elevation. The variable climate types and geographical environment make the vegetation distribution in the basin have obvious latitudinal and vertical zonality characteristics (Figure 1).

**Figure 1.** The location elevation and vegetation type maps of the study area. (**a**) the location of TPRR at the eastern part of Hindu-Kush-Himalayan Region; (**b**) Elevation gradients across the TPRR; (**c**) the vegetation types of TPRR.

#### *2.2. Data Sources and Preprocessing*

#### 2.2.1. Remote Sensing Data

NDVI data were provided by the National Aeronautics and Space Administration (NASA) MODIS Terra (MOD13Q1) 16-day vegetation index product (2000–2019) at a spatial resolution of 250 m (available at https://ladsweb.Modaps.eosdis.nasa.gov/search/order/, accessed on 12 August 2020). MODIS Reprojection Tool (MRT) was used to pre-process the original data (e.g., format and projection conversion, clipping, etc.). Besides, the maximum NDVI value was calculated using the maximum value composition method (MVC), and we obtained the optimal vegetation cover for the study area and further reduced the influence of cloud cover and atmospheric scattering. Finally, regions with NDVI value greater than 0.05 were identified as vegetated areas to avoid the interference of underlying surface information on NDVI in areas with low vegetation cover.

#### 2.2.2. DEM Data

DEM data at a spatial resolution of 90 m × 90 m in the TPRR were collected from the SRTM dataset (https://srtm.csi.cgiar.org/srtmdata/, accessed on 12 August 2020). The primary topographic features based on DEM data, such as elevation, slope, and aspect layers were computed using the "Slope/Aspect Tool" in the ArcGIS 10.8 Spatial Analyst module. According to the Chinese classification standards for potential hazard of soil erosion and some other relevant references [51–53], they were further divided into several categories (Table S1), and three maps of elevation, slope and aspect were generated (Figure 2). In addition, spatial analysis tools were also used to overlay both the topographic factors layers and NDVI data, to obtain the spatio-temporal distribution of NDVI for each topographic factor over the past 20 years.

**Figure 2.** The hierarchical spatial distribution of topographic factors in the study area. (**a**) elevation gradients, (**b**) slope distribution, and (**c**) aspect classification across the TPRR.

#### 2.2.3. Vegetation Type Data

ESA CCI-LC land cover types were provided by the European space agency (http://maps. elie.ucl.ac.be/CCI/viewer/, accessed on 12 August 2020) with a 300 m spatial resolution. With more than 70% overall accuracy, the global ESA CCI-LC data has been widely used in land use/land cover dynamic studies, around the world monitoring [54–56]. The vegetation types of the region comprise of cropland, grassland, shrub grass, evergreen broad-leaved forest, deciduous broad-leaved forest, evergreen coniferous forest, and others (mainly water, bare land, snow, and ice, etc.), a total of seven categories [57,58]. Referring to relevant literature [59], the areas with unchanged vegetation types from 2000 to 2015 were extracted here to represent the vegetation cover status of the study area in recent years, so to minimize the influence of land use change on vegetation type distribution.

In this paper, we resampled all the remote sensing data with 300 m resolution to facilitate correlation analysis based on pixel scale.

#### *2.3. Methods*

#### 2.3.1. Trend Analysis and Significance Test

Firstly, we obtained the spatial distribution of the regional annual mean NDVI using MVC method for trend analysis and significance test. Then, we classified NDVI into six classes according to previous literature and actual vegetation conditions in the region [51], i.e., non-vegetation cover (NDVI < 0.05), low vegetation cover (0.05–0.30), medium and low vegetation cover (0.30–0.45), medium and high vegetation cover (0.45–0.60), medium and high vegetation cover (0.60–0.75), and high vegetation cover (NDVI > 0.75).

Theil-Sen + Mann-Kendall method has been widely used in meteorology and hydrology to get the variation characteristics of NDVI analysis [30,60], which has advantages to avoid outlier data disturbance or measurement error. The formula is:

$$S\_{NDVI} = median \frac{\mathbf{x}\_j - \mathbf{x}\_i}{j - i} \qquad (j > i > 1) \tag{1}$$

where, *xi* and *xj* denote the *i-th* and *j-th* year time series data. *SNDVI* >0 indicates that the vegetation shows an upward trend and vice versa indicates a downward trend. Meanwhile, Mann-Kendall method can test whether the trend is significant and find the abrupt point of NDVI changes, which formula is as follows:

$$Z = \begin{cases} \frac{\frac{S}{\sqrt{Var(S)}}}{\sqrt{Var(S)}} & (S > 0) \\ 0 & (S = 0) \\ \frac{S+1}{\sqrt{Var(S)}} & (S < 0) \end{cases} \tag{2}$$

Among them,

$$S = \sum\_{i=1}^{n=1} \sum\_{j=i+1}^{n} \text{sign}(\mathbf{x}\_{j} - \mathbf{x}\_{i}), \quad \text{sign} = \begin{cases} 1 & (\theta > 0) \\ 0 & (\theta = 0) \\ -1 & (\theta < 0) \end{cases} \tag{3}$$

where, *Z* is a standard normal distribution, and *S* is the test statistic. *xi* and *xj* are time series data; *n* is the sample number; When *n* ≥ 8, *S* is approximately normal distribution, and the variance calculation formula is:

$$Var(S) = \frac{n(n-1)(2n+5)}{18} \tag{4}$$

*Z* is a standard normal distribution, when |*Z*| > *Z1-α*/2, it showed significant change trend. Where, *Z1-α/2* is the corresponding value of the distribution table of the standard normal distribution function at the confidence level α. The trends could be estimated by combining Sen's slope and the Mann–Kendall test in light of the NDVI series. Given confidence level *α* = 0.01 and *α* = 0.05, whether the trend changed significantly depended on |*Z*|, which can be divided into six classes: at 0.01 confidence level, (1) very significant degradation (*SNDVI* < 0, |*Z*| > 2.58); (2) very significant improvement (*SNDVI* ≥ 0, |*Z*| > 2.58); at 0.05 confidence level, (3) significant degradation (*SNDVI* < 0, 1.96 < |*Z*| < 2.58); (4) significant improvement (*SNDVI* ≥ 0, 1.96 < |*Z*| < 2.58); (5) non-significant degradation (SNDVI < 0, |*Z*| ≤ 1.96); (6) non-significant improvement (SNDVI ≥ 0, |*Z*| ≤ 1.96).

#### 2.3.2. Geographically Weighted Regression Model

GWR model can analyze the spatial nonstationary characteristics of the data and explore the spatial heterogeneity between vegetation changes and its driving factors [6,61,62]. The kernel and bandwidth of the GWR model, which are key parameters affecting the accuracy of the model, are now generally determined by Gaussian kernel function, and the optimal bandwidth is generally determined based on the estimated distance of the spatial autocorrelation of the dependent variable [63]. In addition, variance Inflation Factor (VIF) was generally used to eliminate highly correlated variables from the model.

Ordinary Least Squares regression (OLS) model is used to check for the presence of multicollinearity between variables to improve the accuracy of GWR model [34,64]. Compared with OLS model, GWR model is an extension of the traditional global logistic regression model, which includes spatial factors and geographic location information. The weighted least square method is used to estimate the parameters of each sample point, so that each of them has a corresponding estimation coefficient. In GWR4.0 software, the GWR model was constructed with the same variables as the OLS model [65,66], and the optimal bandwidth of 102 was determined according to the minimum *AIC* value. Besides, we divided the region into a 10 km × 10 km grid and extracted NDVI trend slope (*SNDVI*), elevation, slope, and aspect values for each grid-cell [66,67]. Assuming that the dependent variable (Slope of NDVI variation trend from 2000 to 2019 is represented by *SNDVI*) is y, and the elevation, slope and aspect value are x1, x2, and x3 respectively, the GWR model can be expressed as:

$$y = b\_0(\mu\_i, \upsilon\_i) \ + \sum\_{i=1}^k b\_i(\mu\_i, \upsilon\_i) \mathbf{x}\_{ij} + \varepsilon\_i \tag{5}$$

where (*ui, vi*) is the geographic coordinate center of a certain region, b0 is a constant, b1 is the regression coefficient of independent variable, and *xij* is the independent variable. Figure 3 summarizes the analysis procedures.

**Figure 3.** Conceptual framework of study design.

#### **3. Results**

#### *3.1. Spatial and Temporal Distribution Patterns of NDVI*

Inter-annual variation analysis results revealed that regional NDVI value fluctuated between 0.56 and 0.61 in the past 20 years, with an overall growth rate of 0.007/10a (*p* > 0.05), and only the Lancang River basin showed a significant upward trend (*p* < 0.05) (Figure 4). Additionally, the results of Mann-Kendall test showed an abrupt change of NDVI in the whole basin occurred around the year 2002, 2004, and 2008, respectively. Combined with the abrupt changes of NDVI in each sub-basin, these changes occurred almost simultaneously around 2004.

The mean NDVI in the basin tended to increase and then decrease with increasing elevation, and its value reached the maximum at 3000 m (about 0.80) and the minimum at 5500 m (about 0.05) (Figure S1a). Vegetation in the area below 2500 m was dominated by high vegetation cover, mainly in the southern part of the region, and the vegetation types were mainly cropland and evergreen coniferous forest (Table S2). Besides, the spatial distribution of mean NDVI in different elevation zones showed that regional vegetation was affected by elevation showing significant north-south differences (Figure S2a). The proportion of area occupied by each class indicated that the vegetation cover of all elevation zones is dominated by high vegetation cover, especially in the southern regions of the TPRR and in its core area (Figure 5). However, as the elevation continued to rise, NDVI started to decline from 4500 m, and the vegetation cover gradually transitions from low vegetation cover to unvegetated snow-capped mountain tops.

**Figure 5.** Area statistics of multi-year mean NDVI for topography factors during 2000–2019.

With the increase of slope, NDVI showed an increasing trend with gradually slower growth rate (Figure S1b). The mean NDVI with slope <5◦ was the lowest (0.40) and the value was between 0.05 and 0.3, which was mainly low vegetation cover type dominated by grassland and distributed in the river source area of the basin (Figure 5, Figure S2b, Table S2). The NDVI for regions within the slope range of 5–8◦ ranges from 0.6–0.75, which was mainly middle and high vegetation types. NDVI increased with the increasing of slope grades and each slope is gradually dominated by high vegetation cover types, whose distribution gradually shifted to the dry-hot valley and eventually reached its maximum value (about 0.80) at slope >35◦.

Furthermore, Figure S1c showed the NDVI changes with aspect from shady slope to sunny slope that the mean NDVI values were uniformly distributed in all aspects (approximately 0.60), and the smallest (about 0.01) existed on flat land. The proportion of high vegetation cover in each aspect was the largest (Figure 5), which distribution pattern was relatively consistent (Figure S2c). The brief descriptions were the proportion of high vegetation cover in semi-shady slope was the largest (35.25%), followed by semi-sunny slope, while the proportion of sunny slope was the smallest (27.83%) (Table S2).

#### *3.2. Influence of Different Topographic Factors on the Variation Trend of NDVI*

There was a significant spatial variation in NDVI of the TPRR from 2000 to 2019, showing that the vegetation cover in the basin had slightly improved (Figure 6). However, 31.74% of the region showed no changes and 20.28% areas indicated slight degradation. Overall, the areas with improvement were mainly concentrated in the source area of Jinsha river basin in the north and the southern part of the TPRR. Meanwhile, degradation trends were primarily observed in north of the Jinsha river basin and the valley of the TPRR as well as other areas with significant response to vegetation growth changes.

**Figure 6.** Spatial distribution of (**a**) NDVI trends and (**b**) significance during 2000–2019.

Most of NDVI in the TPRR showed increasing trends under different topographic factors (Figure S3). Among them, NDVI value at all elevation zones showed increasing trends except for regions above 5500 m (Figure S3a), and we also observed that both slope and aspect did not play any significant role in determining NDVI trend distributional patterns (Figure S3b,c).

More specifically, except for the three elevation zones, 3500–4000 m, 4000–4500 m, and above 5500 m, with the largest proportion of no changes at 34.36%, 35.04%, and 43.81%, respectively, and the majority of elevation zones (mainly 1000–3500 m) showed a slight improvement trend (Figure 7a, Table S3).

**Figure 7.** NDVI trends area statistics for different topography factors from 2000 to 2019. (**a**–**c**) represent elevation, slope and aspect, respectively.

For different slopes, approximately 35% of the area under all slope grades showed a slight improvement (Figure 7b, Table S3). The proportion of vegetation improvement area increased initially and eventually decreased with the steeper slope. Among them, the slope in the northern part of the basin and the source area of TPRR was relatively small (<5◦), while the vegetation improvement trend was obvious, with grassland as mainly vegetation type (Figure S4a). The slope in the valley was mainly between 8–35◦, which dominated by grassland and evergreen coniferous forest. The slope more than 35◦ in a few river valleys, especially in the Yunnan section of the TPRR, was mainly evergreen coniferous forest. Besides, the proportion of improved areas (37.12%) was the largest in regions with 15–25◦ slope. However, we observed that the proportion of area with no changes decreased gradually with the increasing slope (Figure S4b).

According to different aspects, although there was little difference in the proportion of areas between vegetation improvement, stability, and degradation, on an average we could witness slight improvement (Figure 7c, Table S3). Given that aspect had little influence on the vegetation variation trend in the basin, particularly in flat land, the area proportion of vegetation improvement had an inconspicuous increasing trend from shady slope to sunny slope (Figure S4c).

#### *3.3. GWR-Based Driving Force Analysis of Vegetation Change among Different Basins*

Taking *SNDVI* as the dependent variable, and topographic factors (elevation, slope, and aspect) as independent variables, the traditional multiple linear regression model based on OLS could not well explain the spatial heterogeneity of vegetation changes and topographic factors (Table S4; Figure S5; OLS model, *R*2 = 0.03, *F* = 60.76), given variance inflation factors (VIF) of each variable were all less than 10.

However, GWR model provided a better overall fitting than OLS (Table 1, *R*2 = 0.33), which is more proper than OLS in explaining vegetation changes. Moreover, the *AICc* difference between GWR and OLS was much larger than 3 (Δ*AICC* = 1458.77), indicating that the fit degree of GWR model is significantly higher than that of OLS, but the Residuals SS of the GWR model also decreases, which indicates that the fitting effect of the GWR model is greatly improved compared with the OLS model (Residuals *F* = 4. 31) (Table S5).

**Table 1.** Parameter estimation and test results of the GWR model.

The elevation in the headwaters of the TPRR (in particular, Nujiang river and Jinsha river basins), the middle reaches of the Lancang river, and a few areas in the south of the TPRR (mainly in the south of the Jinsha river) has a significant positive correlation with *SNDVI*, while the remaining vast majority of areas showed a negative correlation, e.g., the middle reaches of Nujiang river and Jinsha river basin (Figure 8a).

Slope also has a certain regularity on vegetation change since the slope steepness affects the soil and water conservation so that will directly affect the plant growth and development on the earth surface. Slope has location-specific effects positively or negatively on vegetation change in the study area, given whose influence intensity is much smaller than that of elevation. Among them, positive correlation trend mainly existed in the headwaters of the Jinsha river basin and a few areas in the middle reaches, and the lower reaches of the Jinsha river and Nujiang river basins. However, there is a significant negative correlation trend between slope and vegetation variation trend in the junction of the Jinsha river and Nujiang river basins and a small part of the northwest of the source area in Jinsha river basin (Figure 8b).

**Figure 8.** Spatial distribution maps of regression coefficients and significances for all independent variables. (**a**–**c**) are the spatial distribution maps of the local regression coefficients and the significances for all independent variables, e.g., (**a**) elevation, (**b**) slope and (**c**) aspect respectively.

In terms of aspect, some significant positive correlations exist in a few source areas of the Nujiang river basin, a small part of the TPRR, as well as the junction of the south lower reaches between Jinsha river basin and Nujiang river basin given a significant negative correlation presented in some part of the Nujiang river basin and the source areas of the Jinsha river basin (Figure 8c).

#### **4. Discussion**

#### *4.1. Temporal and Spatial Distribution of Vegetation Patterns*

Our results show that from 2000 to 2019, the vegetation in the TPRR shows an improving trend in a majority of the areas, while a few areas remain degraded, especially in the source areas of the entire Three Parallel Rivers basin where degradation is more obvious. Compared with the degraded areas, the proportion of areas with improved vegetation is much larger (45.79%). Our findings are consistent with those of other scholars in the QTP or TPRR, who also emphasized that vegetation dynamics in most areas of the QTP have an upward trend, while some areas such as the Nagqu, a few areas in Qinghai Province and several areas on the southeastern edge of the QTP including the Hengduan Mountains are at degradation trends. Natural environmental changes and anthropogenic disturbances jointly increased vegetation degradation [18,24]. Previous studies of vegetation pattern changes in the TPRR and showed that the significant reduction of vegetation cover was mainly in the source area of the Nujiang river and Lancang river as well as some areas in the central part of the study area, which corroborates with our study's results [20,43]. Different vegetation types of mountainous vegetation are also significantly influenced by spatial scales [47,68,69]. With higher elevations in its northern region and relatively lower elevation in the southern regions the geographical structure in TPRR supports better vegetative growth in the southern region, with cropland and evergreen coniferous forest being main vegetation types. Most of the vegetation is distributed within the elevation range of 1500 to 5500 m, with a myriad of landscape and vegetation types overlaid and intertwined with one another [23]. Along an increasing elevational gradient, the vegetation types transition from subalpine forests, low scrub, alpine grassland to unvegetated areas. This variation in vegetation types is driven by topographic factors that alter regional hydrothermal conditions which in turn indirectly affect vegetation distribution [20,46]. Peng et al. [24] showed the zonal characteristics of vegetation dynamics of different vegetation types, specifically the severe degradation of forests in the southeastern Tibetan Plateau, the increase of alpine grassland vegetation cover in some humid areas of southern Tibet, and the rising trend of cropland under the influence of agricultural activities in the Ali region and the plateau hinterland, all of which undoubtedly confirmed the complex vegetation changes in TPRR.

#### *4.2. Spatial Heterogeneity of Topographic Factors on Vegetation*

The results of GWR model strengthen the intuitive perception of the spatial heterogeneity of different topographic factors on NDVI variation in the TPRR, with results indicating that elevation itself has the greatest influence on local vegetation patterns, followed by slope and aspect, which is consistent with previous studies [24,26]. Altitude mostly influences vegetation patterns by affecting mountain climate conditions [70,71]. Our study area is in the QTP region with large elevation differences, and the distribution of heat conditions is successively controlled by the law of decreasing elevation and solar radiation differences, while precipitation distribution mainly affected by the trend of water vapor channel orientation, which is mainly related to the distribution of mountains. Therefore, analyzing the influence of topographic factors on vegetation change can better understanding mechanism of the climate-vegetation interaction, which is often ignored in current research [72]. Notably, the relative importance of elevation is generally reflected in its effect on the gradient of species' ecological niches, which determines the distributed of vegetation types in the region and the vegetation growth by representing laterally such environmental factors as slope or the amount of soil moisture [73,74]. Besides, slope can be one predictor variable to investigate plant diversity or capture potential changes in elevation in a region [73–76]. Slope can affect species richness by influencing soil moisture, soil erosion rates, and litter accumulation. However, related studies have shown that alpine treeline elevation changes in the Three Parallel Rivers region of southwest China are less affected by slope, further confirming the validity of our findings that slope has little effect on NDVI trends [39]. Moreover, Yirdaw et al. (2015) suggested that aspect affected woody plant species diversity, especially for woody plants at higher elevations, and that sunny slopes tended to receive higher amounts of solar radiation, which could accelerate greater air and soil temperatures [27]. Nevertheless, our results indicated the effect of aspect on NDVI trend was not significant, and the highest vegetation cover (35.25%) was found on the semi-shady slopes. Few scholars have explored the effect of aspect on vegetation changes and even the distribution of some major species in our study area, which may require further research.

#### *4.3. Analysis of the Other Drivers Influencing Changes in Vegetation Dynamics*

Beyond the influence of topographic factors on the vegetation pattern mentioned in this study, climatic conditions are the most well discussed factors for vegetation change in the region. Studies have shown that temperature and precipitation indirectly affect the productivity of terrestrial ecosystems by altering nutrient availability [17,77], and several scholars have reported the influence of temperature and precipitation these two climatic factors on the photosynthesis and respiration of vegetation changes by influencing soil moisture and microbial activities in the QTP. Obviously, their short-term effect on vegetation degradation is not significant [26,78,79]. However, the direct influence of topography on climatic conditions in mountainous areas was mentioned in the previous subsection, and our study confirms the strong influence of elevation on vegetation distribution, but the relationship between topographic and climatic factors and the joint influence of both on vegetation distribution is lacking in this study area and needs further analysis.

In addition, the impact of human activities on the vegetation pattern should not be underestimated. Hengduan Mountain area in the southeast QTP acts as a critical ecological barrier due to its unique ecological location and unique characteristics. In this region, several key national ecological projects have been implemented in the past to protect its natural forests, returning cropland to forest, biodiversity conservation, and ecological restoration, which have shown positive outcomes, and can also be the main reason for the improvement in NDVI observed in this study. Lu et al.compared the effects of three largescale ecological programs (National Nature Reserves, Three North Shelter Forest Program, and the Natural Forest Protection Program) on Chinese vegetation change and revealed that the effectiveness of ecological restoration projects can vary geographically even under the same incentive policy context [80]. Moreover, several national and provincial protected

areas, including the Three Parallel Rivers National Park, have been established in the TPRR. These protected areas can play a positive role in maintaining biodiversity, enhancing the regulation function of ecosystem, and building a scientific and reasonable ecological space. Grazing intensity is also one of the main driving factors affecting vegetation cover in grasslands in the QTP. Grasslands in QTP region have been considered as natural pastures for grazing by herders for many years [81,82]. Studies have shown that long-term grazing directly affects the vegetation cover and productivity of grasslands, but with the implementation of a series of ecological measures related to grassland restoration, e.g., grazing prohibition, and pastureland rehabilitation, the degradation of grasslands has been mitigated, which may also be one of the main reasons why the vegetation trends in our study area tends to improve [83,84].

#### **5. Conclusions**

In this paper, we explored the influence of topographic factors on the vegetation change patterns in the TPRR of southeastern QTP over the past 20 years. We analyzed the spatial and temporal distribution of vegetation patterns and the magnitude of the driving forces of each topographic factor with the GWR model, and the main conclusions are as follows:

During 2000–2019, the NDVI in the basin showed an overall increasing trend, with the most significant upward trend in the Lancang River basin, and the year of sudden change in NDVI in each basin was around 2004. Besides, the spato-temporal distribution of regional multi-year NDVI was influenced by elevation with significant north-south differences, while slope and aspect had little influence. Moreover, NDVI trends in the basin showed an overall improvement trend, mainly with slight improvement. The improvement areas are generally distributed in the southern section and the source area of the Jinsha River basin, and there is a degradation trend in the river source area of the TPRR. Specifically, the vegetation in all elevation zones was slightly improved, and the area of vegetation improvement increased and then decreased with increasing slope, while aspect has the least effect. The GWR model showed that elevation and NDVI were negatively correlated in most areas, and the effect of slope and aspect on vegetation changes was smaller, which were significantly correlated in the middle reaches of the Nujiang River basin and some parts of the Jinsha River basin. The interactions between topographic and climatic factors and their joint effects on vegetation distribution is lacking and needs further study in TPRR. Furthermore, the impact of human activities on the vegetation pattern should also be highlighted and integrated with climate change and other drivers for comprehensive analysis.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/rs14010151/s1, Table S1: Classification of topographic factors in the study area, Figure S1: Variation characteristics of mean NDVI with different topographic factors, Figure S2: Spatial distribution of mean NDVI with different topographic factors during 2000–2019, Table S2: Annual mean NDVI area statistics for different topographic factors during 2000–2019, Figure S3: NDVI trends during 2000–2019 with different topographic factors in the basin (a, b, c show elevation, slope and aspect respectively), Figure S4: Spatial distribution of NDVI trends for different topographic factors from 2000 to 2019, Table S3: Area statistics of NDVI trends for different topographic factors from 2000 to 2019, Table S4: Parameter estimation and test results of the OLS model, Figure S5: Estimation result of regression coefficients obtained through OLS modeling. (Notes: for the entire regression equation, the adjustment R2 is 0.03, F value is 60.76, and P < 0.01; the two dashed lines denote that the t-value is equal to −1.96 and 1. 96 respectively; at 0.05 significance level, t < −1.96 means a significantly negative correlation, while t > 1.96 represents a positive correlation.), Table S5: Comparison of model performance between the GWR and OLS models.

**Author Contributions:** Conceptualization, J.W., N.W. and C.W.; methodology, C.W., X.C. and Y.W.; software, C.W.; validation, J.W. and C.W.; formal analysis, C.W.; investigation, C.W.; resources, C.W.; data curation, C.W.; writing—original draft preparation, C.W.; writing—review and editing, C.W., J.W., N.N., N.W., X.C., Y.W. and Q.C.; visualization, C.W.; supervision, J.W.; project administration, J.W.; funding acquisition, J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Science Foundation of China (No. 31971436 and<PHONE_NUMBER>5); State Key Laboratory of Cryospheric Science, Northwest Institute of Eco-Environment and Resources, Chinese Academy Sciences (SKLCSOP-2018-07) and China Biodiversity Observation Networks (Sino BON)".

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors would like to thank the researchers who have provided the opensource algorithms, which have been extremely helpful to the research in this paper.

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

#### **References**

## *Article* **Analysis on the Spatio-Temporal Changes of LST and Its Influencing Factors Based on VIC Model in the Arid Region from 1960 to 2017: An Example of the Ebinur Lake Watershed, Xinjiang, China**

**Nigenare Amantai 1,2 and Jianli Ding 1,2,\***

**Abstract:** LST (Land surface temperature) is an important indicator for monitoring dynamic changes in the earth's resources and environment. However, the complexity of obtaining long-term, continuous LST data hinders the development of research on LST responses to meteorological factors or LUCC in areas where data is lacking. The objective of this research was to use the VIC-3L (Variable Infiltration Capacity) based on multi-source remote sensing data to simulate and explore spatio-temporal changes in the LST, to analyze the relationship between the LST and meteorological elements by using cross-wavelet transform (XWT) and wavelet coherence (WTC), the relationship between the LST and LUCC by using three-phase remote sensing images of LUCC. The following results were obtained. The annual average LST of the study area is increasing at a rate of 0.027 ◦C per year. The annual average LST level is relatively high in the central and eastern regions. The average temperature has an important influence on LST, which is mainly reflected in the period scale of 1~4a in 1963–1972, 1980–1996, and 2004–2010. The sharp decline in open shrubs may have exacerbated the increase in LST in the study area. This study provides a scientific reference for studying LST in arid areas.

**Keywords:** land surface temperature; VIC model; wavelet analysis; multiscale analysis; LUCC

#### **1. Introduction**

Land surface temperature (LST) is a typical indicator of the energy flow in the land surface-atmosphere-biosphere interactions. The UN Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report pointed out that the global average surface temperature increased by 0.89 ◦C (0.69~1.08 ◦C) from 1901 to 2012 [1]. The change of LST is an intuitive regional climate response to global climate change. On the one hand, the LST variations will directly affect the soil moisture content and the soil evaporation, thereby indirectly affecting regional vegetation growth, agricultural production, and the water cycle [2,3]. On the other hand, the LST is an important element for predicting drought, climate change [4–6]. Therefore, it is of greater research significance in agriculture, hydrology, ecology, environment, climate, and biogeochemistry [4,7–11].

Research on changes in LST is more important in arid regions. One of the most prominent features of arid regions is the low atmospheric precipitation [12], which directly affects the soil moisture content and indirectly affects the vegetation conditions in these areas [13,14]. An increase in the LST can lead to reduce soil moisture and increase dust emission [15]. LST can affect wind erosion through changes in soil moisture content, as temperature increases, the evapotranspiration from the land surface increases, and the soil moisture content decreases [16]. As a result, the adhesion between the soil particles is reduced, and the soil particles are more easily separated from the soil, thereby increasing

**Citation:** Amantai, N.; Ding, J. Analysis on the Spatio-Temporal Changes of LST and Its Influencing Factors Based on VIC Model in the Arid Region from 1960 to 2017: An Example of the Ebinur Lake Watershed, Xinjiang, China. *Remote Sens.* **2021**, *13*, 4867. https://doi.org/ 10.3390/rs13234867

Academic Editors: Zhongwen Hu, Jingzhe Wang, Yangyi Wu and Jie Zhang

Received: 18 October 2021 Accepted: 24 November 2021 Published: 30 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the concentration of the particles in the atmosphere [17]. Moreover, LST has shown a positive relationship with changes in the frequency of wind erosion events in Asia [18]. Therefore, studying the temporal and spatial changes of the LST and its driving factors is helpful to better understand the evolution of drought and water resources and provide theoretical support for sustainable development.

Under the global warming background, experts and scholars in different fields are paying more and more attention to changes in LST. The evolution of the global LST trend in the past century is diagnosed using the spatial–temporally multidimensional ensemble empirical mode decomposition method, the fastest warming in recent decades (>0.4 K per decade) occurred in northern mid-latitudes [19]. It also showed a particularly enhanced trend of LST in semi-arid regions [20]. However, the existing studies on LST dynamics mainly focused on the depiction and quantitative temporal comparison of LST spatial distribution or LST time series analysis over time [21]. The relationship between LST and other climate elements is less involved. Moreover, there are still some difficulties in the long-term large-scale continuous simulation of LST in data-scarce areas.

Exploring the mechanism of interaction between LUCC and local climate has become a popular research topic [22–24]. In the past few years, many studies based on the impact of LUCC on LST have focused on the urban heat island (UHI) [25–27]. The increase of the impermeable layer increases the absorption of solar radiation to a certain extent, thereby increasing the LST [28]. However, each land cover type possesses unique qualities in terms of energy radiation and absorption, in addition, LUCC also alters greenhouse gas emissions, such as CO2, in the atmosphere [29]. Although the urbanization process in some areas is not fast, the effect of LUCC on LST is worth studying, especially for arid areas.

The Land surface model (LSM) is a key tool to predict the terrestrial hydrosphere and its atmospheric coupling and to understand this in-depth. Distributed LSM predicts hydrological states and fluxes, such as LST, at each grid cell [30]. The LSM based on multisource remote sensing data can carry out large-scale simulations of LST to compensate for the discontinuity in time and space caused by only using remote sensing images and station data. The VIC (Variable Infiltration Capacity) model is one of the core application models of the GLDAS and NLDAS [31]. Therefore, it is effective to use VIC to carry out large-scale and long-term continuous LST simulations in scarce data areas.

The main objectives of this study are (1) to simulate daily LST from 1960 to 2017 by using multi-source data and to verify the accuracy of the model from both space and time; (2) to understand temporal and spatial evolution characteristics of LSTs and explore the multiscale, abrupt and periodic characteristics of LSTs in the time series; and (3) to analyze the impact of LUCC on LSTs. The results of this study can provide a reference for hydrological simulation research in arid areas and drought monitoring and management for meteorological delays and other aspects.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Ebinur Lake Watershed (43◦38 ~45◦52 N and 79◦53 ~85◦02 E) is located in the arid area of Northwest China (Figure 1). The study area has a typical continental arid climate, characterized by dramatic annual temperature change, scarce precipitation, and strong potential evaporation [32,33], and the lake basin is surrounded by mountains on the south, west, and north sides, which block outside airflow [34]. Climate change easily affects the basin's water resources and ecosystems. The Ebinur Lake watershed is an important ecological barrier for environmental changes in the Junggar Basin in Xinjiang Uygur Autonomous Region (XUAR), which is a typical inland lake basin in an arid region [35]. Some studies found that the LUCC of the watershed has changed greatly [36]. However, the influence of LUCC on LST in the Ebinur Lake Watershed is not yet clear and needs to be further investigated.

**Figure 1.** Location of the study area.

#### *2.2. Data Sources*

The data used in this study can be described as follows: (1) The elevation data are based on the 90-m digital elevation model (DEM) of Shuttle Radar Topographic Mission (SRTM) data and were obtained from the Geospatial Data Cloud Platform of the Computer Network Information Center of the Chinese Academy of Sciences (http://www.gscloud.cn/ accessed on 27 November 2021). (2) The global land cover data, with a resolution of 1 km, were obtained from the University of Maryland (http://app.earth-observer.org/data/basemaps/ images/global/LandCover\_512/LandCoverUMD\_512/LandCoverUMD\_512.html/ accessed on 27 November 2021), and the Terra and Aqua combined Moderate Resolution Imaging Spectroradiometer (MODIS) Land Cover Type (MCD12Q1) Version 6 data from the USGS EROS Center (http://doi.org/10.5067/MODIS/MCD12Q1.006/ accessed on 27 November 2021), has a resolution of 500 m and uses two images (2001 and 2017). (3) The soil data were taken from the 0.5-degree soil data available for the VIC model released by the University of Washington (https://vic.readthedocs.io/en/master/Datasets/Datasets/ accessed on 27 November 2021) and based on the World Soil Database (HWSD) 1:1 million Chinese soil dataset (v1.1). (4) The meteorological data were obtained from the daily dataset of surface climate data in China (V3.0), the time resolution is daily and the time scale is from 1 January 1960 to 31 December 2017. The data were obtained from the China Meteorological Data Network (http://data.cma.cn/ accessed on 27 November 2021). (5) The verification data of LST were obtained from a combined Terra and Aqua MODIS LST and meteorological station data product for China from 2003 to 2017 [37] with a resolution of 5.6 km. The model resolution is 8 km\*8 km, and the basin is divided into 767 grids (Figure 2).

**Figure 2.** Schematic diagram of the model grid.

#### *2.3. VIC Model*

The VIC hydrological model is a large-scale hydrological model with a variable infiltration capacity. The VIC model is an open-source large-scale distributed conceptual hydrological model based on one kind of water and heat balance and physical dynamic mechanisms (VIC's homepage: http://vic.readthedocs.org/en/master/ accessed on 27 November 2021). Specific theoretical references include [38,39].

The energy balance calculation of the VIC-3L model mainly consists of four parts: LST, sensible heat flux, latent heat flux, and geothermal change. The other three parts are determined by LST, and the connecting factor of water balance and energy balance is latent heat flux. In an ideal situation, the energy balance calculation formula of each land cover type in the grid unit is as follows:

$$R\_n = H + \rho\_w L\_\epsilon E + G \tag{1}$$

where, *Rn* is net radiation (W/m2), *H* is sensible heat flux (W/m2), *ρwLeE* is latent heat flux (W/m2), and *G* is surface heat flux (W/m2). The net radiation is calculated by the following formula:

$$R\_n = (1 - a)R\_s + \varepsilon \left(R\_L - \sigma T\_s^4\right) \tag{2}$$

where, *α* is the surface reflectance (%) of a certain land cover type, *RS* is the downward short-wave radiation (μm), *ε* is the specific emissivity(%), *RL* is the downward long-wave radiation (μm), and *σ* is the Stefan-Boltzmann constant (*σ* = 5.67 × 10−8 W/(m2· ◦C4)). The calculation formula of sensible heat is as follows:

$$H = \frac{\rho\_a c\_P}{r\_{\text{fl}}} (T\_s - T\_a) \tag{3}$$

where, *Ts* is the LST (◦C), *Ta* is the atmospheric temperature (◦C), and *rh* is the aerodynamic impedance (s/m). The surface heat flux is estimated by the heat of two layers of soil. For the upper layer (the first layer + the second layer), the total depth is assumed to be *D*1 (m), then its calculation formula is:

$$G = \frac{\kappa}{D\_1}(T\text{s} - T\_1) \tag{4}$$

where, *κ* is the soil thermal conductivity(W/m· ◦C), *T*1 is the temperature of the upper soil ( ◦C). Assuming that the depth of the underlying soil is *D*2 (m) and the temperature at the bottom of the soil is constant, then:

$$\frac{C\_s \left(T\_1^+ - T\_1^-\right)}{2\Delta t} = \frac{G}{D\_2} - \frac{k(T\_1 - T\_2)}{D\_2^2} \tag{5}$$

where, *CS* is the coefficient of soil heat conductivity (J/m3· ◦C), *T*+ 1 and *T*− 1 is the soil temperature (◦C) at the beginning and the end of the period in the upper soil respectively, and *T2* is the constant temperature (◦C) of the lower soil. Currently, it is believed that *κ* and *CS* do not change with the change of soil moisture content. According to Equations (4) and (5), the calculation formula of surface heat flux can be deduced as follows:

$$G\_{(n)} = \frac{\frac{\chi}{D\_2}(T\_s - T\_2) + \frac{C\_s D\_2}{2 \cdot \Delta t} \left(T\_s - T\_1^{-}\right)}{1 + \frac{D\_1}{D\_2} + \frac{C\_s D\_1 D\_2}{2 \Delta t k}}\tag{6}$$

Combined with the above formula, the LST, sensible heat flux and surface heat flux can be calculated through iteration.

#### *2.4. Wavelet Analysis*

#### 2.4.1. Morlet Continuous Wavelet

The long-term LST changes have obvious periodic characteristics, and there is less analysis of the multi-scale time characteristics between LST and meteorological elements. The wavelet analysis method has a powerful multi-scale resolution function, which can detect the response period of meteorological elements and LST, and can also show multiple change cycles hidden in the time series, reflecting the changing trend of the system within the same time range. Therefore, the wavelet analysis is more advantageous when studying the multi-scale characteristics of LST and meteorological elements. In this paper, we used the toolbox provided by Grinsted et al. [40], you can download it from this address: (http://www.glaciology.net/wavelet-coherence/ accessed on 27 November 2021).

The Morlet continuous wavelet function is a commonly used complex wavelet function. It has the advantages of rich phase information, good time aggregation, and highfrequency resolution. Therefore, it is widely used to study the correlation between two non-stationary sequences [41,42]. Therefore, this study selected the Morlet wavelet to analyze the degree of correlation between the LST and each meteorological element. Morlet wavelet is a complex wavelet, and its expressions in the time domain and frequency domain are:

$$
\Psi(t) = \pi^{-1/4} e^{i\omega\_0 t} e^{-t^2/2} \tag{7}
$$

$$
\hat{\psi}(\omega) = \pi^{-1/4}\sqrt{2}\varepsilon - \frac{1}{2}(\omega - \omega\_0) \tag{8}
$$

In the formula, *ω*0 is a constant, and *i* is an imaginary number. Morlet wavelet is a single-frequency negative sine function under the Gaussian envelope. When using wavelets for feature extraction purposes the Morlet wavelet (with *ω*0 = 6) is a good choice, since it provides a good balance between time and frequency localization.

#### 2.4.2. Cross-Wavelet Transform (XWT)

Similar to the Fourier transform, according to the description of the transformed sequence energy in Perseval's theorem, that is, after the signal is transformed from the time domain to the frequency domain, the total energy remains unchanged, then the local time sequence *x*(*t*) The expression of Wavelet Power Spectrum (Wavelet Power Spectrum, WPS) is:

$$(\mathcal{WPS})\_x(a,b) = |\mathcal{W}\_x(a-b)|^2 \tag{9}$$

The above-mentioned continuous wavelet transform is only a one-dimensional sequence of the wavelet transform, and this research mainly involves the correlation of twotime series. In the early research, it is necessary to perform a one-dimensional wavelet transform on each sequence separately, while Hudgins et al. [43], the XWT can not only directly analyze the correlation between two time series at different frequencies, but also reflect their phase information and local characteristics in the time domain and frequency domain. For any two series of time *x*(*t*) and *y*(*t*), the expression of the XWT is:

$$\mathcal{W}\_t^{XY}(\omega) = \mathcal{W}\_t^X(\omega)\overline{\mathcal{W}\_t^Y(\omega)}\tag{10}$$

In the formula, the wavelet transforms of the two time series *x*(*t*) and *y*(*t*) are *WX t* (*ω*) and *WY t* (*ω*) respectively. The corresponding wavelet cross power (Cross Wavelet Power, XWP) is (*XWP*)*xy* = *WXY t* (*ω*) [44].

#### 2.4.3. Wavelet Coherence (WTC)

WTC is to identify the possible relationship between two time series by finding the common points of two time series in different time scales and frequency ranges. It is a kind of enhanced time-frequency correlation. Recognition tool that can express the local coherence of the XWT in the time-frequency domain [45]. Even if they correspond to the low-energy region of the XWT, their correlation in the WTC spectrum may be significant. The expression of complex wavelet coherency (Complex Wavelet Coherency) *ρxy* is:

$$\rho\_{xy} = \frac{S\left(\mathcal{W}\_t^{XY}(\omega)\right)}{\left[s\left(\left|\mathcal{W}\_t^X(\omega)\right|\right)^2 S\left(\left|\mathcal{W}\_t^Y(\omega)\right|^2\right)\right]^{1/2}}\tag{11}$$

In the formula, *S* is a smooth spectrum operation in the time domain and the frequency domain. WTC is the absolute value of complex WTC. The mathematical expression is:

$$R\_{xy} = \frac{\left| S\left( W\_t^{XY}(\omega) \right) \right|}{\left[ s \left( \left| \left| \mathcal{W}\_t^X(\omega) \right| \right)^2 S \left( \left| \left| \mathcal{W}\_t^Y(\omega) \right|^2 \right) \right) \right]^{1/2}} \tag{12}$$

#### **3. Results**

#### *3.1. Accuracy Verification*

The selection of the best model parameters is the main factor that affects the simulation accuracy of the VIC model. There are seven main parameters that need to be input to the VIC model: the storage volume curve index B; the maximum base flow velocity Dm (m3/s); the nonlinear base flow velocity Ds (m3/s); the soil moisture content Ws (%); and the three layers of soil thickness d1, d2 and d3 (m) during the nonlinear base flow. The selection of model parameters was based on traditional parameter calibration methods [46]. The above parameters were adjusted according to the empirical interval of the hydrological parameters to make the peak discharge value match the measured data. The above steps were repeated several times until a satisfactory result was obtained. The parameters are calibrated by repeatedly calculating the NSE and R2 of the simulated and the measured values of runoff. The NSE and R2 of the monthly runoff of Wenquan hydrological station during the verification period (2002–2017) are 0.44 and 0.47, respectively. The Final selected parameters are described as follows (Table 1):

**Table 1.** Seven parameters of the VIC model.

#### 3.1.1. Accuracy Verification in Terms of Time

By correctly inputting various parameters and running the model, a dataset of daily values of the LST from 1960 to 2017 can be obtained. To ensure the accuracy of the simulation, the average monthly and yearly LSTs of 6 meteorological stations (Jinghe Station, Tuoli Station, Kelamayi Station, Wusu Station, and Alashankou Station) from 1960 to 2017 were calculated as the actual measured values for this time. The results were verified, the determination coefficient R2, the Root Mean Squared Error RMSE, the Mean Absolute Error MAE, the Nash-Sutcliffe efficiency coefficient NSE were selected to evaluate the simulation accuracy of the model. The R2, RMSE, MAE, NSE values of annual simulated LSTs are 0.51 2.21, 1.85, and 0.17, respectively (Figure 3b). The R2, RMSE, MAE, NSE values of monthly simulated LSTs are 0.98, 3.93, 3.31, and 0.93, respectively (Figure 3d). It can be seen from the respective boxplots (Figure 3e,f) that the difference between the simulated and measured monthly data is not very large. The medians are 11.11 ◦C and 11.65 ◦C respectively. The measured values are slightly higher. The simulated and measured annual values have a large gap. They are 8.56 ◦C and 10.07 ◦C, the difference is 1.5 ◦C, our simulated values are underestimated.

**Figure 3.** (**a**). Comparison of LST changes between measured and simulated annual values. (**b**). Verification accuracy of measured and simulated annual value LST. (**c**). Comparison of LST changes between measured and simulated monthly values. (**d**). Verification accuracy of measured and simulated monthly value LST. (**e**,**f**). Boxplot of monthly and annual simulated and measured value, respectively.

#### 3.1.2. Accuracy Verification in Terms of Space

In addition to verifying the LST over time, its spatial accuracy should be verified as well. Therefore, according to the remote sensing LST inversion dataset in China, the simulated LST based on the VIC model and its spatial distribution were compared and studied. However, through verification, it was found that the accuracy of the LST simulated by the VIC model for cold months is inaccurate, especially in January and December, where the LST is even greater than 0 in high-altitude mountainous areas in winter. According to Koch's research [30], who found that the warmth of the VIC model in cold months is obvious, VIC has a distinct seasonality in its bias with a warm bias in winter (~3 ◦C) as opposed to a slight cool bias in summer (~−1 ◦C), which is related to an underestimation of the actual evapotranspiration and altitude. Our study area is located in an arid area, and the actual evapotranspiration in winter is small, which does not affect the accuracy of the LST simulated by the VIC model for areas other than high-altitude mountainous areas. To ensure the reasonableness of the data, the simulated LST values of the high-altitude mountainous area in November, January and December were deleted and converted to a no-data area. Taking 2017 as an example, the spatial resolutions were 0.05◦ and 1/12◦ respectively, as shown in Figure 4. It can be seen that the LST simulated by the VIC model

has a high spatial consistency with the measured LST. The LST in the higher altitude area is significantly lower than that in other areas. As the time changes, the LST inside and outside the basin both "increase first and then decrease", and the simulated LST is slightly lower than the measured data.

**Figure 4.** Spatial comparison map of simulated and measured LSTs in the Ebinur Lake Watershed.

Of course, the accuracy needs to be improved, but for the data-scarce area, the simulation results can provide continuous and long-term serial data support, which has certain advantages.

#### *3.2. Temporal and Spatial Variation Characteristics of the LST*

The overall trend of the daily average LST in the study area was a unimodal distribution. The daily average LST (Figure 5, left) began to increase at the end of February, reached the maximum in mid-July, began to decline at the end of August, and reached the minimum at the end of January. The peak in mid-July (July 13) was 25.92 ◦C, and the LST dropped rapidly to the lowest value in late January (January 22) at −18.99 ◦C. The daily average LST of the entire study area was 6.95 ◦C.

The monthly average LST in the study area also showed a single-peak trend (Figure 5, middle). The peak monthly average LST appeared in July. The order of the July peaks in each decadal period was 2010–2017 > 80s > 90s > 2000–2009 > 70s > 60s, and the peak values were 26.06 ◦C, 25.60 ◦C, 25.58 ◦C, 25.41 ◦C, 25.07 ◦C and 24.89 ◦C, respectively.

A histogram was used to represent the distribution of the annual average LST from 1960 to 2017 (Figure 5, right). The figure shows that the average annual LST is on the rise, with a rising rate of approximately 0.0279 ◦C/a and an increase of approximately 1.62 ◦C in the past 58 years. The annual average LST is 6.40 ◦C. The annual average LST reached a trough in 1960, with a value of 4.54 ◦C; the figure clearly shows that the annual average LST after 2000 was significantly higher than that before 2000, which may be caused

by increasing human activities, the industrial and agricultural activities in the Ebinur Lake Watershed.

**Figure 5.** Temporal changes of daily (**a**), monthly (**b**) and annual average (**c**) LST in the Study area.

Based on the values of the daily LST simulated by the VIC model, the spatial distribution map of the inter-decadal and multi-year average LSTs from 1960 to 2017 was calculated. It can be seen that there are no obvious differences in the LST in the basin, but there are varying degrees of variation between different decadal periods.

The LST changes within the basin are not obvious, and the degree of change is different (Figure 6). The main manifestation is the gradual decrease of the high LST areas in the east and the gradual increase of the central LST. The LST range changed from −11.48–9.04 ◦C in the 1960s to −11.14–9.56 ◦C in the 1980s, and then to −3.5–10.51 ◦C in 2010–2017, with the average values being 3.28 ◦C, 5.24 ◦C, and 6.79 ◦C, respectively. These results show that the LST in this area has an increasing trend. Although the range of high-value areas in the east has decreased, gradually shrinking to the low vegetation coverage area near Lake Ebinur, the LST value has increased, and the multi-year average LST in other regions has increased compared with the values in the 1960s. Therefore, the LST of the basin is increasing.

**Figure 6.** Interdecadal changes and annual mean LST spatial distribution in the Ebinur Lake Watershed.

#### *3.3. Delayed Correlation Analysis of LST and Meteorological Elements*

The XWT and WTC of the dual time series were applied to the coefficients after continuous wavelet transform, and the XWT spectrum between the LST of the Ebinur Lake Watershed and other meteorological elements was obtained (Figure 7) and crossed with the WTC spectrum (Figure 8). The XWT and WTC can well reflect the "correlation" between changes in two different time series. The XWT generally reflects the strength of the "shared period" between sequences. The WTC generally reflects the consistency of the periodic "change trend" between sequences but does not directly reflect the intensity relationship of the change period. As shown in Figures 7 and 8, the meteorological elements and LST of the basin are bounded by 1989. The two figures together show the following:

**Figure 7.** XWT analysis of mean wind speed and LST (**a**), XWT analysis of mean temperature and LST (**b**), XWT analysis of annual precipitation and LST (**c**), XWT analysis of relative humidity and LST (**d**), XWT analysis of duration of sunshine and LST (**e**). (Note: The horizontal axis of the wavelet period diagram corresponds to the time of the original time series, the vertical axis represents the period of change, and the color represents the intensity of the change period. In this figure, yellow represents the high intensity of the change cycle. The thin solid line is the Cone of Influence (COI), and only the energy spectrum within this solid line needs to be considered. In order to eliminate the interference of the boundary effect, the thick solid line represents the key value that passes the significance test with 95% confidence; the arrow indicates the relative phase difference: → indicates that the two changes in phase are consistent, ← indicates that the two changes in phase are opposite, ↑ indicates that the phase change of the meteorological element is 90◦ ahead of the LST phase, and ↓ indicates that the phase change of the meteorological element is 90◦ behind the LST phase [47,48]). (The same below).

**Figure 8.** WTC analysis of mean wind speed and LST (**a**), WTC analysis of mean temperature and LST (**b**), WTC analysis of annual precipitation and LST (**c**), WTC analysis of relative humidity and LST (**d**), WTC analysis of duration of sunshine and LST (**e**). (Note: The right axis generally reflects the consistency of periodic "change trends" between sequences, similar to the correlation coefficient. The larger the coefficient, the higher the correlation and the more consistent the change.).

In the power spectrum of XWT, there was a significant area with a period of 1~2a in the average wind speed and LST from 1962 to 1966, the correlation is approximately positive; there are 4 significant areas between the average temperature and the LST, which are approximately 1.5a (1964–1967), 1~3a (2003–2008), 5.5–7a (2005–2009) and 7~9a (1984–1990). The two showed a positive correlation between the average temperature and the LST; the annual precipitation and the LST had significant areas of 0.5~2a and 6~8a in 1964–1972 and 1984–1994, respectively. Both have an approximately negative correlation; relative humidity and annual precipitation are the same, and there are significant areas of 1~2a and 6~9a in 1964–1971 and 1978–1994, respectively, and the two also show a similar negative correlation. Sunshine hours and LST only existed in a significant area of 1~2a from 1964 to 1970, and the two are in an approximately positive correlation.

In the coherence spectrum of WTC, the average wind speed and surface temperature have a significant area greater than 15 years from 1980 to 1996, and the average wind speed in this area is ahead of the surface temperature, that is, 1/4a ahead of time; the average temperature and the surface temperature in the area with the 95% confidence interval accounted for more than 60%, indicating that the average temperature is the main controlling factor of the surface temperature, and the two existed in 1963–1972, 1980–1996 and 2004–2010 on the time scale of 1~4a. For significant periods, most areas of the two are approximately positively correlated. Among them, there may be abrupt changes around 1980 and 1996 on the time scale of 1~2a; that is, there may be a "high-low" transition of surface temperature, and there are still two significant cycles of ~10a (1970–1994) and 5~7.5a (2000–2009); the annual precipitation and surface temperature only have a small cycle of 1~1.5a (1964–1966), and the two are negatively correlated; the other periodic feature is not obvious; the relative humidity and surface temperature have a significant area of 1~2.5a from 1964 to 1968 and a negative correlation, and the other two periodic features are not obvious; the sunshine hours and surface temperature have only one area where the periodic characteristics are not obvious.

#### *3.4. The Impact of LUCC on LST*

The statistics were classified according to the three-phase LUCC in this study. The specific changes in land cover types in the 1990s, 2001, and 2017 are shown in Figures 9 and 10, and the range distribution of the LST in the corresponding time period is shown in Figure 11.

The order of land cover types in the Ebinur Lake Watershed in the 1990s (Figure 9a) was open shrubland > grassland > bare ground > woodland grassland > crop land > water > woodland > evergreen needleleaf forest > mixed forest > urban and built-up > deciduous broadleaf forest. Among them (Figure 10), the largest area of open shrubland was 18,587.6 km2, accounting for 35.1% of the total area, followed by grassland, with an area of 18,476.3 km2, accounting for 34.9% of the total area. Bare ground accounted for 13.9% of the total area, and cropland accounted for 3.8% of the total area. Deciduous broadleaf forest occupied the smallest area, accounting for only 0.14% of the total area.

**Figure 10.** Changes in LUCC in the Ebinur Lake Watershed.

**Figure 11.** Box plots of LST.

The order of land cover types in the Ebinur Lake Watershed in 2001 (Figure 9b) was grassland > non-vegetation land > cropland > water > woodland > urban and built-up > permanent wetland > evergreen needleleaf forest > deciduous broadleaf forest > mixed forest > woodland grassland > open shrubland > deciduous needleleaf forest. Among them (Figure 10), the area of grassland has increased, reaching the largest area of 34,714.9 km2 and accounting for 65.6% of the total area, followed by non-vegetation land with an area of 11,904.9 km2, accounting for 22.5% of the total area. The two accounted for 88.1% of the total area. The area of cropland decreased to 4184.4 km2, accounting for only 7.9% of the total area. The area of open shrubland dropped sharply, most of them were transformed into grassland or non-vegetation land.

The order of land cover types in the Ebinur Lake Watershed in 2017 (Figure 9c) was grassland > non-vegetation land > cropland > water > woodland > urban and built-up > permanent wetland > evergreen needleleaf forest > deciduous broadleaf forest > mixed forest > woodland grassland > open shrubland > deciduous needleleaf forest. This is consistent with the ranking in 2001, and the area has not changed much. Grassland and non-vegetation land are also the main land cover types in the study area (Figure 10).

In the 90s, the LST of grassland ranged from 5.2 to 8.4 ◦C, with a median of 6.3 ◦C; the LST of cropland ranged from 5.3 to 9.5 ◦C, with a median of 8.1 ◦C; the LST of nonvegetation land was in the range of 6.5–9.7 ◦C, with a median of 8.4 ◦C.

In 2001, the LST of grassland ranged from 5.7 to 8.5 ◦C, with a median of 6.6 ◦C; the LST of cropland ranged from 7.9 to 9.7 ◦C, with a median of 8.7 ◦C; the LST of non-vegetation land was in the range of 6.6–9.7 ◦C, with a median of 8.1 ◦C.

In 2017, the LST of grassland ranged from 5.5 to 8.0 ◦C, with a median of 6.2 ◦C; the LST of cropland ranged from 7.1 to 8.7 ◦C, with a median of 8.3 ◦C; the LST of non-vegetation land was in the range of 6.7–9.6 ◦C, with a median of 7.8 ◦C.

During the period, the area of open shrubland, most of which was converted to grassland, cropland, and non-vegetation land, declined sharply; the area of grassland increased most obviously, and the LST was lower. Although the area of cropland was small, the LST was higher; compared with land types with higher vegetation coverage such as grassland, the LST of non-vegetation land was higher, which may have been because a large amount of solar radiation was absorbed by vegetation, lowering its LST.

#### **4. Discussion**

#### *4.1. Harmful Effects Caused by Changes in LST in Arid Regions*

On the global scale, arid regions cover approximately one-third of the earth [49], and are one of the most sensitive regions to wind erosion phenomena [50]. In some recent research, a positive association was observed between the LST and wind erosion events [16,18]. The soil system has a synergistic relationship between various physical and chemical processes and environmental variables that affect its development [51]. Increasing the LST and decreasing the soil moisture, especially in the warm seasons over the arid regions, can be a factor in exacerbating soil erosion due to their effect on the

salt accumulation on the soil surface [16,52]. Meanwhile, Different LUCC trends have different effects on the LST [53]. Heat sources (including bare land, urban and built-up land, and cropland) increase the LST, while cold sources (woodland, grassland, and water) reduce the LST [54]. Cropland and bare ground have a great impact on the regional LST, which is mainly reflected in the area and intensity of the impact [55]. This study found that between 1960 and 2017, the amount of open shrubs in the Ebinur Lake Watershed decreased sharply, and most of these areas were converted to grassland, cropland, and non-vegetation land. The vegetation coverage in the above-mentioned areas was low, and the LST was relatively high. This finding reasonably explains the variation in the LST of the Ebinur Lake Watershed from 1960 to 2017. However, the increased LST will lead to the reduction of regional soil moisture, vegetation degradation, and even increase the frequency of wind erosion and dust emission events. We believe this should attract the attention of local society.

#### *4.2. Evaluation of VIC Model for LST Simulation*

The LST is an important and complex hydrological state variable of land-atmosphere interactions. Land surface models (LSMs), such as the VIC model, are a key tool for enhancing the understanding of the process and providing predictions of the coupling between the terrestrial hydrosphere and its atmosphere. Distributed LSMs can be used to understand the hydrological status and flux of each grid cell, such as the LST. By comparing and analyzing the simulation accuracy of three different LSMs, Koch et al. [30] found that they all have certain spatial defects. The VIC has a distinct seasonality in its bias with a warm bias in winter (~3 ◦C) as opposed to a slight cool bias in summer (~−1 ◦C). The VIC model has a better simulation effect on the LST in warm months; in cold months, underestimation of the actual evapotranspiration will result in a high LST simulation. In this study, the same situation also appeared, but it may be because the study area is located in an arid area, that is, the actual evapotranspiration in autumn and winter, which are cold months, is inherently low, so except for high-altitude mountainous areas, the accuracy of the LST in cold months in plain areas is still good. However, the LST in high-altitude mountainous areas is obviously overestimated, only in November, December and January there is a big deviation in the mountainous area, and the LST in the areas close to the Tianshan Mountains is greater than zero in winter, which is unlikely to happen. This lack of spatial accuracy may be related to the occurrence of snow in mountainous areas, but further analysis is needed to study this in more detail.

#### **5. Conclusions**

This study used the VIC model to simulate the LST of the Ebinur Lake Watershed from 1960 to 2017, analyzed its temporal and spatial characteristics, conducted a multiscale crossanalysis of various meteorological elements and LST, analyzed the correlations between LUCC and LST, and obtained the following conclusions:

**Author Contributions:** N.A. conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents/materials/analysis tools, prepared figures and/or tables, authored or reviewed drafts of the paper, approved the final draft. J.D. conceived and designed the experiments, contributed reagents/materials/analysis tools, authored or reviewed drafts of the paper, approved the final draft. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by National Natural Science Foundation of China (no. 42171269); Xinjiang Academician Workstation Cooperative Research Project (no. 2020.B-001); National Natural Science Foundation of China (no. 41771470); National Natural Science Foundation of China (no. 4196105).

**Acknowledgments:** Sincere thanks are extended to the teachers and students of the Key Laboratory of Oasis Ecology Department of Xinjiang University for their support and help. Furthermore, the authors would like to thank the reviewers and editors for providing helpful suggestions to improve this manuscript.

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

#### **References**

## *Article* **Mapping Aquaculture Areas with Multi-Source Spectral and Texture Features: A Case Study in the Pearl River Basin (Guangdong), China**

**Yue Xu, Zhongwen Hu \*, Yinghui Zhang, Jingzhe Wang, Yumeng Yin and Guofeng Wu**

MNR Key Laboratory for Geo-Environmental Monitoring of Great Bay Area & Guangdong Key Laboratory of Urban Informatics & Shenzhen Key Laboratory of Spatial Smart Sensing and Services & School of Architecture and Urban Planning, Shenzhen University, 518060 Shenzhen, China<EMAIL_ADDRESS>(Y.X<EMAIL_ADDRESS>(Y.Z<EMAIL_ADDRESS>(J.W<EMAIL_ADDRESS>(Y.Y<EMAIL_ADDRESS>(G.W.)

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Aquaculture has grown rapidly in the field of food industry in recent years; however, it brought many environmental problems, such as water pollution and reclamations of lakes and coastal wetland areas. Thus, the evaluation and management of aquaculture industry are needed, in which accurate aquaculture mapping is an essential prerequisite. Due to the difference between inland and marine aquaculture areas and the difficulty in processing large amounts of remote sensing images, the accurate mapping of different aquaculture types is still challenging. In this study, a novel approach based on multi-source spectral and texture features was proposed to map simultaneously inland and marine aquaculture areas. Time series optical Sentinel-2 images were first employed to derive spectral indices for obtaining texture features. The backscattering and texture features derived from the synthetic aperture radar (SAR) images of Sentinel-1A were then used to distinguish aquaculture areas from other geographical entities. Finally, a supervised Random Forest classifier was applied for large scale aquaculture area mapping. To address the low efficiency in processing large amounts of remote sensing images, the proposed approach was implemented on the Google Earth Engine (GEE) platform. A case study in the Pearl River Basin (Guangdong Province) of China showed that the proposed approach obtained aquaculture map with an overall accuracy of 89.5%, and the implementation of proposed approach on GEE platform greatly improved the efficiency for large scale aquaculture area mapping. The derived aquaculture map may support decision-making services for the sustainable development of aquaculture areas and ecological protection in the study area, and the proposed approach holds great potential for mapping aquacultures on both national and global scales.

**Keywords:** multi-source remote sensing; aquaculture mapping; texture feature; Google Earth Engine; Pearl River Basin

#### **1. Introduction**

Aquaculture has become one of the fastest-growing food industries [1], and the fishery products of China play an important role in the international seafood market [2], with over 60% of the fish farmed in the world [3]. However, according to the 2018 State of World Fisheries and Aquaculture report (SOFIA 2018) by the Food and Agriculture Organization of the United Nations (FAO), the proportion of marine fish resources within a biologically sustainable level showed a downward trend in recent years [4]. In China, many lakes and coastal wetlands were reclaimed in the past few years in order to support the fast development of fisheries [5], putting tremendous pressures on environments and hampering regional sustainable developments [6]. Accurately mapping aquaculture areas is an important support to policy development and implementation at regional, national, and global levels, and to measure progress towards sustainable developments [7].

**Citation:** Xu, Y.; Hu, Z.; Zhang, Y.; Wang, J.; Yin, Y.; Wu, G. Mapping Aquaculture Areas with Multi-Source Spectral and Texture Features: A Case Study in the Pearl River Basin (Guangdong), China. *Remote Sens.* **2021**, *13*, 4320. https://doi.org/ 10.3390/rs13214320

Academic Editor: Qiusheng Wu

Received: 16 August 2021 Accepted: 25 October 2021 Published: 27 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Traditional field survey for aquaculture mapping suffers from low efficiency, and currently the satellite remote sensing technique is one of the most important methods due to its many advantages, such as low cost, wide monitoring range, high efficiency, and high repeated observations [8–10]. Optical and radar remote sensing images have been increasingly utilized to delineate aquaculture areas [11,12], and many methods have been developed for local [13], regional [14], and national scale [15] aquaculture mapping. Meanwhile, due to the periodic repeated observations of satellites, remote sensing images were used not only to map aquaculture areas at a single time [16], but also to map their time series distributions [17].

Previous studies can be roughly summarized into two categories according to basic mapping units: pixel-based and object-based approaches. The pixel-based methods are widely applied to the images with low and medium spatial resolutions. For example, artificially designed spectral and textural features were computed for each pixel, and a supervised machine learning classifier was used to map large scale aquaculture areas [18]; Sakamoto et al. [19] applied a wavelet-based filter for detecting inland-aquaculture areas from MODIS time series images, and deep-learning-based methods were used for aquaculture classification [20,21]. For the object-based classification method [22–24], images are segmented into many homogeneous segments, which are further classified through machine learning classifiers or classification rules derived from expert rules. For example, Wang et al. [25] segmented Landsat images into objects using multi-resolution segmentation method to extract raft-type aquaculture areas.

Considering data sources, most studies used only optical images for they are visually intuitive and easy to be understood. For mapping aquaculture facilities over small areas, high spatial resolution remote sensing images are frequently applied [26,27], whereas medium resolution images are generally used for mapping aquaculture facilities at a regional or national scale due to their wide coverage and better spectral resolution. For example, Ren et al. [28] combined Landsat series images and an object-based classification method to map the spatiotemporal distribution of aquaculture ponds in China's coastal zone. Synthetic aperture radar (SAR) images are also used for aquaculture mapping [14,29]. For examples, Hu et al. [29] detected floating raft aquaculture from SAR image using statistical region merging and contour feature; Ottinger et al. [14] employed time series Sentinel-1 images and object-based approach to map aquaculture ponds over river basins; and Zhang et al. [30] mapped marine raft aquaculture areas using a deep learning approach by enhancing the contour and orientation features of Sentinel-1 images.

Most studies applied satellite images at a certain time to map aquaculture areas, which may affect the accurate mapping of aquaculture areas for their dynamics. Although the influence of weather could be suppressed by carefully selecting images, other factors may also affect the accurate mapping of aquaculture areas. For example, some paddy fields are still water-dominated at the early stage of farming, and inland aquaculture ponds are drained during harvest time. Previous studies demonstrated that more accurate aquaculture maps could be obtained by using time series SAR images [31], and time series optical images was also proved to be effective in improving mapping accuracy [32,33]. Therefore, time series images are an ideal and reliable data source for aquaculture mapping.

The Google Earth Engine (GEE) platform provides a series of free remote sensing images, many kinds of image processing algorithms, and high-performance computing capabilities, and it can process huge amount of time series remote sensing images over a large-scale area [34–36]. Therefore, GEE has been widely used for mapping wetlands [37] and agricultural lands [38], and it also has the potentials to map aquaculture areas [39–42]. For examples, Xia et al. [39] proposed a framework for extracting aquaculture ponds by integrating existing multi-source remote sensing images on the GEE platform; Duan et al. mapped aquaculture ponds over coastal area of China using Landsat-8 images and GEE platform [41] and further analyzed their dynamic changes from 1990 to 2020 [42]. Existing studies were mainly focused on designing artificial image features or training deep leaning models to map a specific type of aquaculture area, such as aquaculture ponds. However, many different aquaculture types are found over a large area, and how to simultaneously extract multiple types of aquaculture areas over a large area has not been well studied.

This study proposed a novel approach for mapping aquaculture areas with multiple types over large areas. With the Pearl River Basin (Guangdong) of China as a case study, time series Sentinel images were used as a data source to overcome the accidental factors of single-time observation. The spectral indices (including normalized difference vegetation index (NDVI), normalized water index (NDWI), and normalized built-up index (NDBI)) derived from Sentinel-2A multispectral images, the VV and VH polarized data of Setinel-1 SAR images, and their derived texture features were used to map aquaculture areas using machine learning algorithms implemented in Google Earth Engine. The proposed method holds great potentials in simultaneously mapping different types of aquaculture areas over a large area.

#### **2. Materials**

#### *2.1. Study Area*

In this study, the Pearl River Basin (Guangdong) refers the part of the basin that located in Guangdong Province. The study area (Figure 1) covers the Pearl River Basin (Guangdong) and its 40 km buffer seaward from coastal line. The Pearl River with 2320 km long is the third-longest river in China, and it covers a region of subtropical maritime monsoon climate. The Pearl River Basin has been one of the most economically dynamic regions of China since ancient times. Local sea reclamation and marine aquaculture started very early, and especially the urbanization in the Pearl River Basin (Guangdong) is very fast since the launch of China's Reform Programme in 1979. Currently, this specific region is one of the most developed regions in China; however, limited river networks cannot meet the demands of aquaculture developments [43]. Therefore, the expansions of aquaculture land are mainly concentrated in the ocean and coastal region. With years of developments, this region has become one of the biggest aquaculture bases with massive abundant inland aquaculture ponds and mariculture areas. The aquaculture, including fish farming, shrimp farming, oyster farming, mariculture, algaculture (such as seaweed farming) and the cultivation of ornamental fish, is well-developed in this region.

**Figure 1.** Geographical locations of the Pearl River Basin (Guangdong) and study area (annual median images of sentinel-2A (2020) after cloud removal).

#### *2.2. Data*

Sentinel series satellites are the important branches of the Copernicus Programme satellite constellation conducted by the European Space Agency (ESA) and an essential part of Global Monitoring for Environment and Security (GMES). Sentinel-1 satellite constellation is the first of this specific program, and it carries C-band synthetic-aperture radar (SAR) instrument with a temporal resolution of 6 days [44]. The Sentinel-2A/B satellites equipped with Multispectral Instrument (MSI) are capable of acquiring high spatial resolution (up

to 10 m) optical images in 13 spectral bands every 5 days [45]. The ESA and European Commission's policies makes Sentinel's data easily accessible [46]. More detailed technical specifications of both Sentinel-1 SAR and Sentinel-2A MSI data are available in previous studies [47,48] and will not be repeated here.

The standard Level-1 (ground range detected, GRD) data product of Sentinel-1 in Interferometric Wide Swath (IW) mode and the calibrated Bottom of Atmosphere (BOA) products of Sentienel-2 used in this study were offered by GEE. Considering the data availability, this study was only focused on 2020. One arc-minute global relief model of the Earth's surface that integrates land topography and ocean bathymetry (ETOPO1) dataset was also collected to extract ocean bathymetry [49]. The Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) that comprises World Vector Shorelines (WVS) [50] was applied to extract shorelines of the study area. All these datasets were re-projected to the uniform coordinate reference system (WGS\_1984\_UTM\_Zone\_50N) and resampled to a same ground spatial resolution of 10 × 10 m.

#### **3. Methods**

The approach proposed in this study for large-scale aquaculture mapping with multisource remotely sensed images and the GEE platform includes seven main steps (Figure 2): (1) data preprocessing, (2) feature extraction, (3) water surface extraction (including aquaculture land), (4) sample selection, (5) Random Forest classification in the GEE platform, (6) morphological post processing, and (7) accuracy assessment. These steps were repeated several times to find the optimal parameters.

**Figure 2.** Approach proposed in this study.

#### *3.1. Data Preprocessing*

A bitmask band with cloud mask information (QA60) for the Sentinel-2A image is provided on the GEE platform. Cloud and cloud shadow affecting observation quality were identified by the code based on the Function of Mask (CFMask) [51], and they were marked in the QA60 band. The clear-sky pixels were selected according to QA60, and the Sentinel-2A images in 2020 with less than 20% cloud and cloud shadow cover were used in this study.

#### *3.2. Feature Extraction for Each Image Pixel*

#### 3.2.1. Image Features of Aquaculture Areas

Only aquaculture ponds show regular shapes and textures on the Sentinel-2A multispectral images (Figure 3), because they are usually very large compared with mediumresolution image pixels. The spectral signatures of water ponds are also clearly presented in NDWI image, and texture features are clearly presented in all optical, spectral indices and SAR images.

**Figure 3.** Images and features of typical aquaculture areas: (**a**) Sentinel-2A multispectral image, (**b**) NDWI, (**c**) NDBI, (**d**) NDVI, (**e**) VV Polarization image, and (**f**) VH Polarization image.

Mariculture facilities do not show obvious signatures in Sentinel-2A multispectral optical images, because they are relatively small compared with medium-resolution image pixels, or some are submerged in water. However, these aquaculture types show significantly different characteristics in NDWI, NDVI and NDBI images, because normalized difference image features amplify the differences. More specifically, the textures are very clear in Sentinel-1 SAR images, because the backscatter signature of the radar is sensitive to the roughness and structure of water surface.

Based on the above analyses, we found that NDVI, NDWI, and NDBI could better highlight the differences between aquaculture areas and other water bodies and the textures derived from index images; VV and VH images could further enhance the differences.

3.2.2. Normalized Difference Spectral Indices

In this study, NDWI, NDVI [52] and NDBI [53] were used as spectral features.

Normalized Difference Water Index (NDWI) (Equation (1)): NDWI was proposed by McFeeters [54], based on the fact of strong spectral absorption of water body in the near-infrared (NIR) band, and it is a robust method for detecting aquaculture water surface and has proved effective in segmenting muddy aquaculture ponds [39]. Correspondingly, NDWI can be used to distinguish land and water bodies.

$$\text{NDWI} = (\rho\_{\text{green}} - \rho\_{\text{rir}}) / (\rho\_{\text{green}} + \rho\_{\text{rir}}) \tag{1}$$

where ρgreen is the reflectance of green band (B3 of MSI), and ρnir represents the reflectance of near-infrared band (B8 of MSI).

In this study, the images with annual cloud cover less than 20% in 2020 were used to calculate a series of NDWI values. The NDWI series were sorted from small to large, and the 25% maximum NDWI values (located in the 75–100% of the sequence) were selected to calculate the average NDWI image. This preprocessing takes full use of all the image features for the whole year and eliminates the seasonal differences of aquaculture areas.

Normalized difference vegetation index (NDVI) (Equation (2)): Generally, aquaculture ponds usually have higher nutrients than natural water bodies, and the algae concentration is also higher. The waterbodies with high algae concentration generally have high reflectance in NIR spectrum compared with natural waters (without aquacultures). Therefore, NDVI may separate aquaculture waters from natural waters at a certain degree.

$$\text{NDVI} = (\rho\_{\text{nir}} - \rho\_{\text{red}}) / (\rho\_{\text{nir}} + \rho\_{\text{red}}) \tag{2}$$

where ρred represents the reflectance of red band (B4 of MSI), and ρnir represents the reflectance of near-infrared band (B8 of MSI). Referring to the processing of NDWI, we also sorted and filtered the NDVI values (75–100%) to compose the final average NDVI image.

Normalized difference built-up index (NDBI) (Equation (3)): NDBI is originally designed to detect built-up and barren areas. In practice, the extracted water mask often contains shadow noise with low reflectivity, which may be caused by tall buildings over built-up areas. The surrounding pixels of shadows usually have high NDBI values, and thus it is able to distinguish building shadows by analyzing their surrounding pixels.

$$\text{NDBI} = (\rho\_{\text{swir1}} - \rho\_{\text{nir}}) / (\rho\_{\text{swir1}} + \rho\_{\text{nir}}) \tag{3}$$

where ρswir1 represents the reflectance of the shortwave infrared band (B10 of MSI), and ρnir represents the reflectance of the near-infrared band (B8 of MSI). To further reduce the potential effects of shadow noise, we also sorted and filtered the NDBI values (75–100%) to compose the final average NDBI image. The morphological operation is an effective way to examine the neighbors of a pixel. Since a higher NDBI indicates a higher possibility of a built-up pixel, an inflation operator was conducted on the NDBI image to eliminate lathy and small building shadows, and a corrosion operator was further adopted to restore the image.

#### 3.2.3. Backscatter Features

A single Sentinel-1 SAR image is usually contaminated by speckle noises, and previous studies demonstrated the effectiveness of SAR image enhancement using time series images [55]. Following Xia's work [39], we selected the mean values of the VV and VH time series images of 75–100% in ascending order at the pixel level to form the final VV and VH image over the whole study area.

#### 3.2.4. Texture Features

Image texture provides relevant characteristics of spatial structures. In this study, the aquaculture lands hold relative regular shapes and arrangements, which are reflected in the regular change of gray levels of images pixels. Considering the results of previous studies, we combined different kinds of textures to address the problem of insufficient accuracy and to reduce the interference of the surroundings.

Gray Level Co-occurrence Matrix (GLCM) [56] is a method to analyze texture characteristics, and it calculates the correlation between two gray levels to reflect the comprehensive information about directions, intervals, and the changes in magnitude and speed. Some studies suggested that better recognition performance could be obtained by using extreme values to calculate texture features in different directions than in any direction or average direction. NDWI, NDBI, and backscattering coefficient of VV/VH were employed as input images for calculating GLCM, with gray level of 16, 150 × 150 m sliding windows and four directions (Figure 4). The NDVI image was mainly used to distinguish aquaculture ponds and natural waters, and it was not applied to compute texture features in this study.

**Figure 4.** Direction measure of GLCM calculation in this study.

Numerous texture feature bands (*n* = 18) can be generated using GLCM method in GEE. Redundant textural measures are possibly complex and may reduce computing speed, which may influence the performance of classification. Following Yu's work [57], Angular Second Moment, Contrast, Correlation, Variance, Sum Average, Inverse Difference Moment, and Entropy were selected and computed from GLCM (Table 1), and their maximum and minimum values in four directions (Figure 4) were analyzed. Therefore, for each input image, a 14-dimension feature vector was obtained (7 measures × 2 values).

**Table 1.** Textural measures based on GLCM method in this study.

#### *3.3. Water Surface Extraction*

The pixels with aquaculture areas mainly show pure water or water-dominated spectral responses, and the probability of water surface can be simply reflected by NDWI. Inspired by [39], the water surface was first detected from NDWI image by that a pixel with a NDWI value higher than a threshold was classified as water. The aquaculture facilities in the sea, rivers, or lakes show different spectral signatures from pure water bodies. The sizes of aquaculture facilities are relatively smaller than medium spatial resolution image pixels and some aquaculture facilities are submerged in the water, and thus the pixels with aquaculture facilities are mixed pixels. The NDWI values of pixels with aquaculture facilities are usually lower than those of pure water pixels but higher than those of other objects, and thus the mixed objects can be classified as water surface using a threshold.

#### *3.4. Sample Selection*

The aquaculture areas in the Pearl River Basin (Guangdong) hold regular spatial morphological characteristics in remotely sensed images (Figure 5). In this study, the aquaculture areas were roughly divided into inland aquaculture ponds and marine aquaculture. Inland aquaculture ponds are mainly distributed in the junction zone between water body (river, lake and sea) and land. The inland aquaculture ponds are formed through the reclamation of coastal wetlands or inland lakes, and they are usually partitioned by embankments and have regular and compact shapes (Figure 5b,c). Inland aquaculture pond units are usually clustered together, showing regular texture structures.

**Figure 5.** Typical aquaculture samples over the study area in high resolution Google Earth image (**a**): costal aquaculture pond (**b**,**c**), visible raft aquaculture (**d**,**e**), and concentrated-distributed cage culture area (**f**).

Marine aquaculture areas are often located in bays and inshore seawaters, and they include raft aquaculture (Figure 5d,e) and cage aquaculture (Figure 5f). The raft aquaculture facilities in neritic zone are composed of aquatic bamboo raft (for floating) and subaquatic thick rope (for fixing aquatic products). Correspondingly, the raft aquaculture areas are characterized by dull grey stripes (Figure 5d,e). The cage culture areas consist of aquatic plastic frames and suspended net cages (Figure 5f). Comparing with water bodies, cage aquaculture facilities appear as brighter colors in remotely sensed imagery, and their distributions are more concentrated, with regular rectangles.

In this study, a total number of 8317 samples were selected for classifier training (*n* = 7317) (Table 2) and accuracy assessment (*n* = 1000).

**Table 2.** Classification system used in this study.

#### *3.5. Random Forest (RF) Classification on GEE Platform*

GEE provides a serial of classification algorithms, including Minimum Distance (MD) [58], Random Forest (RF) [59], Support Vector Machine (SVM) [60], Regression Trees (CART) [61] and Naive Bayes (NB) [62]. According to previous studies, SVM and RF have better performance than other classification algorithms [63–65]. SVM algorithm performs good with limited training samples, however, RF performs better as training sets get larger [66]. RF is conceptually similar to tree-based learners, and they share same advantages [67]. RF can exhibit higher accuracy and efficiency even facing the unfavorable conditions of noise interference. Moreover, RF classifier shows better performance than other classifiers in land cover mapping using long time series images [68,69]. Therefore, RF classifier was employed in this study.

In this study, the composed NDVI, NDWI and NDBI products were also employed as spectral attributes of each pixel, the GLCM features of NDWI, NDBI, VV and VH were used as texture attributes, and the ETOPO1 value was applied as a topographic attribute. A feature vector of 60 dimensions (3 spectral + 4 images × 14 measures/image texture + 1 topographic) was obtained. A RF classifier on GEE platform was trained according to existing reference [70], in which the number of trees (*ntree*) was set to 100, the maximum depth was unlimited, the minimum sample number of each tree node was 1, and the number of features for each tree was set as the square root of the variable number.

#### *3.6. Morphological Post Processing*

Based on the above process, the aquaculture areas in the Pearl River Basin (Guangdong) were obtained preliminarily. However, the texture difference between aquaculture land and streamway were not clearly detected, which may cause misclassification. Streamway is one of the equipped facilities of aquaculture ponds. The obvious morphological difference between open water and streamway can be observed, and thus the compactness of each extracted aquaculture pond was calculated for modifying misclassification with Equation (4), which reflects the geometric shape of the object. The threshold of compactness was set as 11.5 through a trial-and-error process.

$$\text{Compactness} = \frac{\text{W} \times \text{H}}{\text{Area}} \tag{4}$$

where W and H are the width and height of the minimum enclosing rectangle, respectively, and Area is the number of inner pixels of the object. The small holes (≤1000 pixels) within aquaculture patches were filled, and the isolated small patches (≤3 pixels) were removed. An illustration of post processes is presented in Figure 6.

#### *3.7. Accuracy Assessment*

A total of 1000 samples of aquacultures extracted from Google Earth images using visual interpretation were used as ground-truth references. The confusion matrix was first obtained, and then overall accuracy and Kappa coefficient were applied to assess the accuracy of the proposed approach.

**Figure 6.** Overview of the morphological post processes of extracted aquaculture land in this study: (**a**) raw classification results based on RF, (**b**) exclusion of misclassified river, and (**c**) final classification map with postprocessing.

#### *3.8. Optimal Parameter Selection*

Several key factors affect the final classification performance: feature combination, window size to compute textures and threshold to extract water surface. The steps described in Sections 3.1–3.7 were repeated, and the values of these key factors maximizing overall accuracy were selected.

Feature Combination: Based on different types of features, including normalized difference spectral indices, backscatter features and image textures derived from spectral indices and SAR images, many kinds of feature combinations can be obtained. To investigate the optimal feature combination, six groups of features (Table 3) were used to obtain six aquaculture maps, respectively. The threshold *T* = 0.06 and texture window size with 15 were used. The determinations of *T* and window size are described in Sections 4.2 and 4.3, respectively. The 1000 random selected samples were used to evaluate their performances.

**Table 3.** Classification performances using different feature combinations.

Threshold to Extract Water Surface: Since aquaculture areas are first classified as a part of a water body, the threshold separating water and non-water surfaces plays an important role in the whole process. In this study, water surface was extracted by thresholding an NDWI image, and only one parameter *T* was used. The spectral indices, textures from indices and SAR images, as well as ETOPO1, were used. The window size was set at 15 × 15 pixels. The thresholds between −0.2 and 0.2 were tested, with a step of 0.02. The accuracies were evaluated using 1000 samples.

Window Size to Compute GLCM Features: The window size of gray-level co-occurrence matrix can directly affect the detection of aquaculture in a complex water environment. According to [71], a texture window should cover at least the minimum detectable target. Although inland aquaculture ponds have regular arrangements and clear textures, the ponds can be properly separated by thresholding NDWI image. However, the detection of mariculture facilities mainly depends on image textures. Based on our visual interpretation from Google Earth, the size of minimum detectable mariculture block is about 50 to 150 m. Therefore, a window size of 15 × 15 pixels was preferred. However, the experiments using a window of 5 × 5 pixels, 15 × 15 pixels, and 30 × 30 pixels were also carried out to determine the optimal size.

#### **4. Results**

#### *4.1. Aquaculture Maps Obtained Using Multiple Features*

Six maps were obtained using different feature combinations; their accuracies were assessed with 1000 validation samples (Table 3), and four typical sites of the maps were presented in Figure 7. The performances using only spectral indices were poor (Figure 7A and Table 3), and many water pixels were wrongly classified as aquaculture pixels, for some aquaculture facilities are submerged in waters and their spectral signatures are not obliviously different from water. The combination of spectral indices and SAR image features improved the performance (Figure 7B), especially for the second typical site, for SAR images reflect the backscatter characteristics of aquaculture facilities, which is totally different from a pure water body. However, many misclassifications still existed near shorelines, and the overall accuracy and Kappa coefficient are not satisfactory. The use of spectral indices and index-derived texture features greatly improved the classification performance (Figure 7C and Table 3), and most aquaculture ponds and mariculture areas were correctly detected, and an overall accuracy of 88.2% was obtained.

The results of combination D (Figure 7D) show the effectiveness of using SAR images and their corresponding textures. Comparing the results on the second site with those from combination C, it can be easily found that mariculture areas were better detected. However, there are many commission errors along coastlines, resulting in a lower OA (86.9%). The results in Figure 7 show that the combination of E and F obtained better performances than the combination of A, B, C, and D. However, the results of combination E had less holes, and the overall accuracy (89.5%) and Kappa coefficient (0.719) were both the highest compared with other combinations. Therefore, the feature in combination E was applied to obtain the final map.

#### *4.2. Optimal Threshold for Water Body Segmentation*

A total of 21 aquaculture maps were obtained using thresholds from −0.2 to 0.2; their overall accuracies and Kappa coefficients are presented in Figure 8, and five maps over a typical area are presented in Figure 9.

The overall accuracy and Kappa coefficient were low with a very small threshold (*T* = −0.2), due to many misclassified land pixels. With a small threshold, many mixed pixels containing waters were segmented as water surface, and further classified as aquaculture areas. As shown in Figure 9b, many embankment pixels were wrongly segmented as water surface and further misclassified as aquaculture areas. As the threshold increased, the misclassifications were gradually overcome (Figure 9b–d), the overall accuracy gradually increased, and finally the overall accuracy achieved the best result with *T* = 0.06. As the threshold further increased, the accuracy decreased for many mixed pixels and even some pure water pixels were segmented as land, and thus many aquaculture ponds were missed (Figure 9f). Generally, the optimal threshold is a balance of the two misclassification types, and the threshold 0.06 was preferred in our study.

#### *4.3. Optimal Window Size for Texture Calculation*

Figure 10 shows the maximum sum averaged texture derived from VH and the extracted aquaculture with a window size of 5 × 5 pixels, 15 × 15 pixels, and 30 × 30 pixels, respectively. Some small misclassifications were observed along shorelines with the window size of 5 × 5 pixels. With 15 × 15 windows, the textures better delineated the existing aquaculture areas, resulting in a satisfactory result (Figure 10f). With larger window size (30 × 30 pixels), the texture feature was over-smoothed. Moreover, the computational load dramatically increased with the increase in window size. Therefore, the window size of 15 × 15 pixels was used in this study.

**Figure 7.** Aquaculture maps obtained using different feature combinations. The first-row present images and reference aquaculture areas at four typical sites, and the results of combinations (**A**–**F**) are presented in the second to seventh rows. The details of feature combinations are presented in Table 3.

**Figure 8.** Overall accuracies and Kappa coefficients using different thresholds.

**Figure 9.** Typical classification maps using different thresholds: (**a**) Sentinel-2A MSI image, (**b**) *T* = −0.2, (**c**) *T* = −0.1, (**d**) *T* = 0, (**e**) *T* = 0.1, and (**f**) *T* = 0.2.

#### *4.4. Aquaculture Map and Accuracy Assessment*

With optimal feature combination, threshold, and window size, the final aquaculture maps were obtained (Figure 11). The aquacultures in the Pearl River Basin (Guangdong) had a total area of 1445.91 km2 in 2020. The inland aquacultures were concentrated in the coast of bays, including Zhenhai Bay (Figure 11a), Yamen watercourse (Figure 11b), Pearl River Estuary (Figure 11g) and Niwanmen watercourse (Figure 11h). The area along the coast had the most concentrated aquacultures, which are mainly enclosed sea aquacultures. Moreover, a great deal of inland aquacultures were widely distributed in the inland river network.

**Figure 10.** Maximum sum average texture of VH using different window sizes and the aquaculture maps: (**a**) Sentinel-2A false color images, (**b**–**d**) image features with a window size of 5 × 5 pixels, 15 × 15 pixels and 30 × 30 pixels, and (**e**–**g**) aquaculture maps.

**Figure 11.** Aquaculture maps in the Pearl River Basin (Guangdong) in 2020: (**a**) Zhenhai Bay, (**b**)Yamen watercourse, (**c**) Shenzhen Bay, (**d**) Baisha Bay, (**e**) Baisha lake, (**f**) Histogram of aquaculture area in major administrative regions, (**g**) Pearl River Estuary, and (**h**) Niwanmen watercourse.

Marine aquacultures accounted for 22.8% of the total area of aquaculture in the Pearl River Basin (Guangdong), and they were concentrated in the bays along the coast. The raft cultures were mainly distributed in Yamen watercourse (Figure 11b) and Baisha Bay (Figure 11d), and cage culture in Shenzhen Bay (Figure 11d) and and Baisha lake (Figure 11e). The aquaculture areas of the major administrative regions in the Pearl River Basin (Guangdong) in 2020 are shown in Figure 11h. Jiangmen City had the largest area of aquaculture with about 443.56 km2, which accounts for about 30.7% of the total aquaculture area in the study area. The second was Zhuhai City, with an aquaculture area of 251.52 km2, accounting for 17.4% of the total aquaculture area.

The validation of aquaculture area map with the validation samples held an overall accuracy of 89.5% with a kappa coefficient of 0.776, a Producer's Accuracy of 82.82%, and a User's Accuracy of 89.47% (Table 4). Some commission and omissions errors for the proposed approach were still observed. Using a window of 150 × 150 m, some water bodies and tidal flats adjacent to aquaculture were misclassified as aquaculture because of their similar edge texture values. In addition, the aquacultures with high sediment content have similar spectral characteristics with non-aqueous bodies, and they were frequently misclassified.

**Table 4.** Confusion matrix of proposed approach using 1000 validation samples.

The high user's and overall accuracies indicated the satisfactory performance of the proposed approach. However, the producer's accuracy was slightly lower than user's accuracy, for some mixed pixels containing parts of ponds were classified as water bodies at the stage of NDWI thresholding. This problem could be overcome by using sub-pixel mapping or high spatial resolution images.

#### **5. Discussion**

#### *5.1. Comparison with Other Aquaculture Maps*

The proposed approach in this study was compared with the methods proposed by Duan et al. [42] and Xia et al. [39]. The corresponding reference images were manually labeled. The maps are presented in Figure 12, and the confusion matrixes are presented in Tables 5–7. Duan et al. [42] adopted a spectrum, and spatial and morphological features of 30 m resolution Landsat images to build a Random Forest classifier to implement an automatic extraction of large-scale aquacultures. Xia's et al. [39] extracted the water surface from Sentinel-2A images using multiple thresholds, described the water patches using geometric and spectral features, and finally applied a Random Forest classifier to classify inland aquaculture ponds.

**Table 5.** Confusion matrix of Duan's method in typical sites [42].

#### **Table 6.** Confusion matrix of Xia's method in typical sites [39].

**Table 7.** Confusion matrix of proposed approach in this study in typical sites.

The map provided by Duan et al. [42] presented a general distribution of aquaculture ponds; however, almost all the details were missed (Figure 12b,j), for it was obtained using Landsat TM images with a spatial resolution of 30 m. Moreover, they conducted morphological closing and erosion operations post process, which further smoothed the details.

Xia's method and the proposed approach in this study were applied to two typical regions, and their derived aquaculture maps were compared. Xia's [39] and the proposed approach in this study produced similar aquaculture maps, and most of their details are clearer than Duan's map (Figure 12c,d,k,l) for Sentinel images of 10 m spatial resolution were applied and the post processing did not eliminate the details. The proposed method produced similar results to Xia's method (Figure 12c,d,k,l); however, less embankment pixels were misclassified. Therefore, the User's Accuracy of aquaculture areas (92.40%) was higher than that of Xia's method (77.97%). The main problem of the proposed method is that some isolated ponds were missed, resulting in a lower Producer's Accuracy (Table 7). However, the overall accuracy and Kappa coefficient indicated that the proposed method performed better over the tested areas.

The proposed approach in this study achieved better performance than Duan's and Xia's methods, which might be explained by that: (1) the texture features from normalized difference spectral index images increased the distinguishability of aquaculture area from other objects; and (2) radar images were integrated with optical ones, as they are sensitive to the texture structure of aquaculture areas, resulting in further improvement for aquaculture mapping. More important, the proposed approach was originally designed to simultaneously map aquaculture ponds and mariculture areas, whereas Xia's and Duan's methods can only map aquaculture ponds, and thus the proposed method has a better generalizability.

#### *5.2. Impacts of Mixed-Pixels*

The aquaculture maps in the Pear River Basin (Guangdong) were obtained with medium resolution Sentinel-1 SAR and Sentinel-2A multispectral images, and the results showed the advantages of using medium resolution images for large-scale thematic mapping. However, the accuracy and generalization of the proposed approach might be affected by mixed pixels of medium resolution images. The aquaculture pixels are often mixed with embankment among ponds, which often results in ambiguous boundaries and some inevitable errors. Thus, some parts of aquaculture ponds were classified as land, and some embankments were classified as ponds. It is difficult to overcome such errors using medium spatial resolution images, and very high spatial resolution images may easily overcome this problem.

The step of extracting water surface mask in our proposed approach is also affected by mixed pixels. Our method relies on the assumption that pixels covering mariculture facilities are still water-dominated, because the sizes of these facilities are usually relatively smaller than a pixel. Thus, they were first classified as waterbodies, and further been detected at later stages. However, this assumption may be not appropriate when using high spatial resolution images, because many pixels covering mariculture facilities are no longer water-dominated mixed pixels. In such case, these pixels will be segmented into non-aquaculture, affecting the accuracy of final map. Object-based methods might be an optimal solution for high spatial resolution images to overcome this limitation.

The proposed approach was designed for large scale aquaculture mapping using medium spatial resolution images. Although the boundaries are not detected very accurately, the aquaculture map over a large area can be efficiently obtained. In particularly, the medium spatial resolution images acquired by many satellite sensors (such as Landsat TM series and Sentinel-2A series) for several decades provide great convenience to monitor the development of the aquaculture industry.

**Figure 12.** Aquaculture maps obtained using three methods. (**a**,**i**) present the Sentinel-2A true color images, and their corresponding reference maps are presented in (**e**,**m**). The results obtained using Duan's [42], Xia's [39], and the proposed approach are presented in (**b**–**d**) and (**j**–**l**), respectively. Their comparisons with the ground truth maps are presented in (**f**–**h**) and (**n**–**p**), respectively.

#### *5.3. Selection of Time Series Images*

In this study, time series images for a whole year were used for aquaculture mapping in order to eliminate the influence of accidental factors, such as water-dominated paddy fields and dry ponds during harvest period. Time series SAR images also are effective in suppressing speckle noise of radar image. Although there are many advantages, two problems should be noted.

First, the seasonal characteristics of study areas may affect the selection of time series images. For example, this study area is located in southern China, in which the temperature is usually high, and the water does not freeze; thus, the images with good quality acquired at any time can be used. However, in northern China, the pond water freezes in winter and its optical properties will change. Therefore, the time series images should be selected according to the specific seasonal characteristics of study areas.

Second, the assumption and basis of using time series images are that ground entities are not changed suddenly. If some aquaculture ponds are converted into agricultural land or built-up areas in winter, they may still be recognized as aquaculture areas. To improve the ability to respond to such sudden changes, it is necessary to shorten the time interval of time series images.

#### *5.4. Limits and Future Works*

An effective and efficient approach was proposed in this study for aquaculture area mapping over large areas; however, some issues should be further investigated. Firstly, some narrow rivers adjacent to aquaculture ponds were still misclassified. Thus, more accurate classifier and post-processing methods are still worth investigating, and the combination of post-processing and river vector boundaries may be a potential solution. Secondly, the aquaculture areas were roughly classified into two types in this study: aquaculture ponds and mariculture areas, and more specific aquaculture types are needed to be investigated. Finally, only one global threshold was used to extract waterbodies (containing aquaculture areas); it was not always optimal for different and complex image scenes, and thus a locally adaptive thresholding approach is a potential solution to improve the segmentation of the water surface. Only the aquaculture maps in 2020 were obtained in this study, and they are not sufficient for more potential applications. With historic earth observation images, it is necessary to analyze the long-time spatial–temporal changes of aquaculture areas and their impacts on economy and ecological systems and further to provide supports for the sustainable developments of the study area.

#### **6. Conclusions**

A novel approach was proposed in this study for simultaneously mapping multi-type aquaculture areas over large scale areas by combining spectral and texture features from optical (Sentinel-2A multispectral) and radar (Sentinel-1) images, and a case study in the Pear l River Basin (Guangdong) showed its efficiency. The main contribution of this work could be summarized as follows:

(1) We analyzed the spectral and textural features of aquaculture areas and demonstrated the effectiveness of fusing multiple image features for aquaculture mapping. We found that the use of textural features derived from the spectral indices can greatly improve the mapping accuracy and the use of textural features derived from SAR images can further improve mapping accuracy, as they are sensitive to marine aquaculture facilities.

(2) The proposed approach could generate a more accurate aquaculture map than previous studies. Moreover, the proposed approach was implemented on the GEE platform, and has great potential for national-scale and long-term aquaculture mapping.

**Author Contributions:** Conceptualization, Z.H.; methodology, Z.H., Y.X.; validation, Y.X., Y.Y.; formal analysis, Z.H., Y.X., J.W. and Y.Z.; investigation, Y.X. and Y.Y.; data curation, Y.X. and Y.Y.; writing—original draft preparation, Y.X. and Z.H.; writing—review and editing, Z.H., J.W., Y.Z. and G.W.; visualization, Y.X.; supervision, Z.H.; project administration, Z.H.; funding acquisition, Z.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was jointly supported by the National Natural Science Foundation of China (NSFC) (No. 41871227), the Natural Science Foundation of Guangdong Province (No. 2020A1515010678, 2020A1515111142) and the Basic Research Program of Shenzhen (No. JCYJ20190808122405692,<PHONE_NUMBER>2628001).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The codes and data are available from this linkage: https://github. com/XuBenBen15045812917/Muilti-aquaculture\_extraction.git (accessed on 25 October 2021).

**Acknowledgments:** We thank Xia at Eastern China Normal University for providing their code and Duan at Jiangsu Normal University for providing their aquaculture map for comparison.

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

#### **References**

#### MDPI

St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Remote Sensing* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/remotesensing

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8131-6