A Deep Multitask Semisupervised Learning Approach for Chlorophyll-a Retrieval from Remote Sensing Images

Abstract

This article addresses the scarcity of labeled data in multitemporal remote sensing image analysis, and especially in the context of Chlorophyll-a (Chl-a) estimation for inland water quality assessment. We propose a multitask CNN architecture that can exploit unlabeled satellite imagery and that can be generalized to other multitemporal remote sensing image analysis contexts where the target parameter exhibits seasonal fluctuations. Specifically, Chl-a estimation is set as the main task, and an unlabeled sample’s month classification is set as an auxiliary network task. The proposed approach is validated with multitemporal/spectral Sentinel-2 images of Lake Balik in Turkey using in situ measurements acquired during 2017–2019. We show that harnessing unlabeled data through multitask learning improves water quality estimation performance.

1. Introduction

Besides containing invaluable aquatic habitats and rich biodiversity, lake water also constitutes a significant resource with a wide range of civil and industrial purposes, ranging from urban water supply and agricultural irrigation to fishery and recreation. Consequently, monitoring its quality is of paramount importance, especially given the growing threat from uncontrolled farming practices (e.g., excessive fertilizer and pesticide use) and industrial pollution [1].

Chlorophyll-a (Chl-a) concentration is one of the most commonly used water quality parameters. Its high levels indicate a state of eutrophication, most often due to an abundance of nutrients [2]. It is known to not only cause the death of aquatic life but also constitutes a threat to public health [3]. The traditional measurement of Chl-a concentration is a prolonged and labor-intensive process. It involves arduous and repeated sample acquisitions from the field and their subsequent laboratory analysis. The end results, in addition, only reflect the water quality of a handful of measurement sites at most, at often large temporal intervals, or even fewer, in cases where the lake is geographically inaccessible [4]. Therefore, developing a computer model that employs as input remote sensing images acquired over the area of interest, and that outputs pixel-level Chl-a estimations, constitutes a highly attractive and efficient alternative solution, capable of the systematic monitoring of an entire water body, rendering both field visits and laboratory analysis redundant.

The research question underlying this study is how to develop such a computer model, capable of exploiting the scarce labels as well as abundant unlabeled data.

1.1. Related Work

The assessment of water quality from remote sensing images in terms of Chl-a concentration constitutes a long-standing challenge, with an abundance of published approaches spanning open oceans [5], and coastal [6] and inland waters [7]. Chl-a retrieval algorithms can be aggregated into two broad categories: semianalytical and semiempirical [8]; for a comprehensive survey, the reader is referred to [9].

Semianalytical approaches rely on radiative transfer theory and modeling the propagation of light in water, through the radiative transfer equation that connects inherent optical properties of lake water with its radiance levels [10]. Solving this equation via either look-up tables or inversion methods can lead to the estimation of Chl-a concentration. Even though they can be developed with no available field samples, and are broadly applicable thanks to their physical foundation [11], they still require ample prior information about the lake’s optically active water constituents. However, their main disadvantage is their high sensitivity to atmospheric effects, as the models rely heavily on precise radiance readings, thus limiting their operational capacity [12].

Semi-empirical approaches, on the other hand, are more robust against atmospheric effects and do not require prior information about a lake’s physical characteristics [9]. They rely on feature engineering; in other words, each pixel of an often-optical remote sensing image is described through a numerical feature vector that is subsequently provided to a statistical/machine learning algorithm with the end goal of developing a regression model for Chl-a estimation. Notable examples include Gaussian processes [13], multilayer perceptron [14], and support vector regression [15]. Comprehensive comparative studies of such features have been provided in [16]. As features, various linear and nonlinear band combinations can be encountered in the state-of-the-art [17], involving mostly blue–green and red–near-infrared bands [18], as well as various expertly derived band ratios and indices [19]. Naturally, they require extensive in situ measurements of Chl-a in order to produce models with high estimation accuracy.

The advent of satellites such as Sentinel-2, with shorter revisit times and higher spatial and radiometric resolutions, has not only paved the way for new advances, but has also rendered necessary the use of more elaborate methods, to deal with the greater level of image detail [20]. More precisely, the results of an extensive comparative study across various band ratios, using Sentinel-2 multispectral images with multiple atmospheric correction processors, has been presented in [16]. Jadidi et al. [21] explored spectra-derived features through color space and coordinate transforms, from images belonging to three sources while working on data collected from Central European lakes. The combined use of spectral–spatial features obtained through connected morphological operators has also been investigated [22]. One of the most extensive studies is due to Neil et al. [23], who explored 185 inland and coastal aquatic systems at a global level, through which the performances of 48 distinct Chl-a estimation algorithms were tested.

However, the Chl-a estimation performance of semiempirical approaches depends heavily on the effectiveness of the underlying machine learning method. Consequently, the paradigm-shift that has occurred in the field of machine learning through the advent of deep learning (DL) [24] is of paramount importance in this context. More specifically, DL has led to groundbreaking performances in various highly challenging computer vision tasks, especially via convolutional neural networks (CNNs) [25]. This new family of algorithms enabling the design and efficient training of deep neural networks has rendered, de facto, the process of feature engineering redundant, as it is now delegated to the network itself. Depending on the quality of the provided data, the deep networks often compute features outperforming handcrafted alternatives. For a survey of DL applications in remote sensing, the reader is referred to [26].

It is thus not surprising that DL methods have been already applied to Chl-a estimation. In particular, Peterson et al. [27] have developed an artificial neural network (ANN) with six hidden layers, operating on spectral signatures of Landsat-8 and Sentinel-2 images at 30 m resolution, and applied it to the water quality estimation of lakes in the United States. It constitutes one of the first studies that assesses the suitability of DL for this context, and their results show increased accuracy and robustness with respect to alternative methods. In addition, Pu et al. [28] reported the first application of CNNs to this problem, where small patches of Landsat-8 images, centered on the pixel under study, were used as network input, with the end goal of classifying the water quality of Chinese lakes. Since, however, Chl-a estimation is inherently a regression problem, such networks estimating the direct Chl-a level of its input have appeared as well [29,30]. Moreover, in an effort to exploit temporal correlations across multitemporal data, long short-term memory networks have also been studied, both individually [31] and together with CNNs [32]. Lastly, a highly comprehensive study was presented in [8], on lakes across four continents, using a relatively large dataset of 2943 samples. They proposed a single cross-mission solution for both Sentinel-2 and Sentinel-3 sensors, based upon a five-layer mixture density network, and obtained remarkable results.

1.2. Aim and Contribution

As with many remote sensing applications (e.g., semantic segmentation), overfitting due to labeled data shortage is a critical problem in this context as well. In fact, even more so, as, for each Chl-a sample of the ground truth, a rigorous in situ sample collection and laboratory analysis workflow is required (as opposed to an expert labeling pixels on a screen), and as such, labeled datasets are often relatively small. This problem is further exacerbated by the notorious labeled data need of contemporary deep networks. Consequently, it is imperative to exploit, to the maximum extent possible, whatever precious little amount of data is available.

To this end, this article’s main contribution is a new network structure that learns to estimate Chl-a concentration levels while exploiting both labeled image locations, as well as unlabeled portions of the same visual content. We propose to achieve this via a multitask double-branch CNN, whose main task is Chl-a level regression via training with labeled input, and whose auxiliary task is the classification of unlabeled samples to its month of acquisition. The proposed approach was tested with multispectral Sentinel-2 images of Lake Balik in North Turkey, using samples collected over 3 years. The ablation studies conducted show that the inclusion of unlabeled data improves the correlation coefficient of the resulting model by multiple percentile points.

In the sequel of this article, following an overview of the collected dataset and study area (Section 2.1), we elaborate on the proposed method of Chl-a concentration estimation (Section 2.2), present the results of our regression experiments (Section 3), and discuss our findings, while Section 4 is devoted to concluding remarks.

2. Materials and Methods

2.1. Study Area and Data

The study area for our water quality experiments is Lake Balik, located in the northern part of Turkey, at the Kizilirmak Delta (Figure 1a). It has an area of 13.9 ${\mathrm{km}}^2$ and an average depth of 1.5 m. Maximum depth reaches 2.5 m in summer and 3.5 m in winter. Lake Balik has been a designated Ramsar site since 1998 and is home to 74% of the bird species encountered in Turkey. It is surrounded by agricultural areas, and unfortunately is under significant threat from uncontrolled fertilizer and pesticide use. Surface vegetation is present only at the shallows and extends to, at most, 3 m from the shore [33].

In situ measurements were conducted on 32 distinct dates from April 2017 until June 2019 (7 during 2017, 14 during 2018, and 11 during 2019), once, or at most twice, a month, depending on weather conditions. At each field visit, samples were collected from 10 measurement sites, approximately 500 m apart and 250–300 m from the shore (shown in Figure 1a), with no surface vegetation and minimal water movement. Sample collection was arranged to coincide with the satellites’ visits over the lake, and lasted at most 3 h. Thus, the resulting dataset contains $10\times 32=320$ samples.

Sample analysis was conducted according to EPA-445 [34]. In particular, 1-liter water samples were collected from a depth of 0.2–0.3 m below the lake surface and filtered on 0.7 $\mathsf{\mu }$m Whatman GF/F filter pads, and then stored in the dark at −20 ${}^{\circ }$C pending laboratory analysis. Filtering was performed under subdued light, soon after sampling, since algal populations, and thus Chl-a concentration, are susceptible to change in relatively short periods of time. Chl-a was extracted using 6 mL of $90\%$ acetone for 24 h [34]. The Chl-a concentration was corrected for pheophytin and measured fluorometrically using a TD-7210 Fluorometer (Turner Designs, CA, USA) fitted with a daylight white lamp and properly calibrated chlorophyll optical kit.

Laboratory analysis of the samples shows that the lake is in hypereutrophic state [35], with an annual average Chl-a concentration of 35.82 $\mathsf{\mu }$g/L; the distribution of concentrations across time and measurements sites are shown in Figure 1b.

Atmospheric correction was conducted on the Level 1C Sentinel-2 images with the Polymer method v.4.6 [36], using default parameters. Polymer is an atmospheric correction algorithm for processing oceanic waters with and without sun-glint presence [37], and has been reported to show good performance [8]. Finally, all bands were resampled to 10 m/pixel spatial resolution via bilinear interpolation. After removing noisy pixels close to the shore and the in situ sampling locations, approximately 74 thousand lake surface pixels remain (Figure 2c) for which Chl-a estimation values are unknown, forming an unlabeled dataset of approximately 2.3 million samples across all acquisition dates (

Table 1. The number of samples per month for labeled and unlabeled datasets.

	Month
	1	2	3	4	5	6	7	8	9	10	11	12	Total
Unlabeled	148 K	222 K	222 K	222 K	222 K	296 K	222 K	148 K	148 K	296 K	148 K	74 K	2365 K
Labeled	20	30	30	30	30	40	30	20	20	40	20	10	320

).

2.2. Proposed Method

Chl-a levels of inland waters are known to be positively correlated with water temperature and fertilizer inflow often resulting from rainfall, among other meteorological parameters [38]. In other words, one can expect a certain level of correlation between Chl-a levels and the months of a year, as regular seasonal fluctuations occur (i.e., higher concentrations during sunnier seasons and after rains). Based on this observation, we hypothesize that a network trained to guess the month of a given remote sensing spectral or spectral/spatial sample (with no associated Chl-a information/label) will be able to indirectly learn features helpful for Chl-a level estimation.

To avoid any confusion, we assume that the date of acquisition is always available for the remote sensing images under study, and in the rest of this article "labeled dataset" will refer to the limited image samples of the water body for which the Chl-a concentration levels are known, and “unlabeled dataset” will refer to the rest of the image samples of the water body.

The aforementioned concepts can be implemented in various ways. For instance, one could first train a convolutional network with unlabeled image locations, to learn the month of acquisition until a reasonable level of performance, and then switch the loss function to regression with a smaller learning rate and continue training with the labeled samples as a form of “forward” transfer learning. However, despite its simplicity, this approach risks directing the network to search space regions too far from the optima of the main task of Chl-a estimation, to which the network is only afterwards introduced.

Instead, a superior alternative is to train the network’s feature extractor for both tasks simultaneously as a multitasking network (Figure 3), where the main task is Chl-a regression and the auxiliary task is month classification. The two tasks form the two branches following the feature extraction layers, and direct the network to learn features useful for both tasks, thus reinforcing each other. Moreover, being able to benefit from the often far greater amount of unlabeled data in this way can have a strong regularization effect, thus avoiding overfitting.

2.2.1. Input Shape

Since the availability of spatial information has been shown to enable networks to exploit correlations among neighboring pixels [28,29], it was chosen to provide the input as a 3D tensor of the shape $b\times k\times k$. It is formed by stacking b multispectral image patches of size $k\times k$ pixels centered on the pixel under study.

2.2.2. Network Architecture

We adapted the 2D CNN architecture presented in [29] as the feature extractor shared across tasks. It consists of four convolutional layers without pooling, as the input patches are already small, and two fully connected (FC) layers. Before the convolutional layers, 1-pixel wide padding is applied to be able to use them with 3 × 3 kernel sizes. Each convolutional layer is followed by a $tanh$ activation function. Then, the feature map generated by the 4th convolutional layer is fed into each of the two task-specific branches, each consisting of two FC layers. These task-specific branches are identical except for the number of neurons in the second FC layer. In this layer, the regression branch has 1 neuron to estimate the continuous Chl-a values, while the classification branch has 12 neurons to predict the sample acquisition month. The number of shared layers (four) and the learning rate of each branch were determined through ablation studies (presented in the next section). The overall network architecture is shown in Figure 3.

2.2.3. Loss Functions

More formally, consider a Chl-a estimation dataset ${\mathcal{D}}_L={\left\{(x_i,y_i,m_i)\right\}}_{i=1,\dots ,|{\mathcal{D}}_L|}$, called the labeled dataset, where each patch $x_i$ has Chl-a concentration level $y_i$ and is labeled with the acquisition month $m_i$ for $x_i\in {\mathbb{R}}^{b\times k\times k}$, $y_i\in \mathbb{R}$, and $m_i\in \{1,2,\dots 12\}$. We also have an unlabeled dataset ${\mathcal{D}}_U={\left\{(x_i,m_i)\right\}}_{i=1,\dots ,|{\mathcal{D}}_U|}$ where each patch $x_i$ is labeled only with the acquisition month $m_i$. The labeled dataset ${\mathcal{D}}_L$ is further split into a training set $\mathcal{T}$ and a test set $\mathcal{H}$ with $\left|\mathcal{T}\right|$ and $\left|\mathcal{H}\right|$ samples, respectively, where $\left|\mathcal{T}\right|+\left|\mathcal{H}\right|=|{\mathcal{D}}_L|$.

The goal is to learn a model ${\Phi }^{\theta }$ with parameters $\theta$ to correctly estimate the Chl-a concentration of a given sample x as $y={\Phi }^{\theta }\left(x\right)$. Contrary to existing approaches [28,29,30,31,32,39] that rely on a single network architecture for the task of Chl-a retrieval, our approach consists of a feature extractor ${\Phi }_{fe}^{\theta }$ shared between two task-specific branches, namely ${\Phi }_{reg}^{\theta ,FC}$ and ${\Phi }_{cls}^{\theta ,FC}$, such that the model produces an estimation for the main task, ${\Phi }_{reg}^{\theta }=({\Phi }_{fe}^{\theta },{\Phi }_{reg}^{\theta ,FC})$, and a prediction for the auxiliary task, ${\Phi }_{cls}^{\theta }=({\Phi }_{fe},{\Phi }_{cls}^{\theta ,FC})$.

MSE loss: Following the standard approach in Chl-a estimation studies [29,30,32], the MSE (mean squared error) loss is used for the main regression task. The loss is evaluated only for the training set $\mathcal{T}$ of the labeled dataset ${\mathcal{D}}_L$,

$${\mathcal{L}}_{Reg}(x_{\left(i\right)},y_{\left(i\right)};{\Phi }_{reg}^{\theta })=\frac{1}{|{\mathcal{B}}_{reg}|}\sum _i{(y_{\left(i\right)}-{\Phi }_{reg}^{\theta }\left(x_{\left(i\right)}\right))}^2$$

(1)

where $\left(i\right)$ represents the ith batch of the training set $\mathcal{T}$, $|{\mathcal{B}}_{reg}|$ is the batch size of the regression task, and ${\Phi }_{reg}^{\theta }\left(x_{\left(i\right)}\right)$ and $y_{\left(i\right)}$ denote, respectively, the estimated and target Chl-a concentrations.

Cross-entropy loss: As far as the auxiliary task is concerned, a “standard” classification loss is used to exploit the unlabeled image samples by predicting their acquisition months. It is evaluated for both the training set $\mathcal{T}$ of the labeled dataset ${\mathcal{D}}_L$ and the unlabeled dataset ${\mathcal{D}}_U$, since the samples of both datasets are labeled with the acquisition month:

$${\mathcal{L}}_{Cls}(x_{\left(i\right)},m_{\left(i\right)};{\Phi }_{cls}^{\theta })=-\frac{1}{|{\mathcal{B}}_{cls}|}\sum _i{m}_{\left(i\right)}log\left({\Phi }_{cls}^{\theta }\left(x_{\left(i\right)}\right)\right)$$

(2)

where $\mathcal{D}=\mathcal{T}\cup {\mathcal{D}}_{\mathcal{U}}$, and $(x_{\left(i\right)},m_{\left(i\right)})\in \mathcal{D}$, $\left(i\right)$ represents the ith batch of $\mathcal{D}$, $|{\mathcal{B}}_{cls}|$ is the batch size of the classification task, and ${\Phi }_{cls}^{\theta }\left(x_{\left(i\right)}\right)$ and $m_{\left(i\right)}$ denote, respectively, the predicted and target months.

2.2.4. Multitask Loss

The design of the multitask loss function is of critical importance. A naive multitask loss ${\mathcal{L}}_{MT}$ could be written as a simple addition of the respective task losses, ${\mathcal{L}}_{MT}={\mathcal{L}}_{Reg}+{\mathcal{L}}_{Cls}$. However, the range of values produced by ${\mathcal{L}}_{Reg}$ and ${\mathcal{L}}_{Cls}$ can vary strongly, especially considering that they belong to different loss families, thus possibly leading the network to ignore one of the tasks. A better alternative is instead to express the total network loss ${\mathcal{L}}_{MT}$ as a weighted combination of the task-specific losses ${\mathcal{L}}_{MT}=c_{Reg}{\mathcal{L}}_{Reg}+c_{Cls}{\mathcal{L}}_{Cls}$. It is observed in [40] that by adjusting appropriate values for $c_{Reg}$ and $c_{Cls}$, the optimization of ${\mathcal{L}}_{MT}$ ensures that no task is left behind when the other one is being learned. This can be achieved by manually computing the coefficient $c_{\tau }$ of each task ${\mathcal{L}}_{\tau }$, where $\tau \in \mathrm{T}$ and $\mathrm{T}=\{Reg,Cls\}$, but this is a tedious and error-prone process that also does not take into account changes in loss during the training phase. Instead, these loss coefficients $c_{\tau }$ can be learned by the model by adding them to the trainable network parameters $\theta$ as ${\omega }_{\tau }=(\theta ,c_{\tau })$. To define the multitask objective function with the aforementioned loss coefficients, we used an implementation (https://github.com/Mikoto10032/AutomaticWeightedLoss, accessed on 16 December 2021) of Liebel and Körner’s approach [41]:

$${\mathcal{L}}_{MT}(x,y_{\tau };{\omega }_{\tau })=\sum _{\tau \in \mathrm{T}}\frac{1}{2c_{\tau }^2}{\mathcal{L}}_{\tau }(x,y_{\tau };{\omega }_{\tau })+ln(1+c_{\tau }^2)$$

(3)

where $\tau \in \mathrm{T}$ and $\mathrm{T}=\{Reg,Cls\}$, $y_{\tau }=y$ for $\tau =Reg$ and $y_{\tau }=m$ for $\tau =Cls$.

In the next section, we report our numerical results and detail the effects of multitask learning, patch size, length of the shared feature extractor, number of unlabeled samples, and task loss weighting scheme on the Chl-a estimation performance through ablation studies.

3. Results and Discussion

The experiments presented in this section aimed to test our hypothesis of whether unlabeled data can contribute to the Chl-a estimation performance via multitask learning.

3.1. Data Split

Contrary to the state-of-the-art, where single random splits are common [16,21,27,31,32,42], a 10-fold cross-validation strategy was adopted to avoid overfitting. Subsequently, the training data is further split in order to produce a validation set with the same size as the test set to be used for hyperparameter tuning. Consequently, in each fold we have 256 samples (80%) for training, 32 samples (10%) for validation, and 32 samples (10%) for testing. After the hyperparameters are tuned, the model is trained without generating the validation set, since the size of the dataset is relatively small. Thus, the final split has 288 samples (90%) for training and 32 samples (10%) for testing.

3.2. Setup

The proposed multitask learning model was compared against several alternative approaches from the state-of-the-art, including its single-task regression version [29], support vector regression (SVR) [15] as a traditional machine learning method using a radial basis function kernel, multilayer perceptrons (MLP) with various numbers of hidden layers [30] as nonconvolutional neural networks, where the total number of trainable parameters was adjusted to be equal to the proposed model, a recently introduced morphology-based spectral–spatial approach (EFAL) [22], a convolutional neural network (CNN), as proposed in [28] adapted for regression, and of course a statistical spectral band characteristic known to be effective for Chl-a description [16] (band ratio): $R_{705}/R_{665},(R_{665}^{-1}-R_{705}^{-1})\times R_{740},R_{490}/R_{443}$.

3.3. Training

Since the sizes of the labeled and unlabeled training sets differ by several orders of magnitude, it was not a plausible approach to set the batch sizes to have the same number of iterations. Instead, batches of the labeled training set are reused to ensure the same number of iterations on both tasks. The learning rate is chosen to be ${10}^{-4}$. To optimize the multitask network, the root mean square propagation (RMSProp) algorithm is used. Training lasts for 200 epochs as determined through experimentation. The batch size is chosen as 64 for the labeled set and as 12,318 for the unlabeled set. The same parameters/hyperparameters determined by the validation set were used in all experiments to ensure consistency. Finally, the coefficient of determination ($R^2$), root mean square error (RMSE), and mean absolute error (MAE) were used for performance measurement.

3.4. Results

According to the results presented in

Table 2. Mean Chl-a estimation results across 10 folds of cross-validation with the best results in bold, where $R^2$, MAE, and RMSE denote the coefficient of determination, mean absolute error, and root mean square error, respectively.

Method	Dataset	${\mathit{R}}^2\times 100$	MAE	RMSE	Model Parameters
Band ratio [16]	Labeled	71.37	0.147	0.199	-
SVR [15]	Labeled	51.82	0.193	0.27	-
EFAL [22]	Labeled	77.12	0.151	0.221	-
CNN [28]	Labeled	73.06	0.146	0.197	340 K
MLP (2 layers) [30]	Labeled	64.99	0.176	0.233	267 K
MLP (3 layers) [30]	Labeled	72.47	0.152	0.204	267 K
MLP (4 layers) [30]	Labeled	75.18	0.139	0.194	267 K
MLP (5 layers) [30]	Labeled	72.11	0.144	0.206	267 K
Single-task CNN [29]	Labeled	76.09	0.142	0.19	192 K
Multitask CNN (ours)	Labeled+unlabeled	83.89	0.107	0.154	267 K

, the band ratio approach [16] provides results close to those of the DL models, as it was expertly derived to obtain Chl-a values in water bodies from Sentinel-2 images. As expected, all deep learning models outperformed SVR [15], confirming the reported superiority of deep learning with respect to classic machine learning. EFAL [22] provides the second-best result thanks to its ability to represent pixels with a hierarchy of connected components, which makes it suitable for Chl-a estimation, since it provides a multiscale pixel description and adaptability in the spatial domain. CNN [28] performs similarly to MLP and single-task CNN models in that its performance is limited with the samples of the labeled dataset. For the MLP models [30], even though they have approximately the same number of trainable parameters with the proposed architecture, the estimation performance increases up to four layers and then decreases, possibly due to the model starting to overfit. The multitask learning model trained on both labeled and unlabeled data outperforms the single-task regression approach [29]. Despite the lack of Chl-a measurements in the unlabeled data, capturing seasonal appearances through the auxiliary task’s month classification appears to help Chl-a estimation performance by almost 8 percentile points.

In addition to the quantitative results, Chl-a estimation maps were generated via the proposed method for four dates in 2017 (one from each season) by training with the remaining years’ data (Figure 4).

3.5. Ablation Study

Furthermore, the effects on performance of crucial design parameters were studied and are reported in this section: the number of shared layers within the architecture, the size of the input image patch, the number of unlabeled samples, and the task loss weighting scheme.

Since the number of layers in the proposed architecture is six, we experimented with the number of shared layers ranging from one to five by training five multitask models, starting with the first layer. The results obtained are reported in

Table 3. Chl-a estimation levels of the proposed multitask model across 10-fold cross-validation with different number of shared layers among the two task branches, with the best results in bold, where $R^2$, MAE, and RMSE denote the coefficient of determination, mean absolute error, and root mean square error, respectively (patch size = $3\times 3$).

Shared Layers	1	2	3	4	5
$R^2\times 100$	83.57	82.81	83.38	83.89	83
RMSE	0.156	0.159	0.158	0.154	0.16
MAE	0.108	0.11	0.109	0.107	0.113

. The best performance is obtained by splitting the architecture into two task-specific branches after the fourth layer. This suggests that the multitask model benefits from sharing all of its convolutional features between the main regression task and the auxiliary classification task.

For the sample patch size, experiments were conducted with sizes ranging from 1 × 1 to 9 × 9 pixels. Based on the results of

Table 4. Chl-a estimation levels of the proposed multitask model across 10-fold cross-validation with various patch sizes, with the best results in bold, where $R^2$, MAE, and RMSE denote coefficient of determination, mean absolute error, and root mean square error, respectively (number of shared layers = 4).

Patch Size	1 × 1	3 × 3	5 × 5	7 × 7	9 × 9
$R^2\times 100$	75.44	83.89	81.48	82.19	83.45
RMSE	0.194	0.154	0.167	0.163	0.158
MAE	0.134	0.107	0.114	0.112	0.105

, one can observe the contribution of spatial information with respect to pixelwise analysis. However, as the spatial resolution is low, larger patch sizes hardly improve it.

In addition, to measure the effect of the number of unlabeled samples on performance, four multitask learning models were trained with progressively more unlabeled samples, ranging from none all the way to the entire lake surface. The intermediate lake masks with the employed unlabeled image locations are shown in Figure 2a–c, and the results are reported in

Table 5. Chl-a estimation results of the proposed multitask model across 10-fold cross-validation with different number of unlabeled samples used for its auxiliary task, with the best results in bold, where $R^2$, MAE, and RMSE denote coefficient of determination, mean absolute error, and root mean square error, respectively.

Unlabeled Samples	0	28160	503040	2364928
$R^2\times 100$	76.63	82.41	83.03	83.89
RMSE	0.187	0.161	0.159	0.154
MAE	0.137	0.118	0.118	0.107

. As expected, training the model with the utmost possible number of unlabeled samples by using the entire lake surface achieves the best performance. This is a promising result for lakes that are geographically difficult to access, and it obtains a large number of in situ samples.

The effect of the number of unlabeled samples on the auxiliary classification task was also examined. Similarly to the previous ablation for Chl-a estimation, four models were trained with an increasing number of unlabeled samples using the same lake intermediate masks as in Figure 2a–c. Likewise, the multitask learning model trained using all possible samples achieved the best performance, and their confusion matrices are shown in Figure 5.

Finally, a couple of experiments focusing on task-loss weighting were conducted, where a direct sum of task-specific losses was compared against a task-weighting scheme learned by the model. According to

Table 6. Chl-a estimation scores of the proposed multitask model across 10-fold cross-validation with two task loss weighting schemes: direct sum and the weighting learned by the model during training, with the best results in bold.

Task Loss Weight	Sum	Weight Learning [41]
$R^2\times 100$	83.45	83.89
RMSE	0.156	0.154
MAE	0.113	0.107

, the latter is superior, but with only a small margin.

3.6. Discussion and Limitations

The results presented so far provide relatively clear indications of the proposed method’s practical interest with respect to its counterparts. The multitask network outperformed its single-task-regression-oriented alternative with a significant margin, and exhibited increasing performance as the number of unlabeled samples increased. Moreover, one can also safely presume that the large unlabeled dataset probably has an additional regularization effect. In other words, using such a large dataset helps the model learn a better representation of the data. Thus, the auxiliary task satisfies the notorious data hunger of the deep learning model. However, neither the method nor the undertaken study are without limitations.

First of all, even though utmost care was taken to ensure the reliability of the estimation results and to avoid overfitting, this still does not change the fact that all of the samples originate from a single water body. Hence, it is unknown whether the estimation performance of the method under study can be reproduced with samples of a geographically distinct lake.

Furthermore, even though the samples acquired from Lake Balik demonstrated a significant deviation (from 11.04 $\mathsf{\mu }$g/L to 108.35 $\mathsf{\mu }$g/L), the lake is still in overall hypereutrophic state. Consequently, the validity of the results obtained, across distinct mesotrophic or oligotrophic lakes, remains uncertain.

It is planned to rectify both limitations in the future, via studies in new lakes, pending financial support.

4. Conclusions

In this article, we proposed the first neural network architecture for water quality estimation that can exploit unlabeled satellite imagery. We addressed the scarcity of labeled data with a multitask CNN, where the main task is Chl-a regression and the auxiliary task is month classification for the unlabeled samples. We thus aimed to exploit the semiregular seasonal fluctuations of Chl-a in order to improve the computed features.

The experiments that were conducted on Lake Balik showed that the proposed multitask learning model increased the coefficient of determination by 7.8 percentile points compared to the single-task model, asserting the contribution potential of unlabeled data via the proposed technique. In fact, the progressive improvement of performance relative to the amount of employed unlabeled data was also shown.

Consequently, the proposed method renders possible the development of a Chl-a estimation model with a modest amount of labeled data, thus minimizing the effort required in terms of the number in situ samples, while enabling the exploitation of unlabeled remote sensing images that are abundantly available.

Moreover, since the only assumption was the semiregular seasonal fluctuations of Chl-a levels, the auxiliary month classification task, in our opinion, holds the potential of benefiting other multitemporal remote sensing analysis contexts with similarly behaving variables.

Nevertheless, despite the strong potential of these findings, and the rigorous experimental protocol to avoid overfitting, this study is limited all the same by the use of a single water body. Future work will concentrate on the validation of this approach with water bodies of distinct trophic indexes at alternative geographical locations. Furthermore, we will also focus on how to exploit multitemporal unlabeled data via month classification for land-cover/land-use map estimation.