Retrieving the Lake Trophic Level Index with Landsat-8 Image by Atmospheric Parameter and RBF : A Case Study of Lakes in Wuhan , China

The importance of atmospheric correction is pronounced for retrieving physical parameters in aquatic systems. To improve the retrieval accuracy of trophic level index (TLI), we built eight models with 43 samples in Wuhan and proposed an improved method by taking atmospheric water vapor (AWV) information and Landsat-8 (L8) remote sensing image into the input layer of radical basis function (RBF) neural network. All image information taken in RBF have been radiometrically calibrated. Except model(a), image data used in the other seven models were not atmospherically corrected. The eight models have different inputs and the same output (TLI). The models are as follows: (1) model(a), the inputs are seven single bands; (2) model(c), besides seven single bands (b1, b2, b3, b4, b5, b6, b7), we added the AWV parameter k1 to the inputs; (3) model(c1), the inputs are AWV difference coefficient k2 and the seven bands; (4) model(c2), the input layers include seven single bands, k1 and k2; (5) model(b), seven band ratios (b3/b5, b1/b2, b3/b7, b2/b5, b2/b7, b3/b6, and b3/b4) were used as input parameters; (6) model(b1), the inputs are k1 and seven band ratios; (7) model(b2), the inputs are k2 and seven band ratios; (8) model(b3), the inputs are k1, k2, and seven band ratios. We estimated models with root mean squared error (RMSE), model(a) > model(b3) > model(b1) > model(c2) > model(c) > model(b) > model(c1) > model(b2). RMSE of the eight models are 12.762, 11.274, 10.577, 8.904, 8.361, 6.396, 5.389, and 5.104, respectively. Model b2 and c1 are two best models in these experiments, which confirms both the seven single bands and band ratios with k2 are superior to other models. Results also corroborate that most lakes in Wuhan urban area are in mesotrophic and light eutrophic states.


Introduction
Available water resources, including rivers, reservoirs, lakes, coastal waters, and oceans, are emerging as a limiting factor, not only in quantity, but also in quality, for human development and ecological stability.The decline in water quality has become a global issue of significant concern for hydrological cycles [1,2].Conventional in-situ sampling is expensive, and laboratory work is time-consuming.The sample sizes are also limited and usually cannot cover all types of water bodies over a range of temporal and spatial scales [3,4].Water quality retrieval by remote sensing can quickly and punctually reflect spatial distribution of the water quality of specific lakes or areas [5,6], which makes up for the shortage of conventional monitoring and saves a lot of manpower, materials, and financial resources.
The water environment in inland lakes is largely different from that of the open ocean, such as high concentration of suspended solids and chromatic dissolved organic carbon [7].Inland water usually has a smaller surface area and more complicated spectral features [8]; remote sensing of inland water with the ocean algorithm is only partially successful.Giardino et al. developed a software package that incorporated their bio-optical model-based tool to estimate water quality and bottom properties from remote sensing images, originally intended to retrieve optical and benthic properties for lakes, but also applicable in other aquatic environments [9].Brando et al. [10] presented an adaptive implementation of the linear matrix inversion method, which accounts for naturally occurring spatial or temporal variability in inherent optical properties (IOPs) and mass-specific IOPs.Doerffer and Schiller described the case-II coastal water algorithm based on the artificial neural network (ANN), which relates the bi-directional water leaving radiance reflectance to these IOPs [11].Three ANN processors available as plug-in modules for the Basic ERS (European remote sensing satellite) and ENVISAT (A) (Environmental Satellite) ATSR (Along-Track Scanning Radiometer) and MERIS (Medium Resolution Imaging Spectrometer) Toolbox (BEAM) were validated in Lake Taihu, China [12].ANN algorithms obtain good results in European lakes and ANN with local tuning is significant and necessary and should receive more attention [11,12].Radical basis function (RBF) neural network has a good generalization ability and superiority in avoiding the tedious calculation of backpropagation.It can also overcome problems where the calculation easily falls into the local minimum.The advantages mentioned above indicate that RBF is more suitable than the back propagation (BP) neural network for water quality inversion [13,14].
For water, 90% of the energy received by satellite sensors is affected by atmosphere, while only 10% is radiation energy with water information [15].Atmospheric correction is one of the most limiting factors for an accurate retrieval of water constituents from satellite remote sensing data [16].There have been many atmospheric correction methods developed by researchers since the middle of the 1980s.ATREM [17], for "Atmospheric REMoval program", and ACORN [18], for "Atmosphere CORrection Now", assumed a flat homogeneous ground hypothesis with a Lambertian surface.Hadjimitsis et al. assessed the temporal variations of water quality in inland water bodies using atmospheric corrected satellite remotely sensed image data [19].The bathymetry of the Sarca River was mapped using a WorldView-2 (WV-2) image, for which researchers evaluated the atmospheric compensation (AComp) product [20].Bernardo et al. evaluated the outputs from several atmospheric correction methods (Dark Object Subtraction (DOS), QUick Atmospheric Correction (QUAC), Fast Line of Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH), Atmospheric Correction for OLI 'lite' (ACOLITE), and Provisional Landsat-8 Surface Reflectance Algorithm (L8SR)) to estimate total suspended matter concentrations (TSM) in the Barra Bonita Hydroelectrical Reservoir and demonstrated that L8SR is the most suitable method with the lowest TSM retrieval errors [21].Ibrahim et al. presented and evaluated an improved Atmospheric Correction (AC) for hyperspectral sensors developed within NASA's SeaWiFS Data Analysis System software package [22].Yi Wang et al. refined the MODIS Dark Target algorithm by considering that water-leaving radiance at 2.1 µm to be negligible regardless of water turbidity [8].
All the methods can be classified basically into two categories, namely scene-based empirical method and the radiative transfer modeling approach [19,23].The scene-based methods mainly include the flat field correction method [24], the dark object subtraction [25,26], the empirical line method [27], the internal average reflectance method [28], the cloud shadow method [29], the log residuals calibration [30], and the quick atmospheric correction [31].Most of these algorithms usually eliminate the atmospheric influence from the image itself under the assumption of a uniform atmospheric level, thus avoiding the complex scattering terms in the atmospheric radiation transfer equation [32].The second category comprises atmospheric correction based on the theory of radiative transfer, which establishes the strict radiative transfer model about the radiation transmission process in the range of the visible to short-wave infrared band from the sun to the surface and then to the sensor, and it describes the condition of the radiation source, the atmosphere status, the relationship between the reflection characteristics of the ground objects, and the radiances recorded by the sensor.The atmospheric correction is then conducted on the image by using the measured, or standardized, atmospheric parameters [19,23].Some typical methods of this category are 6S (Second Simulation of the Satellite Signal in the Solar Spectrum), FLAASH (Fast Line-of-sight Atmospheric Analysis Spectral Hypercubes), ATCOR (ATmospheric CORection), HATCH (High-accuracy Atmospheric Correction for Hyperspectral Data), and BRDF (Bidirectional Reflectance Distribution Function) [23].The accuracy of the radiative transfer modelling approach is very high in the case of known atmospheric data, but the application of this atmospheric correction algorithm is often limited by the real-time measurement of atmospheric parameters at the corresponding time of satellite transit [26,33].
The procedures of the above methods were complex and need too much time due to the human intervention in setting many configuration parameters, which greatly limit the applicability for large numbers of inland waters [23,32].Gong et al. proposed that the spatial variation of the gas content and aerosols optical thickness is obvious for a given atmospheric model or aerosols type in a large-scale area [33].Li et al. also mentioned that although some scientists have developed some improved atmospheric correction algorithms, which can be applied to coastal and inland waters based on some other assumptions, these algorithms can only get good results in specific areas [34].Moreover, Palmer et al. highlighted the challenge that the continentality of the atmosphere over inland waters and their proximity to the land surface also introduce additional difficulties for atmospheric and adjacency correction procedures, and this further impacts the performance of in-water algorithms [3].The ideal atmospheric correction can be achieved only by remote sensing image information itself rather than using auxiliary information, which includes in situ data or other remote sensing data [35].
Landsat-8 was successfully launched on 11 February 2013, and deployed into orbit with two instruments on-board: (1) the Operational Land Imager (OLI), with nine spectral bands in the visual (VIS), near infrared (NIR), and the shortwave infrared (SWIR) spectral regions; and (2) the thermal infrared sensor (TIRS), with two thermal infrared spectral bands (band 10: 10.6-11.2µm; band 11: 11.5-12.5 µm; see Figure 1).The relative spectral response of the TIRS bands is presented in Figure 2. The two TIRS bands were selected to enable the atmospheric correction of the thermal data using a split-window algorithm on the basis of the different sensitivities of two thermal infrared bands to atmospheric response [36,37].We can get the water vapor content based on the split-window algorithm or atmospheric difference by the two TIRS bands, which have different sensitivities to the atmosphere.The objective of this work is the following.(1) Take majority of lakes in Wuhan as examples, compare the water quality retrieved based on different inputs: image data with the FLAASH atmospheric correction, band ratios, and introduce the atmospheric parameters into the input layer of RBF neural network.(2) Search for an improved method for inland water quality retrieval.

Study Area
Wuhan, the capital city of Hubei Province, is the confluence of Yangtze River and Hanshui River and known as the "River City".The total area of water is 2217.6 square kilometers (km 2 ), which accounts for 1/4 of the total city.Rivers and lakes are interlaced together, bringing in riverside and lakeside ecological environments in Wuhan (Wang, 2013).The lake area is surrounded by a dense population, developed industry, and agriculture, which has been affected by a series of human activities for years.With the development of the economy, many lakes in Wuhan have been polluted to different degrees, some of them even appear to be black and stinky.We put forward the lake pollution control and the lake wetland environment protection measures for the lake wetlands by water quality monitoring and water quality status analysis, which aim to provide a technical reference for lake development and management.The lakes in Wuhan urban area and its surroundings make the trail area in this study.Figure 1 shows the distribution of most lakes in the study area and 43 samples.

Images
Considering the matching of the remote sensing images and the in situ data, we chose the images of Wuhan area on March 1, 2016.The L8 image come from the United States Geological Survey (USGS), with a resolution of 30 m, which has been systematically corrected for radiation and geometry.The images are marked as L1T and corrected by terrain geometry.In this paper, the pre-processing includes median filtering with a window of 5*5, radiometric calibration, FLAASH atmospheric correction, and water extraction.The pre-processing mentioned above was executed by ENVI5.3.The 5*5 median filtering before water extraction is aimed at dismissing the effect of non-water pixels, such as boats, bridges, roads, and pavilions on lakes.Radiometric calibration is to convert the digital number (DN) of an image into physical quantities, such as radiance, reflectance, or brightness temperature.In the pre-processing, we convert the DN to radiance by ENVI5.3.We extracted water bodies with the normalized difference water index (NDWI).To minimize influence of mixed pixels on the shore, only the lakes with an area larger than 71,100 m 2 have been masked.

In-situ Data
Since there were many monitoring stations, water samples were collected for three days, and we obtained in-situ data from March 1 to 3, 2016.The parameters of water quality include PH value, secchi depth (SD), chlorophyll a (Chla), total phosphorus (TP), total nitrogen (TN), and chemical oxygen demand (CODMn).All laboratory analyses were done in Wuhan Environmental Monitoring Center.With no rain or heavy rain, weather in Wuhan from February 27 to March 3 was sunny or cloudy.Natural factors have not caused large changes in water quality.The in situ data from March 1 to 3 can reflect the water quality of March 1.

Images
Considering the matching of the remote sensing images and the in situ data, we chose the images of Wuhan area on 1 March 2016.The L8 image come from the United States Geological Survey (USGS), with a resolution of 30 m, which has been systematically corrected for radiation and geometry.The images are marked as L1T and corrected by terrain geometry.In this paper, the pre-processing includes median filtering with a window of 5*5, radiometric calibration, FLAASH atmospheric correction, and water extraction.The pre-processing mentioned above was executed by ENVI5.3.The 5*5 median filtering before water extraction is aimed at dismissing the effect of non-water pixels, such as boats, bridges, roads, and pavilions on lakes.Radiometric calibration is to convert the digital number (DN) of an image into physical quantities, such as radiance, reflectance, or brightness temperature.In the pre-processing, we convert the DN to radiance by ENVI5.3.We extracted water bodies with the normalized difference water index (NDWI).To minimize influence of mixed pixels on the shore, only the lakes with an area larger than 71,100 m 2 have been masked.

In-Situ Data
Since there were many monitoring stations, water samples were collected for three days, and we obtained in-situ data from 1 to 3 March 2016.The parameters of water quality include PH value, secchi depth (SD), chlorophyll a (Chla), total phosphorus (TP), total nitrogen (TN), and chemical oxygen demand (CODMn).All laboratory analyses were done in Wuhan Environmental Monitoring Center.With no rain or heavy rain, weather in Wuhan from February 27 to March 3 was sunny or cloudy.Natural factors have not caused large changes in water quality.The in situ data from March 1 to 3 can reflect the water quality of March 1.

Trophic Level Index
The water quality parameter was calculated and categorized according to the trophic level index (TLI) method [38].Equation (1) and Equations ( 3)-( 7) come from the Lakes (Reservoirs) Eutrophication Assessment Methods and Classification Technology Requirements, China National Environmental Monitoring Station [38].
In Equation ( 1), TLI(Σ) is the comprehensive nutritional status index.TLI(j) represents the jth composite indicator with the corresponding weight W j .With Chla as the reference parameter, the correlation weight calculation formula for the normalization of the jth parameter is as follows: where R j implies the correlation coefficient between Chla concentration and jth indicator, and m is the number of all indicators.R j values given in Table 1 are from the survey results of 26 major lakes in China [39].Parameters for the eutrophication assessment of lakes (reservoirs) are Chla, TP, TN, SD, and CODMn.The expressions for calculating indexes are as follows: TLI(TP) = 10(9.436+ 1.624lnTP) ( 4) TLI(SD) = 10(5.118− 1.94lnSD) ( 6) Equation ( 1) and Equations ( 3)-( 7) come from the Lakes (Reservoirs) Eutrophication Assessment Methods and Classification Technology Requirements, China National Environmental Monitoring Station.According to Equation (2) and Table 1, we can calculate the weight Wj of five parameters, Chla, TP, TN, SD, and CODMn.Then, taking the parameter weight Wj and the parameter TLI (Equation ( 3)-( 7)) into Equation ( 1), respectively, we can obtain Equation (8).The 43 measured TLI can be calculated according to Equation (8).We randomly choose 33 data for the RBF model training, leaving the remaining 10 as validation data.

Methods
Figure 2 shows the frame of this work; k1 is the atmospheric water vapor parameter retrieved by the split-window algorithm, k2 is the atmospheric water vapor (AWV) difference coefficient based on the two TIRS channels, which have different response to AWV.The specific algorithms are as follows.Firstly, we calculated the TLI by in-situ data and randomly chose 33 data for the RBF model training, leaving the remaining 10 as the validation data.Secondly, we get the reflectance of L8, from band1 to band7, based on the longitude and latitude of the 33 data.The reflectance comes from the L8 image, which has been both radiometrically calibrated and FLAASH atmosphere corrected.Thus, we get the input of RBF neural network, which include two parts: the station reflectance used to the input of model training, and the pre-processed image (band1-band7) as input of model forecasting.Similarly, we built seven other models after getting k1 and k2.The image information applied to RBF in the seven models is radiance after being radiometrically calibrated.It is remarkable that the image data of the seven models are only radiometrically calibrated, without FLAASH atmospheric correction.Lastly, the eight models were validated and estimated by the 10 data of group B. The input parameter selection and the eight models will be described in detail in Sections 3.2 and 4.1.
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 21 The 43 measured TLI can be calculated according to Equation (8).We randomly choose 33 data for the RBF model training, leaving the remaining 10 as validation data.

Methods
Figure 2 shows the frame of this work; k1 is the atmospheric water vapor parameter retrieved by the split-window algorithm, k2 is the atmospheric water vapor (AWV) difference coefficient based on the two TIRS channels, which have different response to AWV.The specific algorithms are as follows.Firstly, we calculated the TLI by in-situ data and randomly chose 33 data for the RBF model training, leaving the remaining 10 as the validation data.Secondly, we get the reflectance of L8, from band1 to band7, based on the longitude and latitude of the 33 data.The reflectance comes from the L8 image, which has been both radiometrically calibrated and FLAASH atmosphere corrected.Thus, we get the input of RBF neural network, which include two parts: the station reflectance used to the input of model training, and the pre-processed image (band1-band7) as input of model forecasting.Similarly, we built seven other models after getting k1 and k2.The image information applied to RBF in the seven models is radiance after being radiometrically calibrated.It is remarkable that the image data of the seven models are only radiometrically calibrated, without FLAASH atmospheric correction.Lastly, the eight models were validated and estimated by the 10 data of group B. The input parameter selection and the eight models will be described in detail in section 3.2 and 4.1.

Radial Basis Function Neural Network
Most inland lakes belong to case-II water bodies, whose optical properties are complex, and water quality components affect one another.Thus, the relationship between radiance value and water component cannot be accurately expressed by a simple linear relationship [40,41].The RBF network is a type of neural network with a local approximation.The number of neurons in the RBF network is greater than in the BP network.The learning speed of RBF is faster and its ability of approximation and pattern recognition function is better [42,43].Existing theories have proven that a network with deviation and with at least one S function hidden layer plus a linear output layer can approximate any rational functions.Therefore, it is feasible to use a neural network for the remote sensing inversion of water quality parameters [39,44].Experiments in this paper were based on the neural network toolbox provided by MATLAB R2012b.The newrb() function can automatically increase the hidden neurons of the network until the mean square error satisfies the required accuracy or until neurons reaches the maximum.
The following steps are required to retrieve TLI through the RBF neural network: 1.
The input data of training samples is obtained by sampling the preprocessed remote sensing image.

2.
Data normalization includes the training samples and the input parameters.The data were normalized to minimize the influence of shore mixed pixels on retrieved results.Combined with the method of water threshold value and water pixel value statistics by ENVI5.3, the input and output were normalized to [−1,1] according to Equation ( 9), which intends to make the sample at the same magnitude, accelerate the convergence speed, and improve the generalization ability of the network.
where P (i,j) represents the normalized value of jth pixel value b (i,j) in band i, min(b(i,:)) and max(b(i,:)) represents the maximum and minimum value in band i samples, respectively.

3.
The NN was constructed by the newrb() function in the Matlab toolbox.The pixel values came from the L8 satellite image on 1 March 2016, and the normalized TLI were taken as training samples, when goal=0.1, spread=1.2,and the RMSE of the network reached the minimum.4.
The TLI of lakes in the study area was retrieved.On basis of the constructed neural network model, the input layer parameters are band combinations with water information and the output layer is the TLI.Finally, the sim function was applied to simulation and TLI spatial distribution was obtained after the inverse normalization.

Selection of Input Layer Parameters
The atmosphere consists of gases with different components, tiny suspended particles with different properties and scales, cloud droplets, raindrops, ice crystals, and other objects.Atmosphere effect on electromagnetic waves can be summed up in two physical processes, namely, absorption and scattering.These two processes can weaken the electromagnetic waves.When the sensor receives the electromagnetic waves affected by the atmosphere, it is bound to have the information of the atmosphere [45].Adding proper atmospheric parameters to the input layer of RBF may improve the accuracy of water quality inversion.The atmospheric parameters in the study area were retrieved based on split-window algorithm and used together with the b1-b7 bands of L8 as input parameters of the neural network.
The key to achieve AWV by the split-window covariance-variance ratio algorithm is that two bands have different wavelengths, which means a difference in AWV absorption.The algorithm has been mainly applied to an atmospheric window with a wavelength of 11-12 µm, and the atmospheric window can be divided into two bands, one of which is more sensitive to water vapor (WV) than the other [45].Given that the wavelength of thermal infrared channels of L8 is very close to AVHRR-4 and AVHRR-5, or even narrower (Table 4), we attempt to apply the split-window algorithm to L8 thermal infrared channels and obtain the split-window WV algorithm.It is assumed that the atmosphere is homogeneous under a clear sky, so that the radiative transfer equation can be derived from two adjacent pixels with different surface temperatures [45,49].
Here, 11 and 12 respectively refer to the split-window channels of 11 µm and 12 µm; i and j are two adjacent pixels, and it is necessary that they have the same surface emissivity, therefore, the change of surface temperature can be measured.These conditions can usually be found over the ocean area and the surface emissivity of the ocean is very close to a single value, which indicates that the two adjacent pixels can be assumed to have the same surface emissivity.The surface emissivity of the inland water body is also close.Thus, the above hypothesis is tenable.Under the condition of different surface temperatures, R 12,11 can be calculated by assuming that atmosphere is homogeneous in a region covered by N pixels (20*20 pixels in this study is an area of 0.6 km*0.6 km).Thus, this effective ratio can be written as follows [45]: T 11 and T 12 refer to the brightness temperatures of 11 and 12 µm channels, T 11 and T 12 are the average brightness temperature of two channels in the selected area.The numerator and the denominator in formula, respectively, represent the covariance and variance of the brightness temperature of split-window channels.
In this spectral region, the atmospheric absorption is mainly attributed to the continuous absorption of WV.Radiative transmittance is mainly attributed to the absorption of atmospheric constituents with relatively independent and measurable temperatures.According to the total gas absorption and the band average absorption coefficient, we can obtain the following expression [49]: where k λ is the absorption coefficient when the wavelength is λ and the total gas absorption is Q.Kleespies and Jedlovec asserted that, if we assume that the total transmittance due to the other atmospheric gases is constant and that the mean value of surface emissivity at 11 and 12 µm are equal, the ratio of the transmittance in the split-window channels is related to the total WV content in the atmospheric column [48,49].
C 1 and C 2 are constants related to the absorption coefficients of the atmospheric constituents and are obtained by the regression analysis between radiosonde measurements and computed transmittance ratio with radiative transfer models.Their results do not show a good agreement when the variation of estimated WV is higher than 50% [45,52].
Band 10 and 11 of L8 satisfy the requirements of the split-window algorithm.According to Equation (11), we can obtain the following expression, which describes the covariance and variance ratio of two channel brightness temperatures.
where T 10 and T 11 refer to the brightness temperatures of band 10 and 11, respectively, and T 10 and T 11 are the average brightness temperatures of two channels in the study area.R 11,10 is linearly stand for the total atmospheric column WV.We obtain another formula according to Equation ( 13): C 1 [45].To minimize the atmosphere effect, the visible or near-infrared bands and thermal infrared bands of L8 are set to avoid WV absorption bands, carbon dioxide, and ozone as much as possible.These bands are still partly affected by atmospheric absorption, among which the influence of WV occupies an absolute proportion.Compared with that, the influence of carbon dioxide and ozone can be ignored, so the atmospheric influence parameters are mainly related to the total WV content.Therefore, the hypothesis is as follows: where B λ is the radiation value received by satellite, and B T is the surface radiation, ε represents water emissivity, MT is the radiation of black body with the same temperature [53].We get the following expression on basis of the Planck function: h is the Planck constant, 6.626 •10 −34 J•s, K is the Boltzmann's constant, 1.38 •10 −23 J/k. and c is the speed of light in the vacuum, 2.998•10 16 µm/s, T is the absolute temperature (K) [53].We obtained water temperature at 10 a.m. by calculating the average value of temperatures at 8 and 12 a.m., which are measured by the automatic monitoring stations of Wuhan Environmental Monitoring Center.Combined with the above theories, we obtained the k1 and R 11,10 of four sites (Table 5).Through the scatter plot of four sets, we obtained a linear fit between k1 and R 11,10 and concluded the following expression k1 = −0.391* R 11,10 + 1.346. ( The R square (R 2 ) of the linear model was 0.867.According to the model, the AWV parameter k1 at each pixel in the study area was obtained.With respect to the thermal infrared data, the split-window method was widely used to retrieve WV on the ratio of split-window channel (11 µm and 12 µm) brightness temperature differences by assuming that the atmosphere in the split-window channels is spatially invariant.Several studies have validated this method and their results affirmed that the method can obtain reliable WV [34,54].Liang and Goward et al. improved the variance ratio algorithm of Jedlovec and made a regression slope between AVHRR-4 and AVHRR-5 to estimate AWV and applied it to the ecological investigation in a cross-section of Oregon, USA [55].As a function of AWV, the regression slope can represent the changes in the ratio between brightness temperatures of two heat channels.Under dry conditions, the brightness temperatures of AVHRR4 and AVHRR5 channels are almost the same and the regression slope is 1 [54].When air humidity increases, the influence on AVHRR-5 is significant and differences between the two channels occur.It confirms that the two thermal infrared channels, AVHRR-4 and AVHRR-5, have different sensitivities to AWV.The difference of WV content can be expressed by the differences of brightness temperature between the two thermal infrared channels if the sensor works normally [55].Ren et al. developed a modified split-window covariance-variance ratio to retrieve AWV from Landsat-8 Thermal InfraRed Sensor (TIRS) data [40].The key of this method is the sensitivities of the two thermal infrared channels to AWV absorption are different, and the band with a center wavelength of 12 µm is highly sensitive.
The spectral curves of L8 bands 10, and 11, and AVHRR-4 and 5 channels were drawn by ENVI5.3 (Figure 3).Although two small wave valleys emerge in the spectral curves of band11, the thermal infrared channels of L8 and AVHRR still have great similarities.Table 4 exhibits the wavelength of L8 and AVHRR thermal infrared bands.The band 10 and 11 of L8 have the same wavelength range as AVHRR-4 and 5, and the first thermal infrared band of L8 is narrower.We can use L8 to roughly estimate AWV parameter.When the WV content in the atmosphere increases, band 11 (11.5-12.5 µm) is highly sensitive and the bright temperature difference in the two bands becomes large.
On basis of the different sensitivities of two thermal infrared bands to atmospheric responses, we can roughly estimate the WV differences in the atmosphere through two brightness temperature ratios between the thermal infrared channels and describe the AWV difference among the different lakes in the study area.
The ratio between the two thermal infrared bands is k2, k2 = b10 /b11 (20) k2 is the AWV difference coefficient; b10 and b11 refer to the brightness temperatures of the two thermal infrared channels of L8.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 21 thermal infrared channels of L8 and AVHRR still have great similarities.Table 4 exhibits the wavelength of L8 and AVHRR thermal infrared bands.The band 10 and 11 of L8 have the same wavelength range as AVHRR-4 and 5, and the first thermal infrared band of L8 is narrower.We can use L8 to roughly estimate AWV parameter.When the WV content in the atmosphere increases, band 11 (11.5 -12.5μm)is highly sensitive and the bright temperature difference in the two bands becomes large.On basis of the different sensitivities of two thermal infrared bands to atmospheric responses, we can roughly estimate the WV differences in the atmosphere through two brightness temperature ratios between the thermal infrared channels and describe the AWV difference among the different lakes in the study area.
The ratio between the two thermal infrared bands is k2, k2 is the AWV difference coefficient; b10 and b11 refer to the brightness temperatures of the two thermal infrared channels of L8.

Model comparison and evaluation
Figure 4 and Figure 5 show the scatter diagram of measured-retrieved TLI.In the two figures (Figure 4 and Figure 5), b1-b7 represent the bands from band1 to band7; "+" means "together with" Figures 4 and 5 show the scatter diagram of measured-retrieved TLI.In the two figures (Figures 4 and 5), b1-b7 represent the bands from band1 to band7; "+" means "together with" or "and".A is group A, the training sample.B is group B, the validation data.Because the best retrieved TLI should be equal to the measured TLI, here we use function y = x to help describe the model performance."R2" and "RMSE" in figures represent the R2 and the RMSE of validations.We evaluated the eight models by RMSE (Table 6).The seven models built in this experiment have the same number of neurons.However, the eight experiments contain different input parameters, which are aimed at comparing the different results from the different input data.For instance, the input of model(a) are data with atmospheric correction.Model(a) is used for comparison with the other seven models, of which the input data have not been atmospherically corrected, but instead with atmosphere parameters or band ratios.Finally, we can conclude which input contains more reliable atmospheric information-the input with FLAASH atmospheric correction, band ratios, atmosphere parameters, or band ratios together with atmosphere parameters.
Figure 4 shows the different performance of four models.The four models mentioned in Figure 4 all have the seven single bands in input layers.For model(a), the scatter plots are concentrated and far away from the y = x line.They are mostly located below the line.Considering the TLI distribution of model(a) (Figure 6), the scatter distribution demonstrates that most predicted values are lower than measured.The training data is also far away from the y = x line, which also demonstrates that model(a) should be further improved.The scatter points of model(c) are dispersed and deviated from the line.The scattered points of model(c1) are mostly located on the line y = x, and the predicted value is close to the measured.For model(c2), most scatter points are far away the line y = x, with only a few on the line.From Figure 4, it is obvious that model(c1) is the best model among the four models.The RMSE is 5.389.Model(a), which only relies on seven single bands, has the worst results.Model (c), model(c1), and model(c2), which take atmospheric parameters into the inputs, are proved to be better than model(a).However, when adding k1, the atmospheric water vapor parameter, to models (model(c) and model(c2)), the RMSEs have not been greatly reduced.
Figure 5 shows the performance of four models based on band ratios and atmospheric parameters.Model(b) is only based on seven band ratios.A few points deviate from the line, which is below the line.The points in Figure 5 (model(b1)) are on both sides of the line and close to the line, but not concentrated.The points of model(b2) are also on two sides of the y = x line.Compared with model(b1), the distribution is more uniform.Most points of model(b3) are far away from the line.According to the RMSE, the error order of the four models is model(b3) > model(b1) > model(b) > model(b2).The comparison also indicates that the error would increase if we introduce the k1 into the band ratio model.
In total, the model(b2) and model(c1) are the two best models.The errors of the two are 5.104 and 5.389, respectively, much smaller than other models (Figure 6).The RMSEs and results of model(c1) and the model(b2) indicate that the ratio of L8 two TIR bands reflect the spatial distribution information of AWV.Once this difference is excavated by the RBF neural network, it could contribute to modifying the model affected by the atmosphere.Thus, model accuracy can be effectively improved.Having the biggest error, model(a), based only on the single bands with the FLAASH atmospheric correction, is not recommended for retrieving water quality.The result in Figure 6 (model(a)) cannot be used to describe the water quality of Wuhan in March 2016.When adding k1 to the single band models (model(c) and model(c2)), the errors are 8.361 and 8.904, a little smaller than the model(a), at 12.762 (Figure 4 and Table 6).However, in the band ratio models (model(b1) and model(b3)), the errors are 10.577 and 11.274, bigger than model(b), at 6.396.Adding atmospheric water vapor parameter k1 to band ratio models increase the errors.
from the line.The scattered points of model(c1) are mostly located on the line y= x, and the predicted value is close to the measured.For model(c2), most scatter points are far away the line y=x, with only a few on the line.From Figure 4, it is obvious that model(c1) is the best model among the four models.The RMSE is 5.389.Model(a), which only relies on seven single bands, has the worst results.Model (c), model(c1), and model(c2), which take atmospheric parameters into the inputs, are proved to be better than model(a).However, when adding k1, the atmospheric water vapor parameter, to models (model(c) and model(c2)), the RMSEs have not been greatly reduced.7 show the retrieved TLI of eight models.We display the results in segments and colors by ARCGIS10.2.Here, we describe the water quality in Wuhan urban area and surrounding lakes in March 2016 based on the best model(b2).The TLI ranged from 25 to 90, and most of them were between 40 and 60 (Figure 7, model(b2)), in the status of mesotrophic and light eutrophic.The TLI decreased from the lake center to the lake branches.Sewage waters, containing domestic sewage and city rain foul water, are drained into the lakes which are located in the city center, such as the East Lake, Shahu, South, and Moshui Lakes.With less amount of water exchange, accumulated pollution in lakes is serious and the endogenous pollution of urban water bodies appears.With the TLI ranging from 40 to 46, most water bodies of Liangzi lake are in a Mesotrophic state.
The shore mixed pixel effect may still exist.However, the experiment was based on pure water pixel extraction and images at the beginning of March 2016, a season where the lake shore grass had not yet flourished, hence the mixed pixel influence was small.Evidently, human activities, especially some unreasonable methods of sewage disposal, have directly aggravated lake water deterioration, and even resulted in black and odorous water bodies.The scope of remote sensing data includes some parts of Yangtze River and Hanjiang River, where water has relatively good liquidity and high or low sediment concentration in some sections.Therefore, the values when evaluating water quality with trophic level index method would appear too low or high, and results based on the method were unsuitable for Yangtze River and Hanjiang River.Moreover, the retrieved results only reflected the water quality of lakes in a low-water period and may change in a normal and high-water period.

Discussion
Due to the existence of atmospheric effects, errors caused by the seven single bands of L8 are unacceptable and the retrieved TLI also become unreliable.Therefore, the models which consider atmospheric influence would be more accurate.The band ratio can effectively reduce atmosphere influence on electromagnetic waves.If the band ratios, which have a relatively high correlation with water quality parameters, are used as the input terms of RBF neural network, the estimated TLI would be close to the measured value.

Discussion
Due to the existence of atmospheric effects, errors caused by the seven single bands of L8 are unacceptable and the retrieved TLI also become unreliable.Therefore, the models which consider atmospheric influence would be more accurate.The band ratio can effectively reduce atmosphere influence on electromagnetic waves.If the band ratios, which have a relatively high correlation with water quality parameters, are used as the input terms of RBF neural network, the estimated TLI would be close to the measured value.
According to the Figure 4, Figure 5, and Table 6, the model(c) established on the basis of the atmospheric parameter k1 retrieved by the split-window algorithm has improved the prediction accuracy, but the errors remain greater than the band ratio model(b).Moreover, when adding k1 into the input layers, we find the model accuracy has not improved obviously.In the band ratio models (Figure 5), model(b1) and model(b3) even have larger errors than the other two band ratio models (model(b) and model(b2)).These may result from the following factors.(1) Instead of real-time water temperature, the mean value of water temperature at 8 and 12 a.m. was considered the temperature at 10 a.m.The asynchrony of water temperature data and remote sensing data were bound to reduce model accuracy.(2) Using only four measured water temperature data to fit the relationship between k1 and R 11,10 was insufficient, and the accuracy of the AWV parameters retrieved by the linear fitting was also poor.(3) For inland waters, the emissivity between adjacent pixels is approximately the same.However, the atmosphere above the inland area is not as homogeneous as the ocean.For inland waters, the atmospheric water vapor parameter k1 retrieved based on the assumption of atmosphere being homogeneous may not be as reliable as for the oceans.(4) The pre-processing did not eliminate the noise that exists in the remote sensing data.The noise affected the retrieved accuracy of parameter k1 and was even input into the RBF neural network.The false atmospheric information increased the unreliability of inversion results.(5) Moreover, Xu proposed that the stray light outside the view field of TIRS has not been effectively solved, and that band11 of TIRS also has a calibration problem [56].Applying the split-window algorithm to the quantitative inversion of the AWV parameter for L8 remote sensing data is not mature.
It may be inadequate to choose lakes only in Wuhan, most of which are urban shallow lakes [57].Usually, we need several repetitive experiments which can reflect lake tropic states in different seasons to assure the effectiveness.However, the lakes in Wuhan urban area are not only numerous but also diverse in water quality, which include lakes with better water quality, such as Liangzi Lake, and also waters in a state of eutrophication, such as the Shahu, South Lake, and Tangxun lake.The differences in water quality contribute to the model training and the experiment generalization.Moreover, we chose numerous lakes in an area in our research instead of single one, which also increases the reliability of the experiment.Due to the time conflict between remote sensing data and in situ data, we need to accumulate more in situ data and seek for hyperspectral remote sensing data to build the comparison model in the time series and reveal the validity of this method.
TLI is a weighted sum based on the correlations between Chla and other substances [58].The chlorophyll-a concentration is not high because algae in the water had been dead or grew poorly in March in Wuhan.Whether the TLI (chla) retrieved in March is reliable or not remains to be seen from further research.The TLI calculated by this method may affect the accuracy of the retrieved results.However, TLI is a comprehensive index based on large numbers of sample data, which also include data obtained in spring.Thus, the index can be widely used in water quality evaluation of lakes (reservoirs).

Conclusions
The inversion results (model(b2) and model (c1)) show that waters of most lakes in Wuhan central urban area are mesotrophic and light eutrophic, and the TLI decreased from center to lake branches, which indicate the impact of human activities on water quality.With industrial, agricultural or domestic wastewater drained into the lake branches, lakes located in the city center, such as the East Lake, Shahu Lake, South Lake, and Moshui Lake, are containers for sewage.
The FLAASH atmospheric correction cannot eliminate the effect completely.Both models with atmospheric parameters and the band ratio model are better than the single-band model based on FLAASH atmospheric correction data.On the basis of the two thermal infrared channels with different sensitivities to AWV, we can regard the ratio of two thermal infrared channels as the AWV difference coefficient.The AWV difference coefficient model can accurately retrieve the TLI of urban swallow lakes and assess the nutritional state of water-quality.
For landsat8 remote sensing, the atmospheric parameters contribute greatly in improving the retrieval accuracy of inland water quality.Atmospheric water vapor difference coefficient k2 is superior to atmospheric water vapor parameter k1 in our experiment.However, neither adding k1 and k2 together with seven single bands or band ratios has achieved improved results.The error order of eight models is as follows.In term of RMSE, model(a) > model(b3) > model(b1) > model(c2) > model(c) > model(b) > model(c1) > model(b2).Among four models based on single bands, model(a) > model(c2) > model(c) > model(c1).Among the models based on band ratios, model(b3) > model(b1) > model(b) > model(b2).In term of these experiments, the incompleteness of FLAASH atmospheric correction directly leads to the error of model(a).Asynchrony of water temperature data and remote sensing data, insufficient fit in the relationship between k1 and R 11,10 by only four measured water temperature data, the hypothesis that atmosphere above the inland lakes is homogeneous, and the calibration problem of L8-TIRS reduces the accuracy of parameter K1.When applying K1 to the input of model(c), model(c2), model(b1), and model(b3), errors increase.
The neural network can search the abstract information of AWV spatial distribution.Thus, the noise effect of the atmosphere is reduced to a certain extent, and prediction accuracy is significantly improved.With optical complexity of inland waters, the RBF neural network, compared with linear models, can more accurately retrieve the water quality of lakes and provide a scientific reference for environmental protection departments.

Figure 1 .
Figure 1.Study area and stations.

Figure 1 .
Figure 1.Study area and stations.

Figure 2 .
Figure 2. The workflow of this research.

Figure 2 .
Figure 2. The workflow of this research.

Figure 3 .
Figure 3. Spectral of thermal infrared channels of Landsat8 and AVHRR made by ENVI5.3.

Figure 3 .
Figure 3. Spectral of thermal infrared channels of Landsat8 and AVHRR made by ENVI5.3.

Figure 4 .
Figure 4. Scatter diagram of measured-retrieved TLI based on seven single bands and atmospheric parameters.

Figure 4 . 21 Figure 5 .
Figure 4. Scatter diagram of measured-retrieved TLI based on seven single bands and atmospheric parameters.Remote Sens. 2018, 10, x FOR PEER REVIEW 13 of 21

Figure 5
Figure 5 shows the performance of four models based on band ratios and atmospheric parameters.Model(b) is only based on seven band ratios.A few points deviate from the line, which is below the line.The points in Figure 5 (model(b1)) are on both sides of the line and close to the line, but not concentrated.The points of model(b2) are also on two sides of the y=x line.Compared with model(b1), the distribution is more uniform.Most points of model(b3) are far away from the line.

Figure 5 .
Figure 5. Scatter diagram of measured-retrieved TLI based on seven band ratios and atmospheric parameters.

Figure 6 .
Figure 6.Distribution of TLI-retrieved data based on single bands and atmospheric parameters.Figure 6. Distribution of TLI-retrieved data based on single bands and atmospheric parameters.

Figure 6 .
Figure 6.Distribution of TLI-retrieved data based on single bands and atmospheric parameters.Figure 6. Distribution of TLI-retrieved data based on single bands and atmospheric parameters.

Figure 7 .
Figure 7. Distribution of TLI-retrieved data based on band ratios and atmospheric parameters.

Figure 7 .
Figure 7. Distribution of TLI-retrieved data based on band ratios and atmospheric parameters.

Table 1 .
The correlation between Partial parameters and Chla on Lakes (reservoirs) of China.

Table 2 .
Correlation coefficient between the radiance and the measured TLI.

Table 3 .
Input parameters of eight RBF neural network models.

Table 4 .
Two thermal infrared channels of Landsat8 and AVHRR.

Table 6 .
Validation errors of eight models.