Evaluation of the Crosta method for the retrieval of water quality parameters from remote sensing data in the Pearl River estuary

In recent years, global studies have shown that it is possible to pair in-situ data collected at discrete stations with remotely sensed data, with the aim of retrieving water quality parameters. Many algorithms have been developed using different kinds of methods (Miller & McKee 2004; Dall'Olmo et al. 2005; Gitelson et al. 2007, 2008; Gilerson et al. 2010; Shen et al. 2010; Xi & Zhang 2010; Mao et al. 2012; Ondrusek et al. 2012; Song et al. 2012; Wu et al. 2013; Chen et al. 2014; Knaeps et al. 2015; Shi et al. 2015; Wu et al. 2015). In recent years, machine learning-based regression algorithms have been applied in the retrieval of water quality parameters from satellite remote sensing images (Tang et al. 2013; Blix et al. 2018; DeLuca et al. 2018; Peterson et al. 2018; Ruescas et al. 2018; Stevenson et al. 2018). These algorithms ranged from linear relationships between satellite remote sensing reflectance and field measurements of water quality parameters (e.g., total suspended sediment concentration (TSM) and chlorophyll a (Chla)) in the estuary and inland waters to region-specific algorithms, explicitly developed for coastal waters, using in-situ measured reflectance and TSM or Chla. For example, Xi & Zhang (2010) established an empirical two-band model using the ratio of remote sensing reflectance at 629 and 671 nm to retrieve the TSM concentration in the Pearl River estuary (PRE). In addition, they used a MERIS image and in-situ remote sensing reflectance to map the distribution of the TSM in the PRE. Xing et al. (2013) used in-situ remote sensing reflectance and a combination of Hyperion bands using an exponential regression model to estimate the TSM concentration in the PRE, and they obtained a good performance. However, all these algorithms have no uniform model, because they lack a physical foundation, and most of them are therefore geographically dependent and cannot be applied to other areas (Zhang et al. 2010).

In the field of geology, the Crosta method has been used to extract remote sensing alteration information for mineral exploration. The idea of this method is to take advantage of the information contained in Landsat Thematic Mapper (TM) imagery for mineral exploration purposes. The Crosta method is essentially principal component analysis (PCA). An appropriate band combination can be used to conduct PCA, based on the specific spectral characteristics of different minerals in different spectral ranges, with the aim of selecting the principal component (PC) with specific mineral information. In 1989, the Crosta method was proposed and successfully applied to extract iron oxide and mudding anomalies from Landsat TM data (Crosta & Moore 1989). Each PC obtained by PCA using Landsat TM data often has a certain geological significance, which it does not share with any other component, meaning that each PC has unique characteristics. For example, the criterion for judging the PC of iron-stained minerals is based on the four bands of TM1, TM3, TM4 and TM5. The characteristic vector of the PC should be constituted by, and its TM3 coefficient should be the opposite of that of, TM1 and TM4. The coefficient symbol of TM3 is generally the same as that of TM5. Based on these criteria, the information on iron dyeing contained in the PC can be selected. Thus, this selected PC can be called the anomalous PC of iron dyeing. Crosta uses these spectral characteristics and combinations of several TM bands to diagnose remote sensing alteration information for iron exploration (Crosta & Moore 1989).

By the end of 2012, the Landsat TM instrument had obtained a large amount of data on the Earth's surface over 28 years (Wu et al. 2013). The image archive collected by the Landsat TM is historically unique, and it provides excellent opportunities for us to monitor and analyze the long-term spatiotemporal dynamics of the Earth's surface parameters, like TSM and Chla.

The objective of this study is to use Landsat TM imagery to evaluate the Crosta method in the retrieval of water quality parameters (e.g., TSM and Chla) in the PRE. Firstly, the spectral characteristics of TSM and Chla in the visible and near-infrared (NIR) band are described. Secondly, appropriate Landsat TM bands for diagnosing the TSM and Chla of the water components are selected. Thirdly, the effect of the combinations of eigenvectors with reflectance values of spectral bands on PC images is analyzed. Finally, an analysis of the correlation between the feature-orientated PC and corresponding MERIS products (TSM and Chla) is presented.

The Pearl River is the second largest river in China and is considered one of the most complicated fluvial networks in the world. It is comprised of three major tributaries (the western, northern and eastern rivers) and other small rivers, draining into the PRE, which occupies an area of ∼17,200 km². It plays a vital role in supplying fresh water to the large cities in the Pearl River Delta region, such as Macau, Hong Kong, Zhuhai, and Guangzhou. The annual suspended sediment load of the river is 88.7 Mt/y. The annual average flowrate is 2,281 × 10⁸ m³/y, and the annual suspended-sediment discharge is 6,567 × 10⁴ t/y for the western river. The annual average water and suspended sediment discharges for the northern river are 449 × 10⁸ m³/y and 864 × 10⁴ t/y, respectively, while those for the eastern river are 234 × 10⁸ m³/y and 236 × 10⁴ t/y, respectively (Wu et al. 2016). The total water and sediment discharges of these three major rivers account for more than 80% and 95% of the total load entering the sea, respectively. As a result, there has been a seaward extension of the mouth region of the PRE at a rate of 40 m per year (Wu et al. 2016). Soil erosion and soil weathering are very serious in the upstream of Pearl River. Sampling points were randomly selected in the PRE (Figure 1), which can be used for validation and statistical analysis.

Landsat TM scenes with no or minimal cloud cover (Row/Path: 44/122) were downloaded from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/). Standard MERIS FR level 2 products were obtained from the Europe Space Agency (ESA) website (https://merisFR-merci-ds.eo.esa.int/). Detailed information can be found in Table 1.

In this study, atmospheric correction was performed on Landsat TM images using the FLAASH atmospheric correction model with the ENVI 5.3.1 software. The pixel surface reflectance for all Landsat bands was retrieved. Due to the limitation of in-situ measurements, the results of the atmospheric correction were not analyzed and evaluated in this study. The MERIS sensor is designed mainly for ocean and coastal water remote sensing, with 15 narrow spectral bands in a range of 390–1,040 nm and a revisit period of one to three days. The eight bands centered at 412, 442, 490, 510, 560, 620, 665 and 708 nm were used, along with a neural network, to derive the MERIS level 2 products (Kratzer et al. 2008). MERIS FR level 2 products were projected into Geographic Lat/Log (WGS84). The sample points that matched with the feature-orientated PC were extracted using Beam 4.8 software (http://www.brockmann-consult.de/cms/web/beam/welcome). The average PC values that matched with the sample points were extracted using a 3 by 3 window on specific PC images. Detailed information on the bands and spectral range can be found in Table 2.

A single MERIS FR level 2 product has many bands, including TSM, Chla, yellow substance and other parameters. Thus, TSM and Chla images in the PRE could be extracted from MERIS FR level 2 products on Oct. 17, 2003 and Oct. 17, 2009 using the Beam 4.8 software. These datasets are assumed as the ground measurements and are further employed to validate the performance of the feature-orientated PCs.

Water components, TSM and Chla, have unique spectral characteristics in the electromagnetic spectrum, and water reflectance is highly variable over the visible and NIR spectral regions. The main source of TSM is suspended sediment, full of organic and inorganic particles (Xing et al. 2013). The spectral characteristic of TSM are very similar to those of suspended sediment. There are two typical reflectance peaks in the range of 500–600 nm and around 800 nm (Cai et al. 2015). The level of phytoplankton activity within a water body, which can be represented by phytoplankton biomass, is one of the most important considerations in water monitoring (Augusto-Silva et al. 2014). The changes of phytoplankton biomass are related to the concentration of phytoplankton's main and major pigment: Chla. The phytoplankton pigment exhibits a unique spectral absorption signature, with two marked peaks: one in the blue region of the spectrum (∼443 nm) and another one in the red region of the spectrum (∼686 nm) (Bernardo et al. 2016). Its main reflectance peak is in the green of the spectrum (∼550 nm), with a secondary peak at around 715 nm, which is associated with a high scattering of the phytoplankton cells (Gitelson 1992; Gitelson et al. 2007). According to the nature of the Crosta method and characteristics of TSM and Chla, the information on TSM and Chla contained in the PC can be selected appropriately.

This sequence continues for all p PCs. The first few PCs will tend to contain (or explain) a large percentage of the total variance and may be used to describe multi-variance patterns or variance in water quality. Often these patterns are related to specific sources of contamination (Olsen et al. 2012). For example, PCA is carried out in six bands of TM1, TM2, TM3, TM4, TM5 and TM7. The criteria for judging the PC of TSM is that the TM1, TM2, TM3, TM4, TM5 and TM7 coefficients of the PC should generally be positive, and the coefficients of TM1, TM2, TM3 and TM4 should be significantly larger than those of TM5 and TM7. The criterion for judging the PC of Chla is that the TM1 and TM3 coefficients of the PC should be significantly larger than those of the other bands, and the TM1 coefficient of the PC should generally be the opposite of that of the TM3 coefficient symbols.

An analysis of the correlation between the feature-orientated PC and ground truth TSM and Chla was conducted using Pearson's correlation coefficient, with a statistical significance of p < 0.0001. The significance level of 0 (less than 0.0001) indicated, in this study, that there are strong correlations between the feature-orientated PC and satellite-based TSM and Chla.

PCA was conducted using pre-processed Landsat TM images (excluding the TM6 band), with only water body pixels included. Feature-orientated PC selection involves examining the eigenvectors used in calculating the PCs and relating each PC image to the two or three original bands that contribute most data to it (Crosta & Moore 1989). Specific water components (e.g., TSM or Chla) can then be selected based on the major contributions (negative or positive eigenvectors) from the original bands that are likely to display the desired water quality parameters. The eigenvectors calculated for Landsat TM on October 17, 2003 are presented in Table 3. In Table 3, the eigenvectors are expressed as a percentage of the loading from the original bands. The transformation of the original eigenvector values to percentages makes the understanding of the original band contributions to each PC more straightforward. The knowledge on the spectral characteristics of TSM and Chla in the Landsat TM bands was used to define and select the PCs containing spectral information due to these water components.

As shown in Table 3, PC1 reveals a significant spectral feature, since it is composed of a positive mixture of TM1–TM4 bands, which, in this case, has a greater proportion of TM2 (29.91%) and TM3 (33.88%) showing TSM information. In the PRE, the reflectance curves of suspended sediment at different concentrations showed that with the increase of the suspended sediment concentration, an increase of the reflectance in the visible and NIR range can be clearly defined. At around 580 nm, the reflectance curves of the suspended sediment concentration can effectively differentiate different suspended sediment concentrations. Additionally, Xi and Zhang demonstrated that the red band is sensitive to changes of TSM and can be used to retrieve TSM concentrations successfully in the PRE (Xi & Zhang 2010). Thus, PC1 can be regarded as a measure of TSM information, according to the water volume. The distribution of the PC1 image is presented in Figure 2(a). PC2 is dominated by the contribution of TM4 (−25.70%), TM5 (−24.04%) and TM7 (−19.09%), showing a water body spectral response due to the high absorption in the NIR and short-wave infrared range. PC3 is dominated mainly by the contribution of TM1 (33.16%) and TM3 (−24.41%), showing Chla spectral information due to the high absorption in the blue and red range, since absorption features can be found in the blue and red ranges due to the high absorption in Chla (Bernardo et al. 2016). The distribution of the inverse PC3 image is presented in Figure 2(b). PC4 does not reveal any significant spectral feature, since it is composed of a positive and negative mixture of all bands, and several original bands contribute balanced proportions. PC5 is mainly formed by two negative loadings, from TM1 (−13.65%) and TM3 (−23.48%), and two positive loadings, from TM2 (24.76%) and TM4 (22.41%). As Chla has a high absorption in the blue and red range and a high reflectance in the green and NIR range (Bernardo et al. 2016), we are able to consider PC5, as a measure of Chla is left unexplained by PC3. PC6 is mainly formed by the contribution of TM5 (−43.13%) and TM7 (49.56%), which is considered to be spectral information on minerals, according to the water volume (Crosta & Moore 1989).

From the formula of the PC image pixel value, a substantial numeric value of an eigenvector (positive or negative) has a higher influence on the PC image. A low numeric value of the eigenvector means that the particular band has a relatively low significance for that PC. Here, we consider the effects of large (positive or negative) eigenvectors, as combined with a high or low reflectance value of the image (which, in turn, imply a strong reflectance or strong absorption, respectively). Low values of the eigenvectors are not considered, since these would have a minimal influence on the target PC image.

An image data set (Figure 3(a)) and PC1 subset (Figure 3(b)) of Landsat TM, taken on Oct. 17, 2003, are presented in Figure 3. Figure 3(c) shows the eigenvector of the PC1 image, and Figure 3(d) presents the spectral curve for the pixel under consideration. The spectral bands, TM1, TM2, TM3 and TM4, have very high positive eigenvectors, and all of the other bands have small positive eigenvectors (Figure 3(c)). The pixel in the crosshair, presented in Figure 3(a), is characterized by a very strong reflectance in TM2 and TM3 and a weaker reflectance in TM5 and TM7. Thus, in the resulting PC1 image, due to the combination of the large positive eigenvector in TM2 and TM3 and the very high reflectance value in the sub-scene of Figure 3(a), the pixel in the crosshair appears bright (Figure 3(b)). This means that the TSM concentration of the pixel is very high in the crosshair shown in Figure 3(a).

A second pixel, selected in the sub-scene (Figure 4(a)) and PC1 image (Figure 4(b)) of Landsat TM on Oct. 17, 2003, is presented in Figure 4. Figure 4(c) presents the spectral curve for the pixel under consideration. The eigenvectors of the PC1 image can be found in Figure 3(c). The eigenvectors of PC1 are the same as those of the pixel shown in Figure 4, while the reflectance values of all bands are smaller than the pixel shown in Figure 3. Thus, in the resulting PC1 image, the pixel in the crosshair appears dark (Figure 4(b)). In this case, the TSM concentration of the pixel is relatively lower than that of the pixel shown in Figure 3.

To evaluate the performance of the feature-orientated PC, 66 satellite-based samples from the MERIS FR level 2 TSM product were used to directly validate the results of the Crosta method. The distribution of the TSM and Chla extracted from the MERIS FR level 2 product are presented in Figure 5. The results show that the PC1 generated from the Landsat TM, taken on Oct. 17, 2003, and the MERIS TSM product are strongly correlated, with a highly significant linear relationship (R² = 0.83, p-value < 0.0001, RMSE = 4.77 g/m³) (Figure 6(a)). This indicates an acceptable level of consistency for the Landsat TM, taken on Oct. 17, 2003.

As for another water quality parameter, Chla, PC3, selected from the Landsat TM, is regarded as a measure of Chla, according to the water volume. To validate the performance of PC3, 66 satellite-based samples from the MERIS FR level 2 Chla product were used to validate the results of the Crosta method directly. The results show that there is a reliable correlation between PC3 and the MERIS Chla product (R² = 0.74, p-value <0.0001, RMSE = 2.88 mg/m³) (Figure 6(b)).

Considering the distribution of PC1 (Figure 2(a)), selected based on the feature-orientated PC from the Landsat TM, the spatial characteristics of the PC1 image are consistent with the MERIS TSM product (Figure 5(a)). From atmospherically corrected images, taken on Oct. 17, 2003 (Figure 3(a) and Figure 4(a)), the spatial distribution of suspended sediment can clearly be observed from visual inspection, which is also consistent with the spatial distribution characteristics of the PC1 image. The high values for the TSM concentration occurred on the west coast of the PRE, while the low values for the TSM concentration were present in the middle and on the east coast of the PRE. These spatial characteristics of TSM can be observed in both the PC1 image and the MERIS TSM product. The spatial distribution characteristics of the PC3 image are consistent with those of the MERIS Chla product. The validation results show a good correlation between the PC3 image and the MERIS Chla product. Overall, TSM and Chla in the PRE present similar spatial distribution characteristics, high values of TSM and Chla always occur on the west coast of the PRE, and low values are mainly distributed in the middle and on the east coast.

In order to test the stability of the Crosta method in identifying the water quality parameters (TSM and Chla) from the Landsat TM images, another seven Landsat TM images, which matched with the acquisition date of the MERIS FR level 2 product, were tested. The correlation coefficient (r), determination coefficient (R²) and RMSE of the feature-orientated PC and MERIS product (TSM and Chla) were calculated and are listed in Table 4. TSM, identified by the feature-orientated PC (TSM) using Landsat TM images, demonstrated that the Crosta method presented a good performance in the retrieval of the spatial distribution of TSM in the PRE. However, a poor performance (except for Mar. 4, 2008 and Oct. 17, 2009) was obtained for the feature-orientated PC (Chla) and MERIS-based Chla product. Overall, the feature-orientated PC (TSM) can effectively obtain the spatial distribution of TSM, while Chla is unstable using the Crosta method and Landsat TM images. It should be noted that the reflectance in the visible and NIR spectral range on surface water is very small in the remote sensing imagery. Additionally, the effect of the atmosphere makes the reflectivity of the water weaker. In the optically complex estuary and coastal waters, suspended sediment is the strongest constituent of scattering and reflection, which may characterize more information on remote sensing images. On the contrary, a high reflectance of suspended sediment masks the Chla information due to the strong absorption of Chla and pigments in the PRE.

In this study, the Crosta method, essential for PCA, was used to identify and select water quality parameters (TSM and Chla) from Landsat TM images and validate the results using MERIS FR level 2 TSM and Chla products. PCA is a dimension reduction technique, which can transform many original interrelated variables into fewer, uncorrelated variables, called PCs. The first few PCs provide meaningful information on some specific parameters. Each PC can represent the most meaningful parameter, according to the coefficients of the loadings on each original variable. Thus, each PC may contain the most meaningful parameter and other unremarkable parameter information due to the PC being a linear combination of original variables. In the visible and NIR spectral range, spectral absorption and a reflectance signature, with two or three marked peaks in TSM and Chla, exist. In addition, there is also an impact on unmarked peaks for TSM and Chla. Therefore, PC contained the most meaningful TSM and Chla information identified by the Crosta method and also included other information on water quality parameters, i.e., it did not only contain information on TSM. This information on TSM and Chla cannot be used for the spatial and temporal analysis of TSM and Chla in the PRE. Instead, it can only be used to analyze the spatial distribution of TSM and Chla qualitatively, when in-situ measurements are not available for deriving satellite-based water quality parameters. While previous studies demonstrated that the FLAASH atmospheric correction method was successfully applied in the retrieval of surface reflectance on Landsat TM/OLI in the estuary and coastal waters, there is also a need to develop an appropriate atmospheric correction algorithm in the specific estuary waters. An approach for atmospheric correction for the MERIS data collected over land was developed based on dense vegetation targets by Guanter et al. (2007) and can be used to retrieve water reflectance data from the MERIS data over optically complex inland waters (Guanter et al. 2010). In future studies, we could use this atmospheric correction method to retrieve water reflectance data from the Landsat TM/OLI data on the PRE.

In this study, the Crosta method was used to investigate the applicability of water quality parameters (TSM and Chla) identified by Landsat TM images. Feature-orientated PC TSM obtained a good performance, compared with the MERIS-based TSM product. However, the Chla information selected from the feature-orientated PC demonstrated that the Chla information derived from the Landsat TM using the Crosta method is unstable. A poor agreement can be observed between the feature-orientated PC Chla and the MERIS-based Chla product. Chla and pigments in the estuary and coastal waters are dominated by absorption, and the information obtained by remote sensing data is not rich. Additionally, this may be caused by the poor accuracy of the MERIS-based Chla product. Atmospheric correction may be another factor associated with this problem.

Overall, the Crosta method is an alternative method for the retrieval of spatial distribution information on TSM from Landsat TM images in the optically complex estuary and coastal waters, where in-situ measurements of TSM are not available in the study period. The findings presented in this study would be helpful in the retrieval of TSM information from the long-term historical Landsat image archive. In future studies, other atmospheric correction methods should be conducted and validated in relation to the Landsat TM/OLI images using in-situ measurements in the PRE.

This study was supported by the China NSFC (U1901215) and 863 Program (2006AA06A306), and the Scientific and Technological Innovation Programs of Higher Education Institutions in Shanxi (2019L0485). We would like to extend our gratitude to Jamie Robert Beagle for his help in editing the English language of this paper.