1.

Introduction

Chlorophyll-a (Chla) concentration is a key indicator of the biophysical status of water bodies. Numerous satellites programs operated by the National Aeronautics and Space Administration (NASA), the National Oceanic and Atmospheric Administration (NOAA), and the European Space Agency (ESA) have global ocean monitoring capability (Table 1). Accurate and fast estimation of global Chla concentration from these satellites’ images is crucial for various purposes, such as the net primary production derivation,^1,2 harmful algal blooms detection,^3–6 physical and biological interactions,^7,8 and ecosystem models validation.^9,10

Table 1

Summary of primary ocean color satellites.

The twin-moderate resolution imaging spectroradiometer (MODIS), terra-MODIS and aqua-MODIS, have high temporal resolution to view the entire Earth’s surface every 1 to 2 days which have provided big remote sensing data. The big data characteristics of satellite remote sensing calls for an efficient Chla concentration estimation algorithm that can handle the big remote sensing data in an efficient and effective manner. Nevertheless, the relationship between the Chla concentration and the remote sensing observations is a very complex nonlinear function, as a consequence, design of an end-to-end machine learning approach with strong learning capability that can capture the complex nonlinear relationship between the Chla concentration and the remote sensing observations is an important research issue. Support vector machine, a classical traditional machine learning approaches, has limitations due to their relative weak modeling capability compared to neural networks to handle the complexity and uncertainty in the inverse problem and also the low computational efficiency.¹¹ Deep learning, especially convolutional neural networks (CNNs), is well known for its using both spatial and spectral information, strong modeling capability and leveraging GPU for high computational efficiency, by greatly outperforming the other machine learning techniques in a lot of computer vision tasks.^12–16 Remote sensing scientists have seen the value of CNNs when they face the tasks of classification, segmentation, object detection, change detection, and super resolution.^17–19 However, very limited studies focus on using CNNs to solve regression problems. Although CNNs have huge potential for addressing the problems in Chla estimation, to the best of our knowledge, there is no study to explore the capability of CNNs for the estimation of global Chla concentration. Therefore, this paper explores deep CNN for the estimation of global Chla concentration. Nevertheless, deep CNN requires a lot of training samples, but the ground-truth field data are very limited, and as such we adopt the Ocean Color Climate Change Initiative (OC-CCI) images as ground truth.

The OC-CCI project is one of the 13 projects in the ESA CCI program that aim at addressing a particular essential climate variable from satellite observations. The objective of this project is to produce a long-term, multisensor time-series of satellite ocean color data, e.g., Chla concentration, with specific information on errors and uncertainties.²⁰ The OC-CCI Chla concentration dataset is produced by combining the ocean chlorophyll 3 (OC3) algorithm, the ocean chlorophyll 5 (OC5) algorithm, and the OCI algorithm [ocean chlorophyll 4 (OC4) algorithm + color index (CI) algorithm]^21–23 on the merged data of SeaWiFS, MERIS, MODIS, and VIIRS. Implemented by a sophisticated processing chain, OC-CCI intends to provide a global Chla concentration product of the highest quality, particularly on Chla retrievals from case 2 waters in global scale, although these may not be the latest.²⁴ Although inevitable biases exist in the OC-CCI products when being matched against the in situ observations, they are the most accessible and reliable Chla concentration product.

This paper explores the use of the CNN regression method for the end-to-end estimation of the global Chla concentration from the imagery recorded by the MODIS sensor on board the Aqua satellite. Unlike approaches that build the relationship between reflectance and Chla concentration obtained from the SeaBASS dataset or cruise measurements,^25,26 the CNN method here derives the relationship between reflectance and Chla concentration acquired from the OC-CCI Chla image, which can provide abundant training samples for effective training of the CNN model. MODIS remote sensing reflectance (Rrs) images and OC-CCI monthly Chla concentration images are used as input and ground truth, respectively. Nearly 15,000,000 samples are selected for training the patch-based CNN regression model. A total of 12 monthly global Chla concentration images are produced by the CNN and support vector regression (SVR), respectively. Qualitative and quantitative analyses are used for evaluation and comparison analysis, and the results imply that the CNN model can be successfully used to estimate the global Chla concentration. This effort represents an initial attempt to estimate Chla concentration using CNN.

2.

CNN

A CNN usually takes an image as input and outputs its corresponding class label. CNN architectures comprise typically several layers: convolutional layers, pooling layers, and fully connected layers. A slightly different version from the classical LeNet¹² architecture is shown in Fig. 1 to explain CNN.

Fig. 1

Structure of a CNN which is slightly different from the classical LeNet. The order of 3-D $X$ and $a$ are channel × height × width. Kernels $W$ in each convolution layer are four-dimensional and the number before at stands for number of the 3-D kernels, so the order of the four-dimensional kernels $W$ is number × channel × height × width.

The convolution layer 1 takes an RGB image $X$ and kernels $W^{(1)}$ which consist of weights as inputs and outputs feature maps $a^{(1)}$. The size of $X$ is $\text{channel}\times \text{height}\times \text{width}=3\times 28\times 28$ and the size of $W^{(1)}$ is $\text{number}\times \text{channel}\times \text{height}\times \text{width}=20\times 3\times 5\times 5$ which means that there are 20 three-dimensional (3-D) kernels. Kernels should have the same channels as the corresponding input image. In general, there will be multiple filters and multiple feature maps are generated and stacked to form a multilayer feature map.^27,28 Each feature map is obtained by dot product between receptive field and kernel in specific rules of scanning the input image. The weights in the kernel will be learned by training. The process of convolutional layer 1 can be expressed as

Conventionally, each convolutional layer is followed by a pooling layer in order to reduce the number of parameters to learn and reduce the variance of features,^29–31 which can be accomplished by operations like maximum, average, etc. No weights need to be learned in this layer. The process of pooling layer 1 can be expressed as

A fully connected layer takes all neurons in the input layer, no matter a convolutional layer or fully connected layer, and connects them to each neuron in its layer. It is regularly used as the last few layers in a CNN architecture and yields as many outputs as labels to a classifier layer. The process of two fully connected layers can be expressed as

3.

Methodology

In this paper, a patch-based CNN regression method is used to estimate global Chla concentration from the MODIS imagery. The patch-based regression CNN model scans through images and crops a patch from the images at each location of the scanning. For each patch, it uses the regression CNN to predict the Chla concentration value of the pixel at the patch center. Then the model assembles the whole predicted values to form the global Chla concentration map. The complete procedure is shown in Fig. 2. It consists of three parts: the preprocessing of the input data, the training of the CNN model, and the data prediction employing the trained CNN, which has the same procedure as testing. Performance of the CNN method is evaluated on the whole dataset.

Fig. 2

Flowchart of the proposed algorithm.

3.3.

Structure of the Patch-Based Regression CNN

The structure of the patch-based regression CNN is demonstrated in Fig. 3. The main difference from Fig. 1 is the last layer which is specifically designed by outputting only one value to solve regression problem. The CNN consists of two convolutional layers, two pooling layers and two fully connected layers. Strategy of how to tune the parameters to get the structure of the CNN will be explained in Sec. 5.

Fig. 3

Structure of the patch-based regression CNN. The order of 3-D $X$ and $a$ are channel × height × width. Kernels $W$ in each convolution layer are four-dimensional and the number before at stands for number of the 3-D kernels, so the order of the four-dimensional kernels $W$ is number × channel × height × width. The red $T^n$ is the residual error in the $n$’th layer during backpropagation which has the same size as the corresponding $a$.

The input $X$ is image patch with the size of $4\times 15\times 15$. The output of the convolution layer 1 is 3-D feature map $a^{(1)}$ with the size of $20\times 13\times 13$, which is a result of the convolution of input patch with the kernel $W^{(1)}$ with stride 1. A max-pooling layer is followed to reduce the size to $20\times 6\times 6$, and the maximum operation is chosen in the pooling layer owing to its simplicity and effectiveness. Another convolution of kernel $W^{(2)}$ with stride 1, following a max-pooling of process on $a^{(2)}$, yielding feature map $a^{(3)}$ with the size of $50\times 4\times 4$ and $a^{(4)}$ with the size of $50\times 2\times 2$. The first fully connected layer takes input of dimension $(50\times 2\times 2)=200$ and outputs 100 neurons. Unlike employing activation function sigmoids in the classical LeNet, the activation function ReLU, $\max \{0,z\}$, is used in each convolutional layer and the first fully connected layer because it helps accelerate convergence in training.⁴⁴ Different with the classification problem that the output number of the last fully connected layer is the number of the class, to address the regression here, the last fully connected layer takes as input the 100 neurons and yields one neuron which corresponds to the Chla concentration value.

4.

Results

Using the patch-based CNN method, 12 monthly global Chla concentration images are predicted, and the prediction of each image takes about 2 h. A one-month Chla concentration map estimated from the Rrs images by the CNN model and the SVR model, as well as the corresponding OC-CCI ground-truth image are shown in Fig. 3. For both SVR and the CNN predictions, values below 0.001 or above $100\mathrm{mg}{\mathrm{m}}^{-3}$ (representing $<0.001\%$ of the Chla estimates) are removed to limit the product as same as the OC-CCI ground-truth image. From Figs. 4(a) and 4(b), we can see that the prediction map achieved by the CNN model has a strong spatial consistency with the OC-CCI ground-truth Chla concentration map, which is achieved by assimilating various sources of information, such as the proximity to land, the depth of the ocean, and ocean currents.⁵¹ Figure 6(a) shows the scatter plots of the predicted values achieved by the CNN model and ground-truth OC-CCI Chla concentration values. In a scatter plot, good performance means that points are close to the equal-value line. Figure 6(a) shows that the majority of data points distribute around the equal-value line, indicating that the CNN method owns the capacity to predict very accurately the Chla concentration information in OC-CCI data, by learning the complex relationship between the Rrs spectra data and the OC-CCI data.

Fig. 4

Spatial distribution of the monthly global Chla concentration in March 2016 (a) obtained using the CNN method, (b) the OC-CCI data, and (c) obtained using the SVR method. In general, Figs. 4(a) and 4(c) have strong spatial consistency with Fig. 4(b). But in the three red rectangle areas, Fig. 4(c) is much deeper red than Fig. 4(a), which means that heavier overestimation happened in the prediction map achieved by the SVR method than the CNN model.

In Table 4, we can see that the $R^2$ achieved by the CNN model is between 0.879 and 0.918 with the mean of 0.901, indicating a high consistency between the Chla concentration values predicted by the CNN model and the ground-truth OC-CCI values. The other numerical measures also indicate small error and high accuracy of the proposed CNN model, i.e., the overall RMSE of the estimated Chla concentration is 0.129, and the max bias and MAPE are 0.054 and 26.687, respectively. The max and min of MAE are 0.153 and 0.102, respectively, and the mean value is 0.125. We want to highlight the fact that although the proposed CNN model is only trained on the month 1 data, it achieves stably good performance on the other 11 months’ data, indicating a strong generalization capability of the CNN model.

Table 4

Statistical results for 12-month Chla concentration estimated by the CNN model.

Figure 7 shows the RMSE and $R^2$ performance comparison between the CNN model and the SVR method on 12 months Chla concentration estimation. We can see that the $R^2$ curve achieved by the CNN model is above the curve of the SVR method, indicating that CNN can outperform SVR over all the 12-month predictions. Moreover, the RMSE achieved by CNN are also consistently lower than SVR. The mean RMSE and $R^2$ of the 12 months are (0.129, 0.901) and (0.175, 0.828) for CNN and SVR, respectively. Figures 4(c) and 5 shows that the color of some case 2 water areas are much deeper red than Figs. 4(a) and 5, i.e., Yellow Sea and East China Sea (red rectangle 1), English Channel and the south of North Sea (red rectangle 2), and the east of Bering Sea (red rectangle 3), which means that heavier overestimation happened in the prediction map achieved by the SVR method than the CNN model. Figure 6 supports the visual analysis above, i.e., although both the CNN model and the SVR method overestimate in the high Chla concentration areas, the SVR method performs much worse. From Fig. 6, we can also see that the CNN model outperforms the SVR method in the low Chla concentration areas.

Fig. 5

Zoomed in regions of the corresponding red rectangle areas in Fig. 4. (a) CNN, (b), OC-CCI, and (c) SVR.

Fig. 6

The scatter plot between (a) the Chla concentration values estimated by the CNN model and the ground-truth OC-CCI Chla concentration values. (b) The values achieved by the SVR method and the OC-CCI values. The SVR approach tends to overestimate when the Chla concentration is high, which is consistent with Figs. 4 and 5 that the SVR results have higher values than the ground truth in high Chla concentration areas.

Fig. 7

The performance measures comparison of CNN and SVR for 12 months Chla concentration estimation based on RMSE and $R^2$.

5.

Discussion

Based on the above experimental results, we can summarize that the CNN approach can better learn the complex relationship between the Rrs and the OC-CCI Chla values than SVR, demonstrating the possibility that the CNN may be used for building an end-to-end approach for efficient Chla concentration estimation. Song et al.⁵² proposed a novel inversion method for rough surface parameters using CNN, which microwave images of rough surfaces are used for training. In Ref. 53, a CNN-based method is employed to establish the solar flare forecasting model, in which forecasting patterns can be learned from line-of-sight magnetograms of solar active regions. Line-of-sight magnetograms of active regions observed by SOHO/MDI and SDO/HMI, and the corresponding soft x-ray solar flares observed by GOES generate the dataset for training. Previously CNN has been used for rough surface parameters and solar flare inversion. Here we demonstrate that it can also be used for Chla concentration inversion by learning the complex nonlinear relationship between the Rrs and the OC-CCI Chla values.

CNNs have a large number of hyperparameters and in this paper the patch size, kernel size, number of kernel, neuron number in the fully connected layer, number of layers, and parameters for training need to tune to obtain good performance for estimating global Chla concentration. Adopting grid search method to obtain the best combination of hyperparameters is very time-consuming. In this experiment, the hyperparameters are tuned in the flowing strategy. Structure of a CNN mentioned in Fig. 1 is the classical LeNet, and the default hyperparameters in this experiment are based on the parameters in Fig. 1. Number of layers will first to be changed manually from Fig. 1 to generate a few experiment groups with different depths. Then other hyperparameters in each group will be tuned to get the best performance CNN model which owns the structure described in Fig. 2.

In this paper, only January 2016 data are used for training and the whole-year 12-month data for testing. Still, CNN works very well, therefore CNN has strong robustness and generalization capability, indicating that it may be used for predicting long-term Chla concentration without the need to fine tune using the new observations. The mean $R^2$ of CNN is 0.901, which is very close to OC-CCI Chla concentration, suggesting that using MODIS data only, the CNN approach can approximate the OC-CCI products, which is probably due to the high learning capability of the CNN that can fully take advantage of the information in MODIS by exploiting the complex nonlinear relationship between the Rrs and the OC-CCI Chla values. In addition, the OC-CCI Chla concentration dataset is adopted as the ground truth because it is well-recognized Chla concentration product that is produced by assimilating different satellite data through a complex data processing steps. The GPU Nvidia Quadro M2000 is used in this study, with the GPU memory of 4 GB, NVIDIA CUDA cores of 768, and boost clock of 1180 MHz. A more powerful GPU like Nvidia GeForce GTX 1080 Ti with the GPU memory of 11 GB, NVIDIA CUDA cores of 3584, and boost clock of 1582 MHz will help deal with big data. With better GPU, more GPUs, and batch processing during the prediction stage, we can further shorten the processing time.

Although the mean RMSE and $R^2$ of CNN are acceptable, overestimation of CNN can be found from Figs. 4(a), 5, and 6(a) in the high Chla concentration areas, which may be caused by the imbalance quantities of the case 1 and the case 2 waters training samples. During sample selection, we scan the whole global Chla image to train the CNN model, resulting in the case 1 samples are much more than case 2 samples. More comprehensive studies, like combining sea surface temperature and digital bathymetry model data as input, are required to in the next stage to solve the overestimation problem. Although OC-CCI is the most accessible and reliable Chla concentration product, uncertainties may exist that will affect the performance of model. In the future, if more reliable data appear, further studies will be conducted to train and evaluate the CNN model.

6.

Conclusion

In this study, a CNN method was applied to MODIS Rrs images to estimate global Chla concentration. The CNN took the patches of four MODIS Rrs images as input and generated the Chla concentration directly. The OC-CCI Chla concentration image was used as ground-truth for training. A total of 12 monthly global Chla concentration images were produced and the generation of each image takes about 2 h. Qualitative and quantitative analyses were used for evaluation and comparison analysis, and the results implied that the CNN constitutes an accurate, fast, and robust method for the estimation of the global Chla concentration. Considering the big data characteristic of remotely sensed global Chla data, these characteristics of the CNN model is of great importance.

7.

Appendix

All experiments are conducted on a 64 bits Intel Xeon E5-2640 v3 workstation with 2.6 GHz of clock and 32 GB of RAM memory. A GPU Nvidia Quadro M2000 with 4 GB memory and version 8.0 CUDA. Caffe,⁵⁴ a popular deep learning framework especially for CNN which is short for convolutional architecture for fast feature embedding, under Ubuntu 14.04 LTS is used in this paper. Caffe provides a convenient way to implement the CNN models by defining the architecture of CNN such as the number of layers, the type of layer and the strategy for optimization. The preprocessing of the data and the numerically assessment are performed under MATLAB R2014a. The data type of CNN in Caffe for reading image patches and their corresponding Chla values is of the HDF5 format instead of the default Lightning Memory-Mapped Database format to achieve easily storage of multichannel images and float labels.