Neural Network Reflectance Prediction Model for Both Open Ocean and Coastal Waters

Remote sensing of global ocean color is a valuable tool for understanding the ecology and biogeochemistry of the worlds oceans, and provides critical input to our knowledge of the global carbon cycle and the impacts of climate change. Ocean polarized reflectance contains information about the constituents of the upper ocean euphotic zone, such as colored dissolved organic matter (CDOM), sediments, phytoplankton, and pollutants. In order to retrieve the information on these constituents, remote sensing algorithms typically rely on radiative transfer models to interpret water color or remote-sensing reflectance; however, this can be resource-prohibitive for operational use due to the extensive CPU time involved in radiative transfer solutions. In this work, we report a fast model based on machine learning techniques, called Neural Network Reflectance Prediction Model (NNRPM), which can be used to predict ocean bidirectional polarized reflectance given inherent optical properties of ocean waters. This supervised model is trained using a large volume of data derived from radiative transfer simulations for coupled atmosphere and ocean systems using the successive order of scattering technique (SOS-CAOS). The performance of the model is validated against another large independent test dataset generated from SOS-CAOS. The model is able to predict both polarized and unpolarized reflectances with an absolute error (AE) less than 0.004 for 99% of test cases. We have also shown that the degree of linear polarization (DoLP) for unpolarized incident light can be predicted with an AE less than 0.002 for 99% of test cases. In general, the simulation time of SOS-CAOS depends on optical depth, and required accuracy. When comparing the average speeds of the NNRPM against the SOS-CAOS model for the same parameters, we see that the NNRPM is able to predict the Ocean BRDF 6000 times faster than SOS-CAOS. Both ultraviolet and visible wavelengths are included in the model to help differentiate between dissolved organic material and chlorophyll in the study of the open ocean and the coastal zone. The incorporation of this model into the retrieval algorithm will make the retrieval process more efficient, and thus applicable for operational use with global satellite observations.


Introduction
The study of global ocean color measurements can help in understanding the ocean ecology and biogeochemistry which can subsequently help better understand the carbon cycle and its impact on climate change. Ocean polarized reflectances obtained from satellite sensors contain a fast ocean reflectance model that can be used in the joint retrieval algorithms to achieve operational retrieval of aerosol and ocean color information.
Recently, Fan et al. [28] have created a neural network to model the ocean body reflection matrix based on 3 parameters. In this work, we present an ocean reflectance model based on machine learning techniques, hereafter referred to as Neural Network Reflectance Prediction Model (NNRPM), which can be used to replace the radiative transfer simulation of ocean waters and expedite the retrieval algorithm significantly. This model incorporates additional parameters that would enable it to be better used for coastal waters. Moreover, it can also be used in semi-analytical models that retrieve IOPs from water leaving radiance or equivalent remote sensing reflectance. The basic advantages of neural networks are exploited, namely their accuracy performance properties, their robustness to noise, and their rapidity of execution.
This paper is organized as follows: Section 2 describes the methodology and definitions used in this work; Section 3 shows the result of this study; Section 4 summarizes the conclusion reached in this study.

Methods
Our goal is to build a polarized reflectance model for ocean waters. The ocean water model is the same as the one used in the MAPOL algorithm [27,29]. In order to model the optical properties of coastal waters the ocean water is assumed to consist of four components: pure sea water, phytoplankton particles, non-algae particles (NAP), and CDOM. The absorption coefficient a w (λ) and scattering coefficient b w (λ) for pure sea water are from the experimental data [30][31][32]. The backscattering fraction for water is 0.5 [33]. We denote the absorption coefficient for phytoplankton as a ph (λ), the total absorption coefficient for both CDOM and NAP as a dg (λ), the total backscattering coefficient and total backscattering fraction for both phytoplankton and NAP as b bp (λ) and B p (λ). The scattering coefficient of CDOM is assumed to be negligible. The bio-optical model can be written as [29]: where wavelength (λ) is in nm; a ph (λ) is determined by [Chla], A ph , and E ph ; [Chla] is in units of mg/m 3 . Values of the coefficients A ph , and E ph are obtained from Bricaud et al. [34]. S dg is the exponential spectral slope for a dg (λ); S bp and S Bp are the polynomial spectral slopes for b bp (λ) and B p (λ) respectively. The total absorption coefficient is a(λ) = a w (λ) + a ph (λ) + a dg (λ), and the total backscattering coefficient is b b (λ) = 0.5b w (λ) + b bp (λ). The inherent optical properties (IOPs), a ph (λ), a dg (λ), b bp (λ), a w (λ), and b w (λ) are used to describe the ocean color spectrum in ocean bio-optical models [11,15,35]. Generalised IOP (GIOP) model is also used in the earth science community [35,36].
The Fournier-Forand (FF) [37] phase function is used to represent the phase functions of both phytoplankton and NAP, which is in agreement with various in-situ volume scattering function (VSF) measurements [38]. The scattering coefficient of the phytoplankton and NAP particles can be derived from b p = b bp /B p .The determination of the FF phase function [39] is done by backscattering fraction B bp [39]. The total phase function for the coastal water system is obtained by mixing the FF phase function with the pure seawater phase function weighted by the scattering coefficients of the constituents. The normalized Mueller matrix is from measurements from average ocean waters [40][41][42]. The ocean water is assumed to be vertically homogeneous with a depth of 200 m. The minimum and maximum values of the parameters [29] considered are given in Table 1. Here, we are focusing on the coastal waters so that the particle size is primarily large. In addition to the aforementioned parameters, viewing zenith angles range from 0 to 87.5 • with an interval of 2.5 • . The viewing azimuth angles are from 0 to 180 • with an interval of 5 • . The incident zenith angles are taken from 0 to 85 • . The wavelengths used are 385, 410, 440, 470, 555, 670, 760, 863.5, and 870 nm. The aforementioned wavelengths include the ultraviolet, visible and near infrared (VNIR) part of the electromagnetic spectrum. The wavelengths 410, 470, 555, 670, and 863.5 nm are used as part of research scanning polarimeter (RSP) and the rest are adopted based on other existing and new ocean color sensors. The above mentioned bio-optical parameters, wavelengths, incident zenith angles (θ i ), viewing zenith and azimuth angles were provided as inputs to the radiative transfer code based on successive order of scattering for coupled atmosphere and ocean systems (SOS-CAOS) [43,44]. There is no atmosphere in our radiative transfer simulation as we are focused on the ocean reflectance model. Glint from the ocean surface can be predicted analytically, which is excluded in this simulation. The output of the radiative transfer simulations were considered as the true or desired value while training the neural network.
The objective of the neural network is to predict the non-dimensional effective Mueller matrix [45], M as a function of the wavelength λ, viewing zenith angle θ v , viewing azimuth angle φ v , incident zenith angle θ i , and incident azimuth angle φ i : ( We denote the radiance column vector as L = (I, Q, U, V) T , where (I, Q, U, V) are the Stokes parameters representing the polarization state of a light beam and the superscript T stands for transpose.
Once M is obtained, the upwelling radiance vector L(θ v , φ v ) can be calculated from incident downwelling radiance vector L(θ i , φ i ). For an arbitrarily polarized incident light, the reflected radiance can be calculated from the following relation: where we have implicitly assumed the wavelength dependency of both M and L. Hereafter, we neglect the fourth column and row of M, due to their small contribution for oceanic environment [46]. The 3 × 3 Muller matrix M can be obtained by considering the following incident polarizations: (I 1 , Q 1 , U 1 ) as complete unpolarization (1,0,0), (I 2 , Q 2 , U 2 ) as horizontal polarization (1,1,0) and (I 3 , Q 3 , U 3 ) as 45 • polarization (1,0,1) to get: where (I w1 , Q w1 , U w1 ), (I w2 , Q w2 , U w2 ), (I w3 , Q w3 , U w3 ) represents remote sensing reflectance just above the ocean water surface corresponding to the three polarized incident beams. F 0 is the incident irradiance. The elements of the matrix in Equation (4), are obtained by solving SOS-CAOS [43,44]. As part of the neural network training, all the elements (except M 11 ) of the matrix in Equation (4) are normalized by M 11 . We have created a supervised, model based, batch learning, neural network algorithm with the architecture of a multilayer perceptron (MLP) using MATLAB's Neural Network package [47]. An MLP is a feedforward and backpropagated artificial neural network [48]. This neural network approximates the function relating the various inputs of the network to the output [49] and would be used to predict ocean bidirectional reflectance distribution function in this work [50]. The network was trained with 10 inputs (Chla, a dg , b bp , B p , S Bp , S bp , S dg , θ i , θ v and φ v ) and the output reflectance for every wavelength being studied. Figure 1 shows the structure of the neural network used in this work. We arrived at the above architecture by using the general guidelines in [51] as the base and then modifying it as per the requirements of our dataset. The input layer has 10 inputs x 1 to x 10 . The hidden layer consists of 4 layers each having 40 nodes. The output layer has one node. Each node has a summation operator and an activation function within it except the input nodes. The output of every node in a layer is multiplied by a weight and summed with the output of other weighted nodes in that layer and a bias to be fed as input to the activation function. The output of the activation function becomes the input of a node in the following layer. we have used tangent sigmoid as our activation function: where x denotes the sum of weighted nodes of the layer and the bias.  For each training instance this neural network algorithm first makes a prediction on the forward pass, measures the error, then goes through each layer in reverse to measure the error contribution from each previous connection (on the reverse pass), and finally tweaks each of the connection weights to reduce the error as part of Levenberg-Marquardt optimization algorithm [52]. So, essentially, we start off with random weight values and then improve it gradually one step at a time such that each step reduces the error until the algorithm converges to a minimum. The Levenberg-Marquardt optimization algorithm can be represented as: where w k is the weight/bias of previous iteration and w k+1 is the weight/bias of next iteration. J is the Jacobian matrix and J T is its transpose. J contains first derivatives of the network errors with respect to the weights and biases, and e is a vector of network errors. I is the identity matrix.
The parameter µ determines direction as well as the length of the particular iteration step. It is a value that speeds up or slows down how quickly an algorithm learns. Technically, it determines the size of step an algorithm takes when moving towards a global minimum. Levenberg-Marquardt algorithm combines the Gauss-Newton and gradient descent optimizer methods. For larger values of µ, this algorithm will use gradient descent as the primary optimizer and for very small values of µ, Gauss-Newton would act as the primary optimizer. Here we have used µ = 10 −4 . After each iteration of training, the performance of neural network was monitored on the test dataset by calculating the mean square error. In order to avoid over-fitting, the training was stopped when the mean square error continued to increase for 6 consecutive iterations on the test dataset [51]. The model equation can be represented as: where y is the output reflectance, predicted by the network. N i is the number of input parameters, N h is number of neurons in each hidden layer, k] , and w 5 [1,k] are the weights of the four hidden layers and the output layer respectively, and b 1 4[k] , and b 5 [1] are the biases of the four hidden layers and the output layer respectively. f is the tangent sigmoid activation function given in Equation (5). Figure 2 shows the high level flow diagram of the neural network. Feature sequences for the various input parameters were generated using a quasi-random number sequence, the Halton sequence [53]. First introduced by Halton [54], it is constructed using a deterministic method and is of low discrepancy. Training and test datasets corresponding to the above sequences were generated by SOS-CAOS [43,44]. The neural network was trained and validated using the training and test datasets.
The error measurement for the element M 11 is computed using the relation below: In order to avoid the numerical difficulty of dividing by zero, the error measurement for all other normalized elements (M ij /M 11 ) is computed using the relation below: The above errors are used to calculate the empirical cumulative distribution function, which gives the probability of the number of samples within a certain (error) value, which is used to illustrate the performance of the neural network. If {X 1 , ..., X n } is an independently and identically distributed sample from an unknown distribution function, F(x) = P (X ≤ x), then the empirical cumulative distribution function (F n (x) ⇒ ECDF) is defined as [55]: where, Another metric used here to evaluate the performance of NNRPM is the Pearson correlation coefficient (R) [56]: where n is the sample size, x i and y i are sample points,x andȳ are sample means of two variables x and y respectively.

Results
In Figure 3, we present a plot between error percentage and Halton sequences. Here, error percentage represents the percentage of test cases with absolute error greater than 0.004. At 8 K Halton sequences, the error percentage reached its minimum value and then varies little after that. This established the following relationship with error percentage and training Halton combinations (N):   In Figure 9, we have shown ECDF versus PE at wavelengths 385 and 555 nm, which shows that ECDF approaches 1 asymptotically as PE increases and becomes very close to 1 when PE is larger than 5%. In Tables 2-4, the performance of the neural networks corresponding to each element of the Mueller matrix from Equation (4), for every wavelength being studied, has been shown. For evaluating the prediction performance of the degree of linear polarization (DoLP = Q 2 + U 2 /I), we tested the AE for the expression M 2 21 + M 2 31 /M 11 , which is equal to DoLP when the incident light is unpolarized. The results are shown in Table 2. For the first element M 11 , the NNRPM is able to predict values with percentage error (PE) less than 5% for 99% of the test cases. For the rest of the normalized elements, the NNRPM is able to predict values with absolute error (AE) less than 0.004 for 99% of the test cases. DoLP for unpolarized incident light can be predicted with an absolute error (AE) less than 0.002 for 99% of the test cases. The corresponding time taken by NNRPM is 0.0296 seconds which is mostly stable for an arbitrary case.

Discussions
In order to determine the optimum number of training cases required, we experimented with various sets of Halton sequences. We found that accuracy increases as the number of Halton combinations used in the training increases until it reaches 8000. From Figure 3, we can see that error initially decreases as the training data is increased and reaches a minimum for a certain volume of training data. This result is similar to the error's inverse relationship with number of training instances mentioned in Gallant et al. [57]. Figure 4 shows a linear relationship between the true value represented by SOS-CAOS and the predicted value represented by NNRPM. The value of R-squared being close to 1 indicates clear correlation between the true value (SOS-CAOS) and the predicted value (NNRPM). The slope (m) being close to 1 and bias (y) being close to 0 further confirms the agreement between NNRPM with SOS-CAOS.
The Bidirectional Reflectance Distribution Function (BRDF) is the quantity which geometrically characterizes the reflecting properties of a point on a surface [58]. The angular distribution of the water leaving signals are illustrated through the polar plots in Figures 5 to 8. The SOS-CAOS and NNRPM polar plots provide a pictorial representation of how close the neural network predicted value is to the actual SOS-CAOS value for the specific scenario under consideration. In Figures 5 and 6, the accuracy of the model for the element M 11 is illustrated by the Percentage error polar plots which has a maximum value around 3%. For test cases with the smaller values of a dg , the values of reflectance will be higher for both wavelengths, and vice versa. In Figures 7 and 8, the accuracy of the model for the element M 12 /M 11 is illustrated by the Absolute error polar plots which has a maximum value around 0.0009. We can see that the model is able to predict both positive and negative results equally well. Figure 9 is the visual representation of how quickly the cumulative distribution function increases to 1. It shows the probability of an unknown test dataset having a percentage error below any chosen threshold. Here we can say that the probability of getting a percentage error of less than 5% across the entire test dataset is close to 1. This further confirms the accuracy of the model on an unseen test dataset.
The simulation time of SOS-CAOS depends on the order of scattering and the number of Gaussian quadratures used for angular integration, which is in turn determined by the required accuracy. Moreover, it also depends upon the total optical depth of the system. Using the parameters mentioned in the results section, we can see that the incorporation of NNRPM would improve the speed of retrieval algorithms by a factor of 6000. For different ocean reflectance cases, the CPU time varies for the SOS-CAOS solution. However the order of magnitude is about the same. Thus the comparison provides a general guideline on the improvement in CPU time.
The future incorporation of this model into the retrieval algorithms will make the retrieval process more efficient, and thus applicable for operational use with global satellite observations. Future work could also examine the necessity of adding more IOPs to these models to simulate the behavior of different water bodies.

Conclusions
Ocean polarized reflectances obtained from satellite sensors contains information about the constituents of the upper ocean euphotic zone, such as colored dissolved organic matter (CDOM), sediments, phytoplankton types, and pollutants. Traditional ocean color remote sensing algorithms retrieve the constituent information by first performing atmospheric correction, then deriving IOPs from the resultant remote sensing reflectances. However, atmospheric correction algorithms are still a challenging task for scenes involving coastal waters and absorbing aerosols. In order to overcome this difficulty, joint retrieval algorithms of aerosols and ocean color have been designed to use information in polarized measurements. Joint retrieval algorithms need to solve the radiative transfer equation in atmospheric and ocean coupled system repeatedly, which is resource prohibitive, making them impractical for operational production of satellite data.
In order to resolve this problem, we have built a forward model using machine learning techniques, called Neural Network Reflectance Prediction Model (NNRPM), which can be used to predict ocean bidirectional reflectance, given inherent optical properties of ocean waters. The accuracy of the model is verified with the test cases generated from radiative transfer code based on the successive order of scattering technique (SOS-CAOS). The NNRPM model is able to predict the reflectance values for the first element with percentage error less than 5% for 99% of the test cases. For the rest of the normalized elements, the NNRPM is able to predict reflectance values with an absolute error less than 0.004 for 99% of the test cases. For unpolarized incident light, the model is able to predict DoLP with an AE less than 0.002 for 99% of test cases. When comparing the average speeds of the NNRPM against the SOS-CAOS model, we see that the NNRPM is able to predict the Ocean BRDF 6000 times faster than SOS-CAOS. The application of this neural network model in inverse algorithms would expedite the retrieval process and make them more efficient for use in global satellite observations. It can also be used in traditional semi-analytical models of retrieving ocean water IOPs, since our bio-optical model is more general than most existing models.