Quantitative estimation of wastewater quality parameters by hyperspectral band screening using GC, VIP and SPA

Sample preparation and chemical analysis

The water samples were taken from different spots (water inlet, anoxic tank, aerobic tank, sedimentation tank and water outlet) under different treatment methods at a domestic sewage treatment plant in China. The domestic sewage treatment plant gave field permit approval to us (NO. 51409221, 51979234). According to our pre-analysis (Table 1) (CESP, 2002; Kotti et al., 2018), the water samples collected at the inlet were abundantly rich in COD, BOD and NH3-N (COD: 425 mg/L, BOD: 101.15 mg/L, NH3-N: 34.853 mg/L), whilst the water samples collected in other treatment ponds were relatively low in COD of 20–140 mg/L, BOD of 3.1–13 mg/L and NH3-N of 0–1.723 mg/L. Such large variations suited well to represent the range of water quality collected from different locations investigated in this study, but also leaving non-negligible deficiency in representing water quality of moderate concentrations. Therefore, to make the hyperspectral model more robust and responsive to minor changes of water quality at other potential locations, we further diluted the water samples into subsamples to form a full range of distribution within the minimum and maximum values detected in this study. To be specific, following the minimum and maximum COD measured from the inlet and other treatment ponds, the water samples rich in COD were diluted by adding different amount of pure water, forming different COD concentrations varying at unit of five mg/L. In total, 87 water samples were obtained to establish the hyperspectral model and strengthen its applicability. In a real application of hyperspectral model, the water quality of targeted samples often covers a variety of concentrations, and thus can be measured and applied directly without any dilution. The laboratory was kept in a constant temperature during the experiments. The main water quality parameters included NH₃-N, TA, TH, TDS, COD and BOD.

Acquisition and pretreatment of spectral data

The wastewater samples were put into black cylindrical cups with a depth of five cm and a diameter of 10 cm for spectral data collection in the laboratory. The hyperspectral data for wastewater samples were measured using an ASD (Analytical Spectral Devices, Inc., Boulder, CO, USA) FieldSpec^® 3 spectrometer with spectral range from 350 to 2,500 nm. The instrument was equipped with one sensor with 1.4 nm spectral resolution for the range of 350–1,000 nm and the other sensor two nm, 1,000–2,500 nm. The spectral data were collected in a dark room by exposing the water sample surfaces to a halogen lamp of 50 W above with a 30° of incident angle and 50 cm in distance and a fiber-optics probe 5°, 15 cm, so as to minimize the effects of external factors (Wang & Zhao, 2000). Before each measurement, we had fully preheated the spectrometer and light source and checked the spectrometer with a standardized white panel of 99% reflectance in order to reduce measurement errors. Each water sample was measured twice vertically, and at each of which the spectral data were gathered 10 times. Altogether there were 20 spectrum curves for each sample (Hong et al., 2018). From the 20 curves the raw spectral reflectance (R_raw), namely, arithmetic mean, were obtained by using ViewSpecPro software (6.0 version). The fluctuation of R_raw would influence the accuracy of subsequent modeling due to the disturbance of the random error, instrument noise and external environment in spectral data collection. In general, the external noise can be eliminated to some degree with such effective pretreatments as resampling, smoothing and transformation, which can improve the spectral characteristics (Ding et al., 2018). Therefore, two steps were adopted to pretreat the R_raw: (1) removing the marginal wavelengths (2,401–2,500 nm and 350–399 nm) of higher noise in each water sample, then smoothing the remaining spectrum data through filter method (polynomial order = 2; window size = 5) of Savitzky-Golay (Savitzky & Golay, 1964) by Origin Pro software (2017SR2 version); and (2) obtaining the precise R_raw−SNV using the standard normal variable (SNV) to remove the effects of baseline shift and surface scattering on the spectral data (Xiao, Yi & Hao, 2016). The spectral curves of R_raw and R_raw−SNV are shown in Figs. 1A and 1B. As can be seen, the spectral curve in Fig. 1B is much smoother than that in Fig. 1A.

Wavelength selection methods

Variable importance projection

The VIP is a variable screening analysis method based on PLSR model. The VIP value of a spectral variable reflects the importance of this variable in the prediction of the substance to be measured (Oussama et al., 2012; De Almeida et al., 2013). The VIP value reflects the explanatory power of independent variable over dependent variable, and represents the importance of independent variable to model fitting (Chavana-Bryant et al., 2019). If the explanatory ability of each variable to Y is the same, then the VIP value of all independent variables is 1. If the VIP value of an independent variable is less than 1, it means that the variable makes little contribution to the model and has a very low ability to explain the dependent variable, so it can be eliminated. In the PLSR model, VIP_i to an independent variable X_i is defined as: (4) $$\mathrm{V}\mathrm{I}{\mathrm{P}}_i=\sqrt{\frac{p\ast {\sum }_{h=1}^{m}Rd(Y;th)W_{hj}^2}{Rd(Y;t_1,t_2,\cdot \cdot \cdot ,tm)}}$$ where p is the dimension of independent variable; t_h is the component h; m is the number of components selected in the model; W_hj is the component of X_j corresponding to axis W_h; Rd (Y; t_h) represents the variation of Y explained by the component t_h; and Rd (Y; t₁, t₂, ···t_m) represents the cumulative variation accuracy of Y explained by t₁–t_m. In PLSR model, the explanatory power of x_j over Y is passed through t_h, so if t_h bears a strong explanatory power over Y, and meanwhile x_j plays a very important role in the process of t_h construction, t_h can be reasonably regarded as very important variable in explaining Y. In this way of wavelengths selection, the wavelengths with strong explanatory power are reserved while those with weak explanatory power are eliminated (Chemura, Mutanga & Dube, 2017).

Set pair analysis

Set pair analysis integrates uncertainty analysis and determinacy analysis (Zhao & Xuan, 1996). The basis of SPA is set pair whose key is correlation degree, which has been applied as a criterion for the analysis of certainty and uncertainty (Zhang et al., 2019). During the hyperspectral band selection, the GC value Y_i = (y₁, y₂,···, y_n) can be defined as set A, and VIP value Z_i = (z₁, z₂, ···, z_n) as set B, and then the two sets compose a set pair H = (A, B). To study the correlation degree of the set pair, the formulas are as follows: (5) $${\mu }_i=S_i+F_i\ast I_i+P_i\ast J={\displaystyle \frac{s_i}{\mathrm{n}}+{\displaystyle \frac{f_i}{\mathrm{n}}{I}_i+{\displaystyle \frac{p_i}{\mathrm{n}}J}}}$$ (6) $$I_i={\displaystyle \frac{S_i-P_i}{S_i+P_i}}$$ (7) $$v_i=1/n+1/n\ast u_i$$ (8) $$W_i={\displaystyle \frac{v_i}{{\sum }_i^N{v}_i}(i=1,2,\dots ,N)}$$ where S_i, F_i, P_i are the identity degree, discrepancy degree and opposite degree of the two sets in the context of the same problem, respectively. They describe the association of the two sets from different aspects. S_i, F_i, P_i meet the relationship: S_i + F_i + P_i = 1, n denotes the total number of characteristics of the set pair; s_i is the number of the common characteristics of the two sets; f_i is the number of the different characteristics of the two sets (Li et al., 2016); p_i is the number of the opposite characteristics of the two sets. When J = −1, the contact number u_i, the relative membership degree v_i, and the weighted value, W_i, of the VIP value and GC value, were calculated.

Model construction and validation

Correlation between water quality parameters content and spectral reflectance

The Pearson correlation coefficients between each water quality parameter and the R_raw−SNV (400–2,400 nm) were tested with the significance level (P < 0.001, |r| = 0.324 or above). The curves of correlation coefficients of water quality parameters were plotted in Fig. 2.

The change of the correlation coefficients of water quality parameters was complex. There were significant differences in the correlation coefficients between different water quality parameters and the wavelengths but their curve patterns were similar (Peng et al., 2018). There were many concentrated similar sensitive ranges but in general, the curves took “sharp fall, fluctuation, and oscillation” (Fig. 2). The COD, BOD, NH₃-N and TA were all highly correlated and represented a positive correlation (the correlation coefficients > 0.6) at the wavelengths of 900–1,300 nm and 1,500–1,800 nm. In contrast, the correlation between the TDS and spectral data and that between the TH and spectral data were not so satisfied (the correlation coefficients <0.6). The six parameters could pass the test of significance (P = 0.001) at the wavelengths of 400–450 nm and 850–1,800 nm. The band range was wide and the peak value of the correlation coefficients was relatively flat, so it was difficult to screen the characteristic bands based on the correlation coefficients.

Selection of characteristic wavelength

GC-based selection of characteristic wavelength

The curves of GCD for water quality parameters and R_raw−SNV are shown in Fig. 3. An obvious peak wavelength appeared in the curves with the similar patterns of “sharp rise, fall and stabilization” on the whole (Fig. 3). The gray value of each water quality parameter gradually went up from 400 nm to a peak at about 820 nm, and then gradually down to 1,000 nm where it remained flat. The overall gray values of TDS and TH were relatively high, generally above 0.5, while those of COD, BOD, NH₃-N and TA were generally between 0.3 and 0.4 with only a few peak bands greater than 0.4 and above 0.8 at the most. It reflected that the GC method excels in eliminating a great amount of redundant information and selecting a few most sensitive bands.

The sensitive bands of the six parameters were counted based on their GCD analysis (Table 2). Of all the GCD values, the maximum ones were concentrated within the spectral wavelengths of 815–830 nm. The numbers of the sensitive bands were sequenced from large to small as: TH (601) > TDS (381) > TA (93) > BOD (50) > COD (49) > NH₃-N (46), and those of the maximum GCD values as: TH (0.8504) > TDS (0.802) > BOD (0.7974) > NH₃-N (0.7973) > COD (0.7949) > TA (0.7878). The comparison between the correlation coefficients and the GCD indicated a great difference: the higher correlation coefficient the parameters had, the lower the GCD values was, and the smaller number of sensitive bands were selected. This result revealed that for the parameters with strong spectral response the GCD method had better band screening ability.

Table 2:

Maximum gray correlation degree and band intervals of water quality parameters content with standard normal variable reflectance.

Water quality parameters	Sensitive band numbers	Maximum GCD	Maximum GCD band intervals (nm)
COD	49	0.7949	820–830
BOD	50	0.7974	820–830
NH3-N	46	0.7973	820–830
TA	93	0.7878	820–830
TDS	381	0.802	815–825/766
TH	601	0.8504	815–825

DOI: 10.7717/peerj.8255/table-2

VIP-based selection of characteristic wavelength

The curves of VIP scores for water quality parameters and R_raw−SNV are shown in Fig. 4. The curves patterns of the six parameters were similar and the overall scores of the VIP were relatively high (Fig. 4). These curves exhibited a sharp drop in the wavelength intervals of 400–450 nm, a sharp rise in 460–500 nm, a violent fluctuation in 700–1,000 nm, but in 1,100–2,000 nm TH and TDS displayed a less violent fluctuation and COD, BOD, NH₃-N and TA showed flat fluctuation. Many peak points existed in the VIP curves, and their peak intervals were relatively scattered, completely different from the single peak value in the GC curves. The principle of VIP > 1 was used to select and count the spectral sensitive bands of the six parameters (Table 3). The comparison of the maximum spectral response bands showed five parameters were found in the wavelength intervals of 460–480 nm except for NH₃-N in 990–999 nm.

Table 3:

The maximum VIP scores of water quality parameters with standard normal variable reflectance.

Water quality parameters	Sensitive band numbers	Maximum VIP scores	Maximum VIP band interval (nm)	Maximum VIP band (nm)
COD	753	1.634	465–475	466
BOD	770	1.439	464–474	762
NH3-N	768	1.397	990–999	999
TA	709	2.466	460–470	464
TDS	497	4.275	460–470	463
TH	543	4.893	460–470	463

DOI: 10.7717/peerj.8255/table-3

The numbers of sensitive bands selected through VIP method were sequenced from large to small as BOD (770) > NH₃-N (768) > COD (753) > TA (709) > TH (543) > TDS (497), and those of the maximum VIP scores as: TH (4.893) > TDS (4.275) > TA (2.466) > COD (1.634) > BOD (1.439) > NH₃-N (1.397). The comparison between the two sequences indicated a great difference: the larger the VIP score was, the smaller number of selected sensitive bands was. On the whole, the effective bands could be retained as many as possible with the VIP method during bands screening.

SPA-based selection of characteristic wavelength

The curves of SPA scores for water quality parameters and R_raw−SNV are shown in Fig. 5. The curves patterns of TDS and TH were similar and exhibited a gentle flatness except for obvious peaks at the wavelengths of about 460 and 990 nm. Those of BOD, NH₃-N and TA were similar and showed notable valleys at the wavelengths of about 500 and 700 nm and notable peaks at the wavelengths of about 760 and 1,000 nm. Those of COD displayed notable peaks at the wavelengths of about 460 and 1,000 nm and fluctuation at other wavelengths. On the whole, great difference existed among the six curves, and the sensitive interval of each parameter was quite notable and different from each other.

The principles of GCD > 0.5 and VIP > 1 were used to select and count the spectral sensitive bands of the six parameters (Table 4). The comparison of the maximum spectral response bands revealed NH₃-N was in the wavelength intervals of 990–999 nm, BOD in 760–770 nm and the other four in 460–475nm. The numbers of sensitive bands selected through SPA method were arranged in a descending order as: BOD (767) > NH₃-N (765) > COD (750) > TA (696) > TDS (280) > TH (223), the scores of the maximum SPA as: TA (2.409) >TH (1.944) > TDS (1.928) > COD (1.598) > BOD (1.409) > NH₃-N (1.364). The comparison between the two orders indicated a great difference: the larger the SPA score was, the smaller number of the selected sensitive bands had.

Table 4:

The Set pair analysis (SPA) scores of water quality parameter with standard normal variable reflectance.

Water quality parameters	Sensitive band numbers	Maximum SPA	Maximum SPA band interval (nm)	Maximum SPA band (nm)
COD	750	1.598	465–475	466
BOD	767	1.409	762–763	762
NH3-N	765	1.364	990–999	999
TA	696	2.409	460–470	464
TDS	280	1.928	460–470	463
TH	223	1.944	460–470	463

DOI: 10.7717/peerj.8255/table-4

Construction and analysis of PLSR model

The sensitive bands selected by different methods as GC, VIP and SPA were applied to PLSR modeling. The results of PLSR model are shown in Table 5.

In general, the PLSR model had a good spectral prediction of water quality parameters (Wang et al., 2019). The GC-PLSR model was proved to be optimal for COD, BOD, NH3-N and TH with the best prediction. SPA-PLSR model was proved to be optimal for TDS and TA with the best prediction (Table 5). The VIP-PLSR model of BOD had the best modeling effect (R²_c = 0.99). The GC-PLSR model of NH₃-N had the best prediction effect (R²_p = 0.962, RPD = 5.894). The inversion effects of TDS and TH were relatively less satisfying, with the R²_c and R²_p generally around 0.8. Overall, the PLSR models based on the characteristic bands of the six parameters exhibited good modeling and prediction effect as well as a greater improvement than those full band-based models.

Construction and analysis of ELM model

The sensitive bands selected using GC, VIP and SPA methods were applied to build ELM model. The results of ELM model are shown in Table 6.

The ELM models based on the characteristic bands of the six parameters exhibited good modeling and prediction effect as well as a greater improvement than those based on the full band. The GC-ELM model was proved optimal for COD, BOD and NH₃-N with the best prediction. SPA-ELM was proved optimal for TDS, TA and TH with the best prediction. The VIP-ELM model of BOD had the best modeling effect (R²_c = 0.986). The GC-ELM model of NH₃-N exhibited the best prediction effect (R²_p= 0.976, RPD = 6.596). The ELM models of COD, TH and TA had the satisfying modeling and validation effect. The validation effect of TDS was relatively less satisfying (R²_c = 0.820, R²_p = 0.790, RPD = 2.198).

Comparison among the estimating results of different water quality parameters

The optimal wavelength selection methods varied when the optimal modeling methods were different (Tables 5 and 6). The hyper-spectrally estimated value and the chemically measured value of the six water quality parameters were compared under the optimal model (Fig. 6).

Different water components can change the spectral radiation by their own absorption and scattering characteristics of light. The water components (such as SPM, CDOM, phytoplankton, blue-green algae) can be obtained by separating the factors of various parameters from the total water radiation (Dörnhöfer & Oppelt, 2016). However, the problem of “lower ion content and weaker spectral response” often exists in the spectral analysis of water components, and was also reflected in the verification results of this study. All the six parameters produced satisfying inversion results, which indicated the feasibility of quantitative inversion of sewage water quality parameters by hyperspectrum (Pu et al., 2017). The sequence of the predicting power of the water quality parameters was NH₃-N > BOD > COD > TA > TH > TDS. The validation results showed that most data points of the five water parameters, NH₃-N, BOD, COD, TA and TH, were concentrated near line 1:1. This indicated that the hyperspectral analysis values of the five parameters were very close to the chemically estimated values and the optimal models of these five parameters had strong prediction ability (RPD > 3.0) (Tables 5 and 6). Comparatively, due to its lower content, the data points of TDS were relatively discrete, indicating that the model prediction ability was average (RPD = 2.198).

Correlation analysis and inversion performance

The raw spectral reflectance curve of each water sample exhibited different shapes (Fig. 2A). This difference results from the inconsistency between the contents and types of the water quality parameters, and the different characteristics of spectrum obtained through absorption and scattering of light. The results accord with those in previous studies (Abdelmalik, 2018; Abd-Elrahman et al., 2011; Rostom et al., 2017; Wang, Pu & Sun, 2016), which demonstrates that the VIS-NIR spectrum can be used to quantitatively determine water quality parameters.

Traditionally, correlation analysis is helpful to reveal the relationship between water quality parameters and the VIS-NIR spectrum (Peng et al., 2018). In this paper, the numbers of the significant wavebands of the six water quality parameters can be arranged in an order, in which the numbers of NH₃-N, BOD and COD were close to one another but larger than those of the other three parameters. The number of TA followed, and those of TH and TDS came last. This order agreed with that of their correlation coefficient ranges. Therefore, significant correlations exist between NH₃-N, BOD, COD and their reflectance spectrum, which reflects the optimal models of the three parameters have a very good prediction performance (Fig. 6). By contrast, the correlations between the other three water quality parameters and their reflectance spectrum were relatively low, which led to poor model prediction performance, especially for the TDS model. However, by comparing the three screening methods, the numbers of the sensitive wavebands of TA, TH and TDS were not as small as had been expected (Tables 3, 4 and 5). This is because different waveband screening methods use different calculation mechanisms, which leads to different methods to screen optimal response spectrum for different parameters (Zhang, Maoguo & Yongqiang, 2018). From the perspective of modeling methods, ELM model performed better than PLSR model, which also indicated the superiority and strong learning ability of machine learning algorithm (Keller et al., 2018). From the perspective of band selection methods, GC and SPA are superior to VIP. Therefore, the selection of band screening methods will be of practical significance for the quantitative inversion of water quality parameters in the future research.

Effects of wavelength selection on hyperspectral estimation models

In the process of hyperspectral analysis, due to the mass of hyperspectral information, there exists a large amount of redundant information irrelevant to parameters content. This makes it difficult to mine and extract effective information. This problem complicates the hyperspectral analysis and decrease the analysis accuracy, which hinders the development and utilization of related testing equipment. Therefore, characteristic selection of hyperspectral data information is required. The characteristic selection methods are devised to find the smallest subset from the original characteristics without reducing the accuracy (Zhang, Maoguo & Yongqiang, 2018). For the hyperspectral band selection, there are two key aspects, namely, effective information preservation and redundancy elimination, which has a great impact on subsequent modeling and validation.

In this study, GC and VIP methods were used to explore their spectral characteristic response intervals for six main water quality parameters in wastewater. These two methods had their fixed preferences in terms of information preservation and redundancy elimination. Therefore, the two methods had quite different performance in the band selection for the parameters and modeling methods, and their performance was very unstable (Wang et al., 2019). Of the two methods, the GC method had better performance in redundancy elimination. Better spectral response of the water quality parameters means smaller number of bands filtered by the GC method, so the GC method is the best screening method for COD, BOD and NH₃-N with a satisfying screening performance in this study.

However, the VIP method performed better in information preservation, and it could retain as many spectral band intervals with good response to the inversion parameters as possible (Oussama et al., 2012; De Almeida et al., 2013). So, the VIP method retained the largest number of bands during band selection. Meanwhile, compared with the GC method, VIP was more applicable to characteristic band selection of TA and TDS.

As to this condition, some researchers considered it necessary to make the optimal trade-offs between information preservation and redundancy elimination according to the different characteristics of the datasets, and have conducted some explorations (Zhang, Maoguo & Yongqiang, 2018). This study put forward a band screening method based on SPA to optimize the trade-offs between GC and VIP. From its application to the datasets in this study, the SPA method could well balance effective information retaining and redundant information eliminating. Compared with GC and VIP, the SPA method had a more stable ability of band screening. It had a strong applicability to predicting the six water quality parameters as well as building the linear modeling method PLSR and the nonlinear modeling method ELM. Therefore, the SPA method had the superiority in certainty and uncertainty analysis of characteristic bands, and it is applicable to the integrative optimization of GC and VIP.

Research limitations and future research

This study was only the statistical application of characteristic band screening methods to the inversion of different water quality parameters and failed to explore the spectral response mechanism of the various water quality parameters (Hu & Wang, 2017), which would lead to some limitations in the applicability of the model. Meanwhile, the optimized band screening method based on the SPA needs further study, which can focus on further:

1. Exploring the different spectral sensitive intervals and analyzing the spectral response mechanism of each water quality parameter.

2. Probing into the band screening mechanism and optimizing various screening methods so as to propose a band screening method suitable for various models; and applying various screening methods at the same time in the longitudinal study to achieve multiple screening.