Study Area

Nansi Lake, a significant water body in Shandong Province, China, is made up of four interconnected lakes: Weishan Lake, Zhaoyang Lake, Dushan Lake, and Nanyang Lake, collectively referred to as Nansi Lake. This lake complex spans 1,266 Km², the area measures 126 km in the north-south direction and 5–25 km in the east-west direction. Its total storage capacity amounts to 4.7 billion cubic meters, making it the second-largest freshwater lake in the Huaihe River Basin of China. Nansi Lake plays a vital role as a water source and storage lake for the eastern route of the South-to-North Water Diversion Project of China, underscoring the importance of maintaining its water quality. This study focused on Dushan Lake, the uppermost lake in the Nansi Lake system, as the primary research area. In particular, we conducted hyperspectral remote sensing inversion to estimate the concentration of Chla of Nansi Lake. Fig. 1 displays the spatial distribution of the remote sensing images and in-suit points in the study area.

Satellite Images and Measured Data

ZY-1 02D Satellite Image

The ZY-1 02D satellite was launched by China’s Long March 4B carrier rocket at China Taiyuan Satellite Launch Center at 11:26 on September 12, 2019. The satellite was developed by the China Academy of Space Technology and has two payloads: a visible near-infrared and a hyperspectral camera. Table 1 displays the main parameters of the ZY-1 02D satellite’s hyperspectral camera payload. This satellite was located in the medium resolution, wide range observation, and quantitative remote sensing mission, which can provide rich spectral information on ground objects. Its hyperspectral camera subsystem was divided into 166 spectral channels in the spectral range of 400–2500 nm, with a visible near-infrared spectral resolution of 10 nm, a shortwave infrared spectral resolution of 20 nm, and the image covers an area of 60 km in width, with a spatial resolution of 30 m.

Fig. 1.

Study area and sampling points.

Considering the high spatial resolution, spectral resolution, and SNR, the comprehensive performance was at the leading level in the world. Its advantages could significantly improve the inversion effect and accuracy of water Chla. The image used in this study was the ZY-1 02D satellite image of the upper lake part of Nansi Lake, which was nearly synchronized with the field measurement time.

Technical Route

In this study, we employed a combination of satellite images and measured Chla concentration to estimate the Chla concentration in water. Firstly, we performed preprocessing on the acquired ZY-1 02D hyperspectral image, followed by spectral SG filtering to eliminate spectral noise. We compared the signal-to-noise ratio (SNR) of the image before and after noise reduction to assess the effectiveness of the denoising process. Based on the Chla measured concentration, the three-band method was applied to extract the characteristic bands of the original and denoised images. Subsequently, the XGBoost algorithm was utilized to model 70% of the dataset and then validate the entire dataset. A machine learning model was established to estimate Chla concentration. The accuracy of the image inversion before and after denoising was compared, and the technical route is shown in Fig. 2.

Savitzky-Golay

In the plane coordinate system, a curve is used to fit a set of data, assuming that this curve is

y = a₀ + a₁x + a₂x² + a₃x³ + a₄x⁴

here, the curve equation is used to substitute the abscissa of each point, and the objective is to minimize the sum of squared differences between the obtained value and the corresponding ordinate of the point. The curve has the highest fit, so that all coefficients can be determined a_p (p=0, 1, 2, 3, 4). The SG filtering method is a classic least-squares smoothing algorithm that uses a simplified least-squares fitting convolution to smooth the curve. Here, y_i represents the value of the i^th point before smoothing, ŷ_i represents the value after smoothing, the smoothing window size is 2m+1, and the fitting curve is a p-order polynomial:

Fig. 2.

Flowchart depicting the technical route.

The error of this polynomial fit curve is:

To minimize the error, the partial derivatives of ε are taken, and the partial derivatives of ε are zero, and the following linear equations are obtained:

By solving the linear equation to obtain a₀, a₁, ..., a_p, a fitted polynomial curve can be obtained.

Three-Band Combined Parameter Construction of Chla Concentration

For the inversion of Chla concentration in water, three main methods are used at home and abroad: empirical, bio-optical, and semi-empirical/semi-analytical methods. Each of the three methods has advantages and limitations. Among the three methods, the three-band semi-analysis (TBS) method, based on the semi-analytical method proposed by Dall’Olmo [32] is the most widely used and the most common method for obtaining the Chla concentration [33, 34].

Gitelson et al. [35] suggested a method that utilizes a three-band sensitive spectral combination extraction technique to estimate the chlorophyll content in terrestrial plants, the form of which is:

and ≥730 nm as the value ranges of λ₁, λ₂, and λ₃, respectively, and the parameter model based on a three-band method was developed for the inversion of Chla concentration in water. The model’s mathematical expression is as follows:

where C_chl-a is the concentration of Chla; R^-1(λ₁) and R^-1(λ₂) are the reciprocals of the reflectance at wavelengths λ₁ and λ₂, respectively; R (λ₃) is the reflectance at wavelength λ₃; A and B are indeterminate coefficients.

Based on the three-band method, this study used the three-dimensional spectral index construction method to identify the position of the feature bands. This index construction method selects characteristic bands by collectively examining three bands through traversal. It traverses the combined correlations between all bands through iterative operations, which fully considers the synergy between bands and minimizes the influence of irrelevant wavelengths [37]. The correlation between the image spectrum and measured Chla concentration was determined using the three-dimensional spectral index construction method, which resulted in different characteristic band positions compared to those selected by the traditional semi-empirical/semi-analytical three-band method.

XGBoost Model

According to the data characteristics of this study, we selected the XGBoost algorithm to perform multivariate regression and fit the dependent variables. XGBoost was designed by Chen Tianqi. When using CART as the base classifier, a regularization parameter is employed to regulate the model’s complexity [38]. By effectively reducing overfitting of the data, the XGBoost model improves the generalization ability of the data. The model’s mathematical expression is as follows:

where ŷ_i is the predicted value of the XGBoost model of Chla concentration, the parameter value x_i represents the i^th three-band combination, K is the number of trees, f_k is the k^th tree model.

Accuracy Verification

In this study, the Pearson correlation coefficient was used as the basis for selecting the correlation between the characteristic band parameters and measured Chla concentration. The Pearson correlation coefficient, also known as the Pearson product-moment correlation coefficient, was employed to assess the linear correlation between X and Y in the two datasets. Its value ranges between −1 and 1. The calculation formula is shown in (7).

where E is the mathematical expectation; cov(X,Y) represents the covariance of X and Y; σ_X and σ_Y are the standard deviations of X and Y, respectively. The closer the value of ρ_X,Y is to one, the higher the correlation between the two datasets.

In this study, the model fitting accuracy was evaluated using the coefficient of determination (R²), which compares the predicted values with the actual measured values. The R² is obtained by dividing the explained variation by the total variation, which represents the quality of the regression line fit. Its value ranges between 0 and 1. The degree of accuracy of the regression equation is higher when R² is closer to 1 and worse when it approaches 0. The accuracy of the model was evaluated in this study using the root-mean-square error (RMSE), which quantifies the extent of the difference between the predicted and actual values. A lower value of RMSE implies a higher level of accuracy, as there is less discrepancy between the predicted and measured values. Formula (8) was used to calculate R², while (9) was used to calculate RMSE.

Fig. 3.

SNR and RMSE. (a) Curves of the calculation results for different image SNR. (b) RMSE of the denoised image relative to the original image.

where n represent the number of samples, y_i and y_j are measured and predicted values of the Chla concentration, respectively; y¯i and y¯j are the mean values of the measured and predicted Chla concentrations, respectively.

ZY-1 02D Image Denoising

The denoising of the preprocessed ZY-1 02D image using the SG filtering method is discussed. Filter window coefficients m and n can be selected over a relatively wide range. Generally, multiple iterations are required for different data to find the optimal kernel parameter range. The spectral curve used in this experiment was relatively complex. To preserve the original spectral curve characteristics to the greatest extent and achieve the best image denoising effect and inversion result, extensive iterative tests were conducted. Based on the experimental results, the SG filter window coefficients m = 7, n = 5 (7–5D), and m = 9, n = 5 (9–5D) were retained as the image results after noise reduction. Subsequently, they were used to compare the inversion accuracy of Chla concentration with the data obtained from the original image (OD).

To verify the denoising effect of the test image, the SNR was objectively measured. The higher the SNR, the better the denoising effect. The three image datasets were imported into the ENVI5.3 IDL console, and the SNR was calculated using the local standard deviation method based on edge elimination through code. The SNR results are shown in Fig. 3(a). The SNR of 7–5D and 9–5D were slightly improved compared with OD, but the denoising effect varied in different bands.

To further examine the similarity between the spectrum of the denoised image and the original image, their RMSEs in each band of the sampling point spectrum were also calculated. The smaller the RMSE, the closer the denoised image spectrum is to the original image spectrum. The RMSE of the denoised image relative to that of the original image is shown in Fig. 3(b). The RMSE between 7–5D and OD was smaller compared with 9–5D, indicating that 7–5D had less image information loss. The image quality has a certain influence on the inversion of the Chla concentration in water. Then, the characteristic bands of the OD, 7–5D, and 9–5D data were extracted to build the optimal concentration estimation model.

Fig. 4 shows the reflectance curves of sampling points OD, 7–5D, and 9–5D. It can be observed that the spectral curve of the denoised image had a significant improvement in the fine sawtooth compared with the original image.

Characteristic Band Extraction

Specifically, the extraction of feature bands is divided into four steps: (1) Extract spectral data from the original image and the denoised image (7–5D and 9–5D) according to the actual measurement point positions. (2) Calculate the correlation coefficient between a single band and the measured Chla concentration. (3) Follow the method in Technical Route Section, iteratively cycle through all possible band combinations, and calculate the correlation coefficient between all combinations and the measured Chla concentration. (4) Sort all correlation coefficients and select the five largest characteristic bands as input to the model. The correlation coefficient statistics of the first five groups of band combinations between the three groups of spectral data and the Chla concentration are shown in Table 1. Simultaneously, the SNR of the sampling points at the corresponding characteristic wavelengths in OD, 7–5D, and 9–5D were extracted. It can be seen from Table 1 that the SNR of the characteristic bands extracted at 9-Dd and 7–5D was slightly improved.

Fig. 4.

The reflectivity curve of a sampling point.

Hyperspectral Image Retrieval of Chla Concentration in Nansi Lake

The method used in this study was the modified normalized difference water index (MNDWI) proposed by Xu [39]. Water bodies were extracted from OD, 7–5D, and 9–5D images, respectively. Based on the combination of the five bands with the best correlation, the three-band formula was used to calculate the inversion parameters of the three sets of images. Using the XGBoost model, 70% of the sampling points were selected for modeling, and the inversion accuracy of the Chla concentration was verified at all sampling points. Fig. 5 shows the image inversion results. The high and low values of Chla concentration in the original image were relatively concentrated, and the high-value areas accounted for a large proportion. With the adjustment of the image noise reduction parameters, the distribution of high and low values of Chla became more dispersed. Fig. 6 shows a fitting diagram of the image inversion accuracy. The inversion accuracy of the Chla concentration in the image after noise reduction was significantly improved compared to the original image. Among them, using 7–5D image data to invert Chla concentration had the highest accuracy, with an R² of 0.8737 and RMSE of 7.45 μg·L^-1. In addition, the Chla concentration in Dushan Lake is generally high. During the field investigation, it was found that the lake had undergone multiple large-scale drainage just before the satellite transit and actual measurement. The direction was from south to north, causing the water flow to concentrate towards Dushan Lake, causing pollution. In addition, there are many rice and terrace fields around Dushan Lake, and the soaking and rotting of straw after returning to the fields, pesticide, and fertilizer residues, etc., will also affect the water quality of the lake, so the concentration of Chla in Dushan Lake will be higher.

Fig. 5.

Chla concentration image inversion results.

Selection of Feature Bands in the Process of Chlorophyll-a Inversion

Regarding the selection of characteristic bands, in the study of Liu et al [40], the authors refer to the previous research experience and, through the traditional band positioning method, determined the central wavelengths of λ₁, λ₂, and λ₃ from the range of band selection of the three-band semi-empirical model to be 671 nm, 705 nm, and 731 nm respectively. This is different from the three-band parameter model constructed by the three-dimensional spectral index method based on the three-band method model formula in this study to invert the concentration of Chla. The selected three-band positions were quite different because the three-dimensional spectral index method adopted in this study did not consider the selection range of each band in the three-band semi-empirical model, but the overall loop iterative traversal of all the bands to achieve the positioning of the band positions.

Necessity of Hyperspectral Denoising and Selection of Inversion Model

The inversion accuracy of the original image was low, which may be caused by the deviation between the image acquisition time and sampling time. In addition, the direct error between the measured point and the pixel position of the extracted image spectrum, radiometric calibration of the image, atmospheric correction, and other errors affect the inversion results of the image. Therefore, in the study of water Chla inversion using hyperspectral images, it is necessary to optimize further the image preprocessing method to reduce the influence of various errors on inversion accuracy. In terms of the inversion model, XGBoost model was selected to invert the Chla concentration of images to alleviate the overfitting phenomenon caused by the machine learning algorithm due to the small number of sampling points used in this study. This is because previous studies have found that the XGBoost model can be trained on sparse data using the sparse-aware split-finding method. The XGBoost model will always have better results in water quality inversion, which is similar to the results of previous studies, the application of this method holds significant importance in remote sensing image processing, particularly in addressing the challenge of low sampling rates in field data.