Remote Sensing Retrieval of Total Phosphorus in the Pearl River Channels Based on the GF-1 Remote Sensing Data

: Total phosphorus (TP) concentration is one of the indicators for surface water quality evaluation. In this study, an indirect algorithm was proposed to retrieve TP concentration. This algorithm retrieves the TP concentration in urban waters based on Gaofen-1 (GF-1) remote sensing data. The algorithm uses the correlation between remote-sensing reﬂectance, optically signiﬁcant constituents of water (chlorophyll, suspended sediment, and organic matter (excluding algae)), and TP to establish a retrieval model. First, the concentrations of optically active components are retrieved using a semi-analytical model. Second, the correlation between TP and optically active components is used to retrieve the TP concentration in waters. The GF-1 remote sensing data for 7 August 2015 were used to perform remote sensing retrieval of TP concentration in the Pearl River channels in Guangzhou, China. The results show that the TP concentration in most areas of the Front Channel, Western Channel, Guangzhou Channel, and the western part of the Back Channel was higher than 0.2 mg / L, while the TP concentration in the middle and eastern parts of the Back Channel was generally lower than 0.2 mg / L. The mean absolute percentage error of the retrieval is 24.18%. The experimental results show that the model is suitable for remote sensing retrieval of TP in urban waters in Guangzhou.


Introduction
Recently, pollution in coastal waters and inland waters has become increasingly severe. Nitrogen and phosphorus nutrients that enter waters through various channels are important reasons for eutrophication [1][2][3][4]. Conventional water quality sampling and monitoring methods are time-consuming and demand large efforts and costs. Only local water quality information about the monitored section can be obtained. Water quality status and changes in a large area cannot be accurately obtained. Remote sensing technology has the advantages of a broad monitoring range, fast speed, and low costs, which facilitate long-term continuous monitoring. Remote sensing technology can also be used to detect problems such as pollution sources and a pollutant diffusion status, which are hardly possible by using conventional monitoring methods. As a result, remote sensing technology is increasingly employed in water quality monitoring studies [5,6].
Nutrients are not optically active substances. They have weak absorption and scattering energy of light in natural waters. Current remote sensing retrieval studies primarily retrieve nutrients using direct and indirect methods based on empirical models [7]. In general, models are only applicable to waters in a specific study area. Direct methods can be applied to obtain the relationship between nutrients and spectral data by using techniques such as regression analysis and the differential spectrum technique. For example, Huang et al. (2015) employed two bands of the Geostationary Ocean Color Imager (GOCI) remote sensing data to construct a total phosphorus (TP) retrieval model, and discovered that the TP concentration has a certain correlation with chlorophyll and suspended sediment (SS) (R 2 > 0.52). Moses et al. (2014) established a regression equation between the TP concentration and the reflectance of the red band. The retrieval results showed that the predicted values were well correlated with the measured values (R 2 = 0.76). Li et al. (2017) used the reflectance of four bands of Landsat 8 OLI to construct an empirical retrieval model for total nitrogen and TP, respectively. Gao et al. (2015) used corresponding combinations of reflectance of four bands of HJ-1A CCD2 images to retrieve the TP in Chaohu Lake. Indirect methods generally establish retrieval models for nutrients based on the relationship between nutrients and chlorophyll, SS, water temperature, etc. For instance, Goes et al. (1999)  Currently, studies of remote sensing retrieval of TP primarily focus on coastal waters or large lakes and large reservoirs [4,[6][7][8][9][10][11][12][13]. Relatively few studies address remote sensing retrieval of TP in rivers. Most retrieval models are directly established based on the relationship between TP and water reflectance or the relationship between TP and SS and chlorophyll. Previous studies have indicated that the TP concentration in waters has a strong correlation with optically active substances, such as SS and chlorophyll [6,8,13]. However, for rivers that flow through cities in Guangdong Province, China, organic pollution in urban river channels is usually severe. Therefore, when the indirect method is used to retrieve the TP concentration in urban river channels, the effects of SS, chlorophyll, and organic pollutants should be comprehensively considered.
The purpose of this study is to construct an indirect algorithm suitable for TP concentration retrieval in urban rivers. First, based on the measured water leaving reflectance and water quality analysis results, a semi-analytical algorithm was used to retrieve the concentration of suspended sediment (SS), chlorophyll (Chl, in this study, represented by the chlorophyll-a, Chla), and organic matter (excluding algae, organic pollutants in water mainly refer to the complex organic materials dominated by oxygen-consuming organic matter in water, including dissolved organics and suspended organic particles, which are similar to the organic pollution characterized by Chemical Oxygen Demand (COD). Therefore, the permanganate index (COD Mn ) is used to characterize them in this study in urban river channels. Regression equations between TP concentration and other water quality parameters were then constructed, according to the laboratory test results of water quality. Second, a retrieval model of TP concentration was established. Last, Gaofen-1 (GF-1) remote sensing images were used to retrieve the TP in the Pearl River channels in Guangzhou, and the performance of the retrieval algorithm was evaluated.

Measured Water Leaving Reflectance Data
In the two days prior to the GF-1 satellite imaging (5-6 August 2015), simultaneous water spectral measurement, water quality sampling, and an analysis experiment were performed. In the week prior to the satellite imaging, the weather conditions in Guangzhou were satisfactory with a gentle breeze and no precipitation. The sampling sites of the experiment are shown in Figure 2. The water leaving reflectance curves were measured with a spectrometer (ASD FieldSpec3 spectrometer manufactured by the U.S. Company ASD) with an effective measurement spectral range of 350-2500 nm. Note that A1-A14 denotes the sampling sites on 5 August 2015, and B1-B7 denotes the sampling sites on 6 August 2015. At each sampling site, the measurement of the water spectrum by the Above-Water Method [14] was repeated multiple times when performing the water spectrum measurement to eliminate random errors.

Measured Water Leaving Reflectance Data
In the two days prior to the GF-1 satellite imaging (5-6 August 2015), simultaneous water spectral measurement, water quality sampling, and an analysis experiment were performed. In the week prior to the satellite imaging, the weather conditions in Guangzhou were satisfactory with a gentle breeze and no precipitation. The sampling sites of the experiment are shown in Figure 2. The water leaving reflectance curves were measured with a spectrometer (ASD FieldSpec3 spectrometer manufactured by the U.S. Company ASD) with an effective measurement spectral range of 350-2500 nm. Note that A1-A14 denotes the sampling sites on 5 August 2015, and B1-B7 denotes the sampling sites on 6 August 2015. At each sampling site, the measurement of the water spectrum by the Above-Water Method [14] was repeated multiple times when performing the water spectrum measurement to eliminate random errors.

GF-1 Remote Sensing Data
The remote sensing data used in this study is a multispectral remote sensing image of China's GF-1 satellite, with a spatial resolution of 16 meters and an imaging date of 7 August 2015. The GF-1 satellite was launched by the China Aerospace Science and Technology Corporation (CASC) in April 2013. The wide field of view (WFV) imaging system is one of the key instruments operating onboard the GF-1 satellite, which includes four integrated cameras with a 16-m spatial resolution [15]. The satellite parameters of GF-1 are shown in Table 1 [16]. The steps to obtain GF-1 remote sensing data of the Pearl River channels in Guangzhou are shown in Figure 3. First, the remote sensing data were preprocessed. By using Equation (1) [17], radiometric calibration of the GF-1 multispectral remote sensing data was performed to convert the satellite-based observation values into radiance.
where ( ) is the converted radiance, DN is the satellite-based observation value, gain is the calibration slope, and offset is the offset of the absolute calibration coefficient. An image-based automatic inversion algorithm of aerosol optical density proposed by Qin et al. (2015) is employed for atmospheric correction to obtain multispectral reflectance data [18]. Second, water extraction is performed. The Multilayer Perceptron (MLP) neural network model is used to extract the water bodies in Guangzhou from the reflectance data in the study area obtained after data preprocessing [19]. The water extraction results are shown in Figure 4.

GF-1 Remote Sensing Data
The remote sensing data used in this study is a multispectral remote sensing image of China's GF-1 satellite, with a spatial resolution of 16 m and an imaging date of 7 August 2015. The GF-1 satellite was launched by the China Aerospace Science and Technology Corporation (CASC) in April 2013. The wide field of view (WFV) imaging system is one of the key instruments operating onboard the GF-1 satellite, which includes four integrated cameras with a 16-m spatial resolution [15]. The satellite parameters of GF-1 are shown in Table 1 [16]. The steps to obtain GF-1 remote sensing data of the Pearl River channels in Guangzhou are shown in Figure 3. First, the remote sensing data were preprocessed. By using Equation (1) [17], radiometric calibration of the GF-1 multispectral remote sensing data was performed to convert the satellite-based observation values into radiance.
where L e (λ e ) is the converted radiance, DN is the satellite-based observation value, gain is the calibration slope, and offset is the offset of the absolute calibration coefficient. An image-based automatic inversion algorithm of aerosol optical density proposed by Qin et al. (2015) is employed for atmospheric correction to obtain multispectral reflectance data [18]. Second, water extraction is performed. The Multilayer Perceptron (MLP) neural network model is used to extract the water bodies in Guangzhou from the reflectance data in the study area obtained after data preprocessing [19]. The water extraction results are shown in Figure 4.

Water Quality Analysis Data
During the experiment of the Guangzhou section of the Pearl River on 5-6 August 2015, water samples were collected at each sampling site of the experiment after the measurement of the water leaving reflectance spectrum of the water. Five liters of the surface water was collected at each sampling point and stored in a plastic sampling bottle, which was sent to the laboratory within 24 hours for water quality detection and analysis. The detection methods of water quality parameters are shown in Table 2. (Note: The water quality test was performed by the China National Analytical Center, Guangzhou).

Water Quality Analysis Data
During the experiment of the Guangzhou section of the Pearl River on 5-6 August 2015, water samples were collected at each sampling site of the experiment after the measurement of the water leaving reflectance spectrum of the water. Five liters of the surface water was collected at each sampling point and stored in a plastic sampling bottle, which was sent to the laboratory within 24 hours for water quality detection and analysis. The detection methods of water quality parameters are shown in Table 2. (Note: The water quality test was performed by the China National Analytical Center, Guangzhou).

Water Quality Analysis Data
During the experiment of the Guangzhou section of the Pearl River on 5-6 August 2015, water samples were collected at each sampling site of the experiment after the measurement of the water leaving reflectance spectrum of the water. Five liters of the surface water was collected at each sampling point and stored in a plastic sampling bottle, which was sent to the laboratory within 24 hours for water quality detection and analysis. The detection methods of water quality parameters are shown in Table 2. (Note: The water quality test was performed by the China National Analytical Center, Guangzhou).

Optical Parameters of Water Components
To establish semi-analytical models for retrieving concentrations of Chl, SS, and organic matter components, the inherent optical properties (IOPs) of each component are required. In the study of ocean color remote sensing, the optically active components in water are usually divided into suspended particles, photosynthetic pigments, and colored dissolved organic matter (CDOM or Gelbstoff, yellow substances) [20][21][22][23]. The IOPs of pure water, SS, Chl, CDOM, and other components can be obtained from relevant research literature [24][25][26][27][28][29][30][31][32][33][34]. However, for the urban reach of the Pearl River channels in Guangzhou, the organic pollution is usually serious and the composition of organic pollutants is complex. Li et al. (2002) used Gas Chromatography-Mass Spectrometry (GC-MS) to determine the Haiyin section water of the Pearl River channel in Guangzhou [35]. The results show that there were many kinds of organic pollutants in the water with a large concentration. Sixty organic compounds were separated and identified, including alkanes, alcohols, ketones, aromatics, polycyclic aromatic hydrocarbons, phenols, acids, and amines, which were from different pollution sources. Zheng et al. (2016) analyzed the water samples of the Back Channel of the Pearl River in Guangzhou and found that the dissolved organic matter (DOM) was more enriched in carbohydrate, amide, carboxylic, and aliphatic compounds [36]. CDOM (historically referred to as Gelbstoff, humic matter, or yellow substances due to its high humic matter content [33,37,38]) is a part of DOM. For the water body of the Back Channel of the Pearl River channel in Guangzhou, analysis results of Shong et al. (2013) show that only about 10% of DOM is humic matter (CDOM) [39]. The composition of organic pollutants in urban river channels in Guangzhou is complex and cannot be simply characterized by the optical parameters of CDOM. Therefore, it is necessary to measure the absorption coefficient and scattering coefficient of organic matter in urban water bodies with complex organic pollutants such as the water body of the Pearl River in Guangzhou. Brando et al. (2009) proposed a semi-analytical inversion method to estimate the concentration of optically active components (chlorophyll, CDOM, and non-algal particles (NAP)) in water. In this study, considering the serious organic pollution in the water body of the Pearl River in Guangzhou and the low content of CDOM in the dissolved organic matter, a semi-analytical inversion method for estimating the concentration of chlorophyll, organic matter, and non-algal particles (in this study, represented by the suspended sediment, SS) in the water is constructed based on the improved Brando algorithm. In the algorithm, the optical parameters of Chl are obtained from the results of Brando [27,28], and the optical parameters of pure water and SS are obtained from the results of Deng [30] and He [31]. The optical parameters of Chl, pure water, and SS are shown in Figure 5. The absorption coefficient and scattering coefficient of organic matter were measured using the measurement method proposed by He et al. (2011). Water samples from heavily polluted waters in the urban area of Guangzhou were collected, and then the optical parameters of organic matter were measured.

TP Regression Analysis
Currently, the retrieval of TP in waters is primarily carried out by establishing a relationship between TP and remote-sensing reflectance or other water quality parameters (such as Chla and SS) [4,[6][7][8][9]12]. In this study, regression analysis of TP concentration in the Pearl River channel in Guangzhou with the water sample quality analysis results and the measured spectral data shows that TP in the Pearl River channels in Guangzhou is well correlated with other water quality parameters. The correlation coefficient between TP concentration and a combination of CODMn, Chl, and SS is the highest (R 2 = 0.9056, F value for F-test = 47.886). Therefore, the regression in Equation (2) is selected for TP retrieval.
where , , , and 0 represent the concentration of TP, Chl, SS, and CODMn, respectively, in water.

Semi-Analytical Model for Water Quality Retrieval
To retrieve TP concentration using Equation (2), the concentration of optically active components, such as Chl, SS, and DOM (organic pollutant composition, represented by the permanganate index, CODMn), must be obtained by the corresponding water quality retrieval models. Water quality retrieval models include empirical models, semi-analytical models, and analytical models [40]. Empirical models directly establish the relationship between water quality and remote sensing spectra, which cannot be directly used for water bodies with different water quality conditions. Analytical models are based on clear physical concepts. However, obtaining many parameters is difficult. Semianalytical models are based on radiative transfer theory and appropriately simplified for solution. Semianalytical models are highly adaptable and have been extensively investigated. In general, water reflectance can be expressed as a function of optical parameters, observation geometry, and bottom reflectance of the relevant components of water bodies [20].
where ( ) is the water-leaving reflectance, ( ) is the absorption coefficient, ( ) is the volume scattering function, ( ) is the bottom albedo, H is the bottom depth, and is subsurface solar zenith angle.
is the subsurface viewing angle from nadir, and φ is the viewing azimuth angle from the solar plane [20], respectively. Gordon constructed a semi-analytical model of water quality retrieval based on the IOPs of water quality components when studying ocean water color remote sensing [21][22][23]41].

TP Regression Analysis
Currently, the retrieval of TP in waters is primarily carried out by establishing a relationship between TP and remote-sensing reflectance or other water quality parameters (such as Chla and SS) [4,[6][7][8][9]12]. In this study, regression analysis of TP concentration in the Pearl River channel in Guangzhou with the water sample quality analysis results and the measured spectral data shows that TP in the Pearl River channels in Guangzhou is well correlated with other water quality parameters. The correlation coefficient between TP concentration and a combination of COD Mn , Chl, and SS is the highest (R 2 = 0.9056, F value for F-test = 47.886). Therefore, the regression in Equation (2) is selected for TP retrieval.
where C TP , C c , C s , and C 0 represent the concentration of TP, Chl, SS, and COD Mn , respectively, in water. coe f c , coe f s , coe f 0 , and coe f ε are the retrieval coefficients.

Semi-Analytical Model for Water Quality Retrieval
To retrieve TP concentration using Equation (2), the concentration of optically active components, such as Chl, SS, and DOM (organic pollutant composition, represented by the permanganate index, COD Mn ), must be obtained by the corresponding water quality retrieval models. Water quality retrieval models include empirical models, semi-analytical models, and analytical models [40]. Empirical models directly establish the relationship between water quality and remote sensing spectra, which cannot be directly used for water bodies with different water quality conditions. Analytical models are based on clear physical concepts. However, obtaining many parameters is difficult. Semi-analytical models are based on radiative transfer theory and appropriately simplified for solution. Semi-analytical models are highly adaptable and have been extensively investigated. In general, water reflectance can be expressed as a function of optical parameters, observation geometry, and bottom reflectance of the relevant components of water bodies [20].
where R rs (λ) is the water-leaving reflectance, α(λ) is the absorption coefficient, β(λ) is the volume scattering function, ρ(λ) is the bottom albedo, H is the bottom depth, and θ w is subsurface solar zenith angle. θ is the subsurface viewing angle from nadir, and ϕ is the viewing azimuth angle from the solar plane [20], respectively. Gordon constructed a semi-analytical model of water quality retrieval based on the IOPs of water quality components when studying ocean water color remote sensing [21][22][23]41].
where a w (λ) and b bw (λ) represent the absorption coefficient and backscattering coefficient of pure water, respectively, and a j (λ), b bj (λ), and C j represent the absorption coefficient, scattering coefficient, and concentration of the j th water quality component at unit concentration. Substituting Equations (5) and (6) into Equation (4) and disregarding the scattering effect of organic matter [21,28], the expression of the optical water-leaving reflectance R rs can be obtained.
where C comp is the concentration of the water components, The subscript comp can be W, C, S, and O represent pure water, Chla, SS, and COD Mn , respectively.

TP Retrieval Model
After the optical parameters of each water component are obtained, the remote sensing retrieval model of optical water-leaving reflectance can be solved. From Equation (7), in order to express conciseness, the frequency λ in Equation (7) and the subscript rs of the water leaving reflectance R rs are omitted, and the following equation can be established.
where R is the water leaving reflectance at the corresponding band, k is the extinction coefficient, The subscript comp can be W, C, S, and O, which represent pure water, Chl, SS, and organic matter (represented by the permanganate index, COD Mn ), respectively. The subscript i corresponds to remote sensing data at the i th band. For the GF-1 multispectral remote sensing data, four multispectral bands exist. Selecting three of the four bands (in this study, band2-band4, i.e., green, red, and near-infrared bands are selected) and, substituting relevant parameters into Equation (8), a set of equations for the concentration of three water quality components is established, namely, Chl, SS, and COD Mn . The concentration of each water quality component in the water can be solved.
Equation (9) can be rewritten as a linear matrix inversion method proposed by Hoge & Lyon [42,[45][46][47]. The Hoogenboom et al. (1998) study, which used the linear matrix inversion Remote Sens. 2020, 12, 1420 9 of 26 method for inland water study on the Dutch Lake Braassem, is referred to as the Matrix Inversion Method (MIM) [46,47]: To simplify the expression, let where R is a diagonal matrix of water leaving reflectance, and the main diagonal element R i is the water leaving reflectance of remote sensing data at the i th band. K con and F are 3 × 3 matrices, where each column of the K con represents the extinction coefficient (or absorption coefficient for organic matter) of a water component (i.e., Chl, SS, organic matter). The element at the i th row and j th column is the extinction coefficient (or absorption coefficient for organic matter) of the j th water quality component for the remote sensing data at the i th band. The coefficient matrix F is a diagonal matrix, and the main diagonal element in the i th line represents the inversion coefficient f i corresponding to the i th band of remote sensing data, which is calculated by Equation (7). b bw i , b bc i , and b bs i are the backscattering coefficients of the corresponding remote sensing data of the water components (i.e., pure water, Chl, SS) at the i th band. k wi is the extinction coefficient of the corresponding remote sensing data of pure water at the i th band. The subscripts, W, C, S, and O represent pure water, Chl, SS, and organic matter (represented by the permanganate index, COD Mn ), respectively. Taking into account the measurement errors of the relevant optical parameters and the influence caused by other factors, the comprehensive retrieval coefficient G and correction coefficient ε of water components obtained by regression analysis of the measured water spectral in the study area are introduced. The remote sensing retrieval model of the concentration of Chl, SS, and COD Mn is obtained.
Substituting Equation (22) into Equation (2), we obtain the retrieval model for the TP concentration C TP :

Optical Parameters Measurement of Organic Matter
The optical parameters of organic matter were measured using the method proposed by He [31]. The schematic diagram of the measuring device is shown in Figure 6. During the measurement, the light from a 65w constant light source passes through a set of lenses and filters to form a beam of parallel light that enters a certain height of the water sample to be measured in a high-transmittance glass container. After the parallel light passes through the bottom of the glass container, the downward scattered light is filtered by another group of lenses and filters, and the downward parallel light is irradiated onto the standard board. The radiance of the downward parallel light is measured by a spectrometer (ASD FieldSpec3 spectrometer manufactured by the U.S. Company ASD). The whole set of devices is installed in a closed multi-layer black cabinet, and the experiment is performed in a completely dark indoor environment, which makes the environment of each experiment almost unchanged and reduces the influence of external environmental factors.

Optical Parameters Measurement of Organic Matter
The optical parameters of organic matter were measured using the method proposed by He [31]. The schematic diagram of the measuring device is shown in Figure 6. During the measurement, the light from a 65w constant light source passes through a set of lenses and filters to form a beam of parallel light that enters a certain height of the water sample to be measured in a high-transmittance glass container. After the parallel light passes through the bottom of the glass container, the downward scattered light is filtered by another group of lenses and filters, and the downward parallel light is irradiated onto the standard board. The radiance of the downward parallel light is measured by a spectrometer (ASD FieldSpec3 spectrometer manufactured by the U.S. Company ASD). The whole set of devices is installed in a closed multi-layer black cabinet, and the experiment is performed in a completely dark indoor environment, which makes the environment of each experiment almost unchanged and reduces the influence of external environmental factors.
The radiance detected by the spectrometer shown in Figure 6 as follows: where ℎ is the path attenuation factor, i.e., the proportion of the irradiance of the incident light decreasing with the increase of the distance, is the irradiance of the incident light source, , are the transmittance of the water surface and the bottom of the glass container, respectively, is the reflectance of the standard board, − is the transmittance of the water layer to be measured, and is the optical thickness of the water layer:. = * (25) where is the thickness of the water layer and is the extinction coefficient of the measured water body. Equation (24) contains multiple unknown parameters. In order to eliminate systematic errors and avoid directly calculating and , the ratio method is used to solve the extinction coefficient. The concentration of high organic pollution water collected in the field is a certain value. During measurement, the radiance of water bodies of different thicknesses is obtained by changing the thickness of the water body, and then the extinction coefficient of the water body is obtained by comparing the radiance information of different thicknesses. The ratio of the radiance 1 and 2 measured by two water layers with different thicknesses 1 and 2 is as follows: The radiance detected by the spectrometer shown in Figure 6 as follows: where D path is the path attenuation factor, i.e., the proportion of the irradiance of the incident light decreasing with the increase of the distance, E is the irradiance of the incident light source, T w , T g are the transmittance of the water surface and the bottom of the glass container, respectively, R b is the reflectance of the standard board, e −τ is the transmittance of the water layer to be measured, and τ is the optical thickness of the water layer: where d is the thickness of the water layer and k is the extinction coefficient of the measured water body.
Remote Sens. 2020, 12, 1420 11 of 26 Equation (24) contains multiple unknown parameters. In order to eliminate systematic errors and avoid directly calculating T w and T g , the ratio method is used to solve the extinction coefficient. The concentration of high organic pollution water collected in the field is a certain value. During measurement, the radiance of water bodies of different thicknesses is obtained by changing the thickness of the water body, and then the extinction coefficient of the water body is obtained by comparing the radiance information of different thicknesses. The ratio of the radiance L 1 and L 2 measured by two water layers with different thicknesses d 1 and d 2 is as follows: Equation (26) is sorted as follows.
In Equation (27), L 1 and L 2 are directly measured by the spectrometer while d 1 and d 2 can be accurately measured in the experiment. Thereby, the extinction coefficient k of the measured water body can be calculated.
For a typical organic polluted water body, the total extinction coefficient k total = α total + b btotal , which is measured in the laboratory, includes the contributions of various water components such as suspended sediment, chlorophyll, organic matter, and pure water. By substituting the total extinction coefficient k total into Equation (4), the total scattering coefficient b btotal can be obtained.
In Equation (28), the water leaving reflectance R rs of the water body was obtained through a simultaneous measurement by spectrometer during field sampling, and the total extinction coefficient k total of the water body was calculated by Equation (27) after conducting a measurement of the water sample collected from the field in the laboratory. Then, the total absorption coefficient α total of the water sample is shown below.
After the concentration of each component C comp of the water body (subscript comp: c, s, and o represents Chl, SS, and organic matter) obtained through water sample analysis, the absorption coefficient α o and backscatter coefficient b bo of the organic matter in the water sample can be obtained from Equations (5) and (6).

Model Verification and Error Analysis
The accuracy of the retrieval model can be evaluated by calculating the mean absolute percentage error (MAPE) using the following equation [48].
where n is the number of the in situ measured samples, C j and C j are the retrieved value and the measured value of the concentration of the relevant water quality component of the j th sample.

Water Quality Analysis Results
The water quality analysis results of the water samples collected at each sampling site in the Pearl River channel experiment are shown in Table 3. The analysis results show that, among the 21 sampling sites in the simultaneous measurement of the Pearl River channels in Guangzhou, compared with B3-B7 sampling sites in the east part of the back channel, the concentration values of the water quality components in the A1-A14, B1, and B2 sampling sites in the front channel, the west channel, and the west part of the back channel are much higher, among which the concentration differences of Chl and permanganate index (COD Mn ) are the most clear. The concentrations of Chla and permanganate index of A1-A14, B1, and B2 sampling sites were higher than 30 µg/L and 3.5 mg/L, respectively, while those of B3-B7 were lower than 21 µg/L and 3.4 mg/L, respectively. For the typical organic polluted water samples collected from the Pearl River channels in Guangzhou, the absorption coefficient and backscatter coefficient of organic matter calculated by Equations (30) and (31) are shown in Figure 7 after using the water leaving reflectance measured during sampling. The water quality analysis results and the extinction coefficient were measured in the laboratory. The scattering coefficient of organic matter is very low compared with its absorption coefficient, and its absorption spectrum decays exponentially with the increase of the wavelength. The absorption of organic matter includes the absorption of dissolved organic matter and the absorption of suspended organic particles, which has a strong absorption effect on visible light. Therefore, the reflectance of heavily organic polluted water is very low, and the water appears black.
Equations (30) and (31) are shown in Figure 7 after using the water leaving reflectance measured during sampling. The water quality analysis results and the extinction coefficient were measured in the laboratory. The scattering coefficient of organic matter is very low compared with its absorption coefficient, and its absorption spectrum decays exponentially with the increase of the wavelength. The absorption of organic matter includes the absorption of dissolved organic matter and the absorption of suspended organic particles, which has a strong absorption effect on visible light. Therefore, the reflectance of heavily organic polluted water is very low, and the water appears black.  Figure 8 shows the water leaving reflectance curves of the waters at the sampling sites in the Pearl River channel experiment on 5-6 August 2015. The figure shows the measured water leaving reflectance spectra in a wavelength ranging from 400 to 900 nm at each sampling site because the direct sunlight is blocked by the clouds. No measured spectra were obtained for the sampling sites of A7 and A8. As shown in Table A1, by using the GF-1 satellite spectral response function [49], the measured water leaving reflectance spectra are integrated to obtain the water leaving reflectance values, which correspond to the four GF-1 spectral bands at each sampling site. As shown in  Figure 8 shows the water leaving reflectance curves of the waters at the sampling sites in the Pearl River channel experiment on 5-6 August 2015. The figure shows the measured water leaving reflectance spectra in a wavelength ranging from 400 to 900 nm at each sampling site because the direct sunlight is blocked by the clouds. No measured spectra were obtained for the sampling sites of A7 and A8. As shown in Table A1, by using the GF-1 satellite spectral response function [49], the measured water leaving reflectance spectra are integrated to obtain the water leaving reflectance values, which correspond to the four GF-1 spectral bands at each sampling site. As shown in Figure 8

TP Inversion Model
The regression of TP concentration in the Pearl River channel in Guangzhou with the water sample quality analysis results and the measured spectral data are shown in Tables 4 and 5. In the regression analysis of various combinations of water sample quality parameters, which are shown in Table 4, the correlation coefficient between TP concentration and a combination of CODMn, Chla, and SS is the highest (R 2 = 0.9055, F value for F-test = 47.886). Similarly, in the regression analysis of various combinations of measured spectral data shown in Table 5, the correlation coefficient between TP concentration and a combination of four bands is the highest (R 2 = 0.7596, F value for F-test = 11.06). The two regression equations with the highest correlation in Tables 4 and 5 are used to calculate the

TP Inversion Model
The regression of TP concentration in the Pearl River channel in Guangzhou with the water sample quality analysis results and the measured spectral data are shown in Tables 4 and 5. In the regression analysis of various combinations of water sample quality parameters, which are shown in Table 4, the correlation coefficient between TP concentration and a combination of COD Mn , Chla, and SS is the highest (R 2 = 0.9055, F value for F-test = 47.886). Similarly, in the regression analysis of various combinations of measured spectral data shown in Table 5, the correlation coefficient between TP concentration and a combination of four bands is the highest (R 2 = 0.7596, F value for F-test = 11.06). The two regression equations with the highest correlation in Tables 4 and 5 are used to calculate the TP concentration of 19 sampling sites in Table 3 (excluding A7 and A8). The results are shown in Figure 9. The MAPE of the combination of COD Mn , Chla, and SS is 8.77%, which is better than that of the combination of four bands (the MAPE is 17.8%).   Figure 10 shows the retrieved concentration from the in situ measured spectra, measured concentration, and relative retrieval errors of TP, CODMn, SS, and Chla at the 19 sampling sites. First, by using the GF-1 satellite spectral response function, the measured spectra are integrated to obtain the water leaving reflectance values, which correspond to the four GF-1 spectral bands at each sampling site (Table A1). Using Equations (22) and (23) Figure 10 shows the retrieved concentration from the in situ measured spectra, measured concentration, and relative retrieval errors of TP, COD Mn , SS, and Chla at the 19 sampling sites. First, by using the GF-1 satellite spectral response function, the measured spectra are integrated to obtain the water leaving reflectance values, which correspond to the four GF-1 spectral bands at each sampling site (Table A1). Using Equations (22) and (23), the concentration of the relevant water quality components is obtained, and the relative retrieval error at each sampling site and MAPE of the model are calculated. For the in situ measured data, the MAPEs of the retrieved concentrations TP, COD Mn , SS, and Chla relative to the measured values are 17.98%, 30.98%, 54.06%, and 51.2%, respectively.

Water Leaving Reflectance of Remote Sensing Data
The remote sensing retrieval process of TP in the Pearl River Channels is shown in Figure 3. By using Equation (1), radiometric calibration of the GF-1 multispectral remote sensing data was performed to convert the satellite-based observation values into radiance. Then, an image-based automatic inversion algorithm of aerosol optical density proposed by Qin et al. (2015) was used for atmospheric correction. First, partition the entire image to select the dark pixels of shady vegetation, and estimate the phase function and scattering ratio of aerosol in the red and blue channels based on the path radiance. Second, on the basis of Gilabert's algorithm [50], the consideration of diffuse reflection on the surface is added, and the simplified radiative transfer equation is used to calculate the aerosol optical density of shading vegetation and dense vegetation dark pixels. Lastly, Kriging interpolation is used to calculate the aerosol optical density of multiple dark image elements for the distribution of the entire scene image, and then an atmospheric correction is carried out. The remote sensing water leaving reflectance corresponding to each sampling site is shown in Figure 11 and Table A2. As an example, Figure 12 shows the water leaving reflectance curves of the A1 sampling site. According to the atmospheric correction results, the Qin's algorithm is similar to the Fast Lineof-Sight Atmospheric Analysis of Hypercubes (FLAASH) algorithm. In the blue and green bands, the water leaving reflectance of Qin's algorithm is slightly lower than the FLAASH algorithm, while, in

Water Leaving Reflectance of Remote Sensing Data
The remote sensing retrieval process of TP in the Pearl River Channels is shown in Figure 3. By using Equation (1), radiometric calibration of the GF-1 multispectral remote sensing data was performed to convert the satellite-based observation values into radiance. Then, an image-based automatic inversion algorithm of aerosol optical density proposed by Qin et al. (2015) was used for atmospheric correction. First, partition the entire image to select the dark pixels of shady vegetation, and estimate the phase function and scattering ratio of aerosol in the red and blue channels based on the path radiance. Second, on the basis of Gilabert's algorithm [50], the consideration of diffuse reflection on the surface is added, and the simplified radiative transfer equation is used to calculate the aerosol optical density of shading vegetation and dense vegetation dark pixels. Lastly, Kriging interpolation is used to calculate the aerosol optical density of multiple dark image elements for the distribution of the entire scene image, and then an atmospheric correction is carried out. The remote sensing water leaving reflectance corresponding to each sampling site is shown in Figure 11 and Table A2. As an example, Figure 12 shows the water leaving reflectance curves of the A1 sampling site.
According to the atmospheric correction results, the Qin's algorithm is similar to the Fast Line-of-Sight Atmospheric Analysis of Hypercubes (FLAASH) algorithm. In the blue and green bands, the water leaving reflectance of Qin's algorithm is slightly lower than the FLAASH algorithm, while, in the red and near-infrared bands, it is the opposite. The water leaving reflectance of the algorithm is higher than that of the FLAASH algorithm. Compared with the measured water leaving reflectance, the results of the two atmospheric correction algorithms are higher than the measured values. Except for the green band, the MAPE of the other bands exceeds 100%, especially in the near infrared band. The MAPE of the two algorithms are more than 500%.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 27 the red and near-infrared bands, it is the opposite. The water leaving reflectance of the algorithm is higher than that of the FLAASH algorithm. Compared with the measured water leaving reflectance, the results of the two atmosp e.

Remote Sensing Retrieval Result of TP in the Pearl River Channels
By substituting the remote sensing water leaving reflectance of the three bands of water pixels and the IOPs of each water quality component in Equations (22) and (23), the TP concentration of each water pixel is retrieved ( Figure 13). The details of five areas with high TP concentration indicated in Figure 13(a) are shown in Figure 14. the red and near-infrared bands, it is the opposite. The water leaving reflectance of the algorithm is higher than that of the FLAASH algorithm. Compared with the measured water leaving reflectance, the results of the two atmosp Figure 12. Water leaving reflectance of the A1 sampling site.

Remote Sensing Retrieval Result of TP in the Pearl River Channels
By substituting the remote sensing water leaving reflectance of the three bands of water pixels and the IOPs of each water quality component in Equations (22) and (23), the TP concentration of each water pixel is retrieved (Figure 13). The details of five areas with high TP concentration indicated in Figure 13(a) are shown in Figure 14.

Remote Sensing Retrieval Result of TP in the Pearl River Channels
By substituting the remote sensing water leaving reflectance of the three bands of water pixels and the IOPs of each water quality component in Equations (22) and (23), the TP concentration of each water pixel is retrieved ( Figure 13). The details of five areas with high TP concentration indicated in Figure 13a are shown in Figure 14.
Remote Sens. 2020, 12, x FOR PEER REVIEW 18 of 27 Figure 15 shows the retrieved concentration from the GF-1 remote sensing data, measured concentration, and relative retrieval errors of TP, CODMn, SS, and Chla at the 19 sampling sites. First, by using the geographical coordinates of each sampling site, the remote sensing water leaving reflectance of each band corresponding to each pixel of sampling sites is extracted from the GF-1 remote sensing image (Table A2). Second, using Equations (22) and (23), the concentration of the relevant water quality components is obtained, and the relative retrieval error at each sampling site and MAPE of the model are calculated. For the GF-1 remote sensing data, the MAPEs of the retrieved concentrations TP, CODMn, SS, and Chla are 24.18%, 31.45%, 42.21%, and 71.99%, respectively.    Figure 15 shows the retrieved concentration from the GF-1 remote sensing data, measured concentration, and relative retrieval errors of TP, COD Mn , SS, and Chla at the 19 sampling sites. First, by using the geographical coordinates of each sampling site, the remote sensing water leaving reflectance of each band corresponding to each pixel of sampling sites is extracted from the GF-1 remote sensing image (Table A2). Second, using Equations (22) and (23), the concentration of the relevant water quality components is obtained, and the relative retrieval error at each sampling site and MAPE of the   (e)). There are numerous industrial enterprises and a dense urban population around these areas. Industrial pollution, and urban sewage increase the TP concentration of these waters. There are also some vegetable planting areas around the Shuikou Channel and Shiliugang Stream. The phosphorus emission caused by fertilization may also be one of the reasons for the high concentration of TP in these two areas. The possibility of the increase of the TP concentration in the water body caused by such fertilization is based on the fact that, at present, in the vegetable planting in South China, almost no vegetable growers will carry out "soil testing and formula fertilization." Over fertilization is a common phenomenon. Moreover, in Guangzhou, the temperature is very high in August, so the vegetable planting must be frequently watered to ensure the normal growth of vegetables, but part of the phosphorus in the artificial fertilizer will also be leached and lost due to watering and enter the nearby water body. Due to the dilution effect of the Sanzhixiang channel with low TP concentration, the TP concentration of waters in the middle and rear part of the back channel is relatively low. The TP retrieval's statistical results show that the TP concentration of most waters (approximately 66%) of the Pearl River channels in Guangzhou ranges between 0.1 mg/L and 0.3 mg/L, and the TP concentration of approximately 32% of the waters is larger than 0.3 mg/L.   (Figure 14c1,c2,d1,d2,e1,e2). There are numerous industrial enterprises and a dense urban population around these areas. Industrial pollution, and urban sewage increase the TP concentration of these waters. There are also some vegetable planting areas around the Shuikou Channel and Shiliugang Stream. The phosphorus emission caused by fertilization may also be one of the reasons for the high concentration of TP in these two areas. The possibility of the increase of the TP concentration in the water body caused by such fertilization is based on the fact that, at present, in the vegetable planting in South China, almost no vegetable growers will carry out "soil testing and formula fertilization." Over fertilization is a common phenomenon. Moreover, in Guangzhou, the temperature is very high in August, so the vegetable planting must be frequently watered to ensure the normal growth of vegetables, but part of the phosphorus in the artificial fertilizer will also be leached and lost due to watering and enter the nearby water body. Due to the dilution effect of the Sanzhixiang channel with low TP concentration, the TP concentration of waters in the middle and rear part of the back channel is relatively low. The TP retrieval's statistical results show that the TP concentration of most waters (approximately 66%) of the Pearl River channels in Guangzhou ranges between 0.1 mg/L and 0.3 mg/L, and the TP concentration of approximately 32% of the waters is larger than 0.3 mg/L. Figure 13b shows the inversion results on 24 October 2015. The results show that the TP concentration in most water bodies of the Pearl River channels in Guangzhou at the time of imaging is greater than 0.4 mg/L. Because there are no synchronized measured values, the inversion accuracy evaluation corresponding to the water quality monitoring data of Guangzhou in October 2015 queried in the Guangzhou Environmental Protection Geographic Information System was used [51]. The MAPE was 22.34%.

Error Analysis of the Retrieval Results
The retrieval error analysis of the model is shown in Figure 16. The retrieval of TP is affected by a combined effect of Chla, SS, and COD Mn , but is insensitive to the retrieval accuracy of the individual components. As shown in Figure 16, when the inversion error of any water quality component in Chla, SS, and COD Mn increases to 60% or 70%, or the inversion error of any two components exceeds ± 35%, or the inversion error of three components exceeds ± 20%, the MAPE of TP will increase from 8.7% to 20%. Therefore, as shown in Figures 10 and 15, when the retrieval MAPE of Chla, SS, and COD Mn change by 20.79%, 11.85%, and 0.47%, respectively, the MAPE of TP increases by 6.38%. The relative errors of Chla, COD Mn , and TP retrieved from the remote sensing data were larger than those of the in situ measured data, among which the retrieved errors of Chla at B2-B7 sites are all more than 100%, which produces a sharp increase in the retrieval of the MAPE of Chla. When B2-B7 are disregarded, the MAPE of Chla based on the in situ measured data and the GF-1 remote sensing data decreases to 19.18% and 25.01%, respectively. It is found that the measured concentrations of Chla, SS, and COD Mn in B2-B7 sites are much lower than those in other sites in Table 3. However, the comprehensive inversion coefficient introduced in this study will overestimate the concentration of water components in the site with low measured concentration (Figures 10 and 15), which makes the relative retrieval errors of the site with low measured concentration become larger. How to improve the inversion accuracy of water with low concentration of Chla and SS needs further study.
Remote Sens. 2020, 12, x FOR PEER REVIEW 21 of 27 queried in the Guangzhou Environmental Protection Geographic Information System was used [51]. The MAPE was 22.34%.

Error Analysis of the Retrieval Results
The retrieval error analysis of the model is shown in Figure 16. The retrieval of TP is affected by a combined effect of Chla, SS, and CODMn, but is insensitive to the retrieval accuracy of the individual components. As shown in Figure 16, when the inversion error of any water quality component in Chla, SS, and CODMn increases to 60% or 70%, or the inversion error of any two components exceeds ± 35%, or the inversion error of three components exceeds ± 20%, the MAPE of TP will increase from 8.7% to 20%. Therefore, as shown in Figures 10 and 15, when the retrieval MAPE of Chla, SS, and CODMn change by 20.79%, 11.85%, and 0.47%, respectively, the MAPE of TP increases by 6.38%. The relative errors of Chla, CODMn, and TP retrieved from the remote sensing data were larger than those of the in situ measured data, among which the retrieved errors of Chla at B2-B7 sites are all more than 100%, which produces a sharp increase in the retrieval of the MAPE of Chla. When B2-B7 are disregarded, the MAPE of Chla based on the in situ measured data and the GF-1 remote sensing data decreases to 19.18% and 25.01%, respectively. It is found that the measured concentrations of Chla, SS, and CODMn in B2-B7 sites are much lower than those in other sites in Table 3. However, the comprehensive inversion coefficient introduced in this study will overestimate the concentration of water components in the site with low measured concentration (Figures 10 and 15), which makes the relative retrieval errors of the site with low measured concentration become larger. How to improve the inversion accuracy of water with low concentration of Chla and SS needs further study.  It can be seen from the figure that, on 24 October 2015, the TP concentration of the Pearl River channels in Guangzhou was mostly higher than 0.4 mg/L, while on 7 August 2015, the TP concentration was relatively low, with the concentration of most water bodies less than 0.4 mg/L. This is because August is the wet season of the Pearl River. The replenishment of a large number of clean water bodies in the main stream of the Pearl River, Xijiang River, and Beijiang River has diluted the TP concentration of the water body. However, October has entered the normal water period, and the replenishment of clean water bodies in the middle and upper reaches has been continuously reduced and the dilution effect has been weakened, so that the TP concentration in the water body of the Pearl River channels in Guangzhou remains at a high level. 0% 20% 40% 60% 80% 100% -90% -80% -70% -60% -50% -40% -30% -20% -10% 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%   Figure 13a,b show the results of TP inversion of the Pearl River channels in Guangzhou by GF-1 remote sensing image with imaging time of 7 August 2015 and 24 October 2015, respectively. It can be seen from the figure that, on 24 October 2015, the TP concentration of the Pearl River channels in Guangzhou was mostly higher than 0.4 mg/L, while on 7 August 2015, the TP concentration was relatively low, with the concentration of most water bodies less than 0.4 mg/L. This is because August is the wet season of the Pearl River. The replenishment of a large number of clean water bodies in the main stream of the Pearl River, Xijiang River, and Beijiang River has diluted the TP concentration of the water body. However, October has entered the normal water period, and the replenishment of clean water bodies in the middle and upper reaches has been continuously reduced and the dilution effect has been weakened, so that the TP concentration in the water body of the Pearl River channels in Guangzhou remains at a high level.

Applicability of the Model
As shown in Figure 17, compared with the TP inversion model established by directly using the reflectance of four bands, the inversion accuracy of the TP indirect inversion model proposed in this study is higher, and COD Mn , which is another water quality index for the surface water environmental quality assessment, can also be obtained by inversion. The TP indirect inversion model is an empirical inversion model based on the measured data of the water body in the Pearl River channels in Guangzhou. This kind of model is usually only applicable to specific water bodies. When applying it to other types of water bodies, the measured data should be used for model calibration to obtain appropriate inversion parameters. According to the analysis results of Zheng et al. (2016), the proportion of the dissolved organic matter (DOM) in the water bodies of the Pearl River channels in Guangzhou in different seasons is somewhat similar. The category of DOM in the water bodies does not change much, except that the concentration is low in the summer during the wet season. Therefore, it is feasible to use this model to retrieve the TP concentration in the water bodies of the Pearl River channels in Guangzhou in different seasons. The water quality analysis results of each sampling site used to build the model in this study show that the concentration of Chla has a high correlation with the concentration of COD Mn (R 2 = 0.75), but the concentration distribution is uneven, which accounts for 52.63% of Chla > 70 µg/L and COD Mn > 4.5 mg/L, 36.84% of Chla < 40 µg/L and COD Mn < 4 mg/L. The number of samples with medium concentration distribution is small, with only A2 and A6, which accounts for 10.53% (Figure 18). It can be seen that the inversion error of Chla and COD Mn for low-concentration water bodies in the model constructed from these samples will be higher, which, thereby, increases the TP inversion error (Figures 10 and 16). In the following work, attention should be paid to the reasonable selection of sampling sites and the study of constructing the TP inversion model, according to the concentration distribution of water quality components.
It can be seen from Figures 10 and 11 that the atmospheric correction algorithm used in this paper is similar to the FLAASH algorithm in atmospheric correction, and the corrected water leaving reflectance is very different from the measured value. Although Qin et al. applications of this algorithm have achieved good results in retrieving the aerosol optical thickness of Hong Kong, China using the multispectral remote sensing image of Ziyuan-3 [18], to be successfully applied to GF-1 remote sensing image, the selection of dark pixels and the inversion algorithm of aerosol optical thickness in urban areas need to be further optimized.
Remote Sens. 2020, 12, x FOR PEER REVIEW 22 of 27 As shown in Figure 17, compared with the TP inversion model established by directly using the reflectance of four bands, the inversion accuracy of the TP indirect inversion model proposed in this study is higher, and CODMn, which is another water quality index for the surface water environmental quality assessment, can also be obtained by inversion. The TP indirect inversion model is an empirical inversion model based on the measured data of the water body in the Pearl River channels in Guangzhou. This kind of model is usually only applicable to specific water bodies. When applying it to other types of water bodies, the measured data should be used for model calibration to obtain appropriate inversion parameters. According to the analysis results of Zheng et al. (2016), the proportion of the dissolved organic matter (DOM) in the water bodies of the Pearl River channels in Guangzhou in different seasons is somewhat similar. The category of DOM in the water bodies does not change much, except that the concentration is low in the summer during the wet season. Therefore, it is feasible to use this model to retrieve the TP concentration in the water bodies of the Pearl River channels in Guangzhou in different seasons. The water quality analysis results of each sampling site used to build the model in this study show that the concentration of Chla has a high correlation with the concentration of CODMn (R 2 = 0.75), but the concentration distribution is uneven, which accounts for 52.63% of Chla > 70μg/L and CODMn > 4.5mg/L, 36.84% of Chla < 40μg/L and CODMn < 4mg/L. The number of samples with medium concentration distribution is small, with only A2 and A6, which accounts for 10.53% ( Figure 18). It can be seen that the inversion error of Chla and CODMn for low-concentration water bodies in the model constructed from these samples will be higher, which, thereby, increases the TP inversion error (Figures 10 and 16). In the following wo

Conclusions
TP concentration is one of the indicators for surface water quality evaluation. In this study, an indirect retrieval algorithm for TP concentration was proposed. Based on the GF-1 remote sensing data, the TP concentration in urban waters was retrieved. The algorithm establishes a retrieval model that is based on the relationship among remote-sensing reflectance, optically significant constituents of water (chlorophyll, SS, organic matter), and TP. First, the concentration of optically active components was obtained by retrieval analysis using a semi-analytical model. Second, the TP concentration of waters is retrieved by using the correlation between TP and optically active components. The retrieval results of TP using GF-1 remote-sensing water leaving reflectance on 7 August 2015 illustrated the distribution of TP concentration in the Pearl River channels in Guangzhou on this day. The TP concentration of most waters in the Front Channel, the Western Channel, the Guanzhou Channel, and the western part of the Back Channel of the Pearl River was larger than 0.2 mg/L. The TP concentration of waters in the east part of the back channel was generally less than 0.2 mg/L. Using TP concentration as an evaluation index, most waters in Guangzhou's river channels are classified as Class III and Class IV water bodies, according to the Environmental Quality Standards for Surface Water (GB 3838-2002). These waters satisfy the standards for general industrial water and entertainment water, which does not directly contact human bodies. The retrieval MAPE of TP concentration is 24.18%. This model is suitable for the retrieval of TP concentration in Guangzhou's urban waters.