Deep Learning-Based Emulation of Radiative Transfer Models for Top-of-Atmosphere BRDF Modelling Using Sentinel-3 OLCI

Abstract

The Bidirectional Reflectance Distribution Function (BRDF) defines the anisotropy of surface reflectance and plays a fundamental role in many remote sensing applications. This study proposes a new machine learning-based model for characterizing the BRDF. The model integrates the capability of Radiative Transfer Models (RTMs) to generate simulated remote sensing data with the power of deep neural networks to emulate, learn and approximate the complex pattern of physical RTMs for BRDF modeling. To implement this idea, we used a one-dimensional convolutional neural network (1D-CNN) trained with a dataset simulated using two widely used RTMs: PROSAIL and 6S. The proposed 1D-CNN consists of convolutional, max poling, and dropout layers that collaborate to establish a more efficient relationship between the input and output variables from the coupled PROSAIL and 6S yielding a robust, fast, and accurate BRDF model. We evaluated the proposed approach performance using a collection of an independent testing dataset. The results indicated that the proposed framework for BRDF modeling performed well at four simulated Sentinel-3 OLCI bands, including Oa04 (blue), Oa06 (green), Oa08 (red), and Oa17 (NIR), with a mean correlation coefficient of around 0.97, and RMSE around 0.003 and an average relative percentage error of under 4%. Furthermore, to assess the performance of the developed network in the real domain, a collection of multi-temporals OLCI real data was used. The results indicated that the proposed framework has a good performance in the real domain with a coefficient correlation ($R^2$), 0.88, 0.76, 0.7527, and 0.7560 respectively for the blue, green, red, and NIR bands.

1. Introduction

The Bidirectional Reflectance Distribution Function (BRDF) defines the surface reflectance anisotropy as a function of sun-surface-sensor geometry and surface and sensor properties [1]. The BRDF plays an important role to geometrically normalize satellite observations to a standard direction. This is useful to various remote sensing applications [2], to estimate surface energy parameters such as albedo [3], to estimate accurate biophysical and structural surface parameters [4], and to improve land cover classification and segmentation [5]. However, in more remote sensing applications, the surface reflectance behavior is assumed to be Lambertian (isotropic, same reflectance in all directions), resulting in uncertainties in the analysis [6]. Theoretically, a set of high-quality, sufficient, and well-distributed remote sensing observations over a short period of time should exist to describe the BRDF accurately, along with a reliable model [7].

The literature reports several simple widely used empirical and semi-empirical models at national, regional, and global scales for BRDF modeling [8]. For example, Ross-Thick Li-Sparse (RTLS) [9] and Rahman-Pinty-Verstraet (RPV) [10], and Walthall [11] are used for BRDF modeling due to their simplicity and low computational cost. However, semi-empirical models are simplified versions of physically based radiative transfer models (RTMs) and are derived by imposing some prior assumptions on some RTMs parameters and are mainly developed for bottom of atmosphere (BOA) BRDF. However, BOA reflectance requires atmospheric, topographic, and adjacency effect correction that may lead to introducing sources of uncertainty and errors in BRDF modeling [12]. Additionally, the early developed BRDF models were based on observations collected by single satellite sensors such as MODIS [13], MIRS [14], or VIIRS [15] in a long-time window to providing more observations with minimum cloud and cloud shadows. However, long time windows can introduce certain uncertainties in BRDF estimation due to surface and atmosphere variability. Recent growth and progress in geospatial Cloud Computing Platforms (CCPs) such as Google Earth Engine (GEE) [16], Amazon Web Services (AWS) [17], and Microsoft Planetary Computer (MPC) [18] have facilitated access to multi-sensor, multi-angular, and multi-spatial resolutions satellite observations that can lead increasing the number of observation and shorten the data collection time window and then minimizing the variability effect of surface and atmosphere. The CCPs like GEE also can provide a wealth of information about atmospheric and surface parameters such as LAI or aerosol optical thickness (AOT) that can be useful in the TOA BRDF modeling [19] Furthermore, to overcome the issues associated with the BOA reflectance BRDF modeling, the directional TOA radiance/reflectance directly can be used to model the BRDF. In this way, the RTMs that express the propagation and interaction of solar radiation with different mediums such as surface or atmosphere can be used to simulate BOA or TOA reflectance [20,21]. For example, PROSAIL is widely used in remote sensing to link directional ground surface reflectance, i.e., BOA, with structural and biophysical properties of the surface [22,23,24]. Moreover, an atmospheric model such as the 6S (Second Simulation of a Satellite Signal in the Solar Spectrum, version 6SV from 2005) can describe solar radiation propagation through the atmosphere [25]. The coupling of PROSAIL and 6S enables accurate BRDF modeling by simultaneously linking up the satellite sensor TOA directional reflectance with atmospheric and surface parameters and BOA directional reflectance [26]. However, the complexity, high dimensionality, and computationally cost of physical RTMs, particularly in large areas with big datasets, restrict their use in real remote sensing applications, especially for BRDF modeling. To overcome these issues, it has been suggested to approximate the physical RTMs by means of machine learning known as emulation [27]. Machine learning (ML) algorithms such as neural networks have demonstrated a good ability to emulate the nonlinear and complex physically based RTMs [28].

The emulation-based approach can train a ML model for different goals, including, approximating the complex physical models [29], retrieving or quantifying the vegetation biophysical variables [30,31]. For instance, in [32,33,34] Support Vector Machines (SVMs), Random Forests (RFs), and Artificial Neural Networks (ANNs) were trained using simulated data generated by physical RTMs to retrieve surface or atmospheric parameters. In [35], an ANN emulation methodology was developed to simulate the climate forecast system (CFS). The results found that the ANN can be used as an emulator of climate simulations and seasonal predictions with an accurate performance. Furthermore, the emulation approach was used in different study to generate the synthetic spectral dataset [36] or synthetic scene generation [37].

However, traditional ML algorithms based on SVM, RFs, and shallow ANN have a straightforward structure that is incapable of discovering and mapping complex and nonlinear relationships between input and output variables from high-dimensional RTMs. A few recent studies have investigated the potential of new and advanced ML algorithms such as deep neural networks (DNNs) as emulators of physically based model. For example, in [38], a fully connected deep neural network was proposed to emulate the Community Radiative Transfer Model (FCDN_CRTM) using simulations of brightness temperatures (BTs) over ocean surfaces. In another research [39], several deep learning-based networks were trained to emulate the lookup table generated by the RTM of the Multi-Angle Implementation of Atmospheric Correction (MAIAC).

Similar to the aforementioned studies, a new deep learning-based approach can be developed to emulate the physically based BRDF model. Many DNNs, including convolutional neural networks (CNNs), recurrent neural networks (RNNs), and long short-term memory (LSTM), have been developed in various application domains, such as computer vision [40], natural language [41] and remote sensing applications [42]. Among the DNNs, CNNs have received more attention and success including in remote sensing applications. Two-dimensional (2D) CNNs are helpful for image-based processing tasks, such as image classification [43], object detection [44] and image semantic segmentation [45], multivariable regression [46]. However, one-dimensional (1D) CNNs are commonly used for 1D signal data such as audio, text, and 1D time-series data, and similarly can also be used for regression data analysis [47]. Since the input and output variables of physical BRDF—obtained by coupling PROSAIL+6S—can be treated as 1D-regression problems, a 1D-CNN can be trained using a representative BRDF dataset generated by PROSAIL+6S. Once the network training, testing, and validation stages are accomplished based on a large dataset in the simulation domain, the network can be transferred to the real domain using the transfer learning (TL) [48]. The TL technique is helpful when a limited amount of data for one task is available in one domain (i.e., real domain), but a large volume of data can be accessible for a similar task in another domain (i.e., simulation domain) [49]. This paper mainly focuses on developing an efficient DNN framework for BRDF characterization in the simulation domain. However, to assess the performance of the proposed network in the real domain, the network was transferred to the real domain and applied to a series of observations collected during two weeks from Ocean and Land Colour Instrument (OLCI) sensor carried by Sentinel-3A (S3A) and Sentinel-3B (S3B) satellites. Thus, this study first begins by simulating a wide range of directional spectral reflectance for a variety of surface, atmospheric and geometric (illumination and viewing angles) configurations using the coupled PROSAIL+6S model. Second, the generated dataset is used to design, train and validate a 1D-CNN architecture for BRDF modeling. The third step consists of evaluating the trained 1D-CNN using an independent dataset also created by the coupled PROSAIL+6S based on unknown configurations to the network in the simulation domain. In the last step, the proposed network is applied to real data consisting of OLCI observations.

The remaining sections go as follows: Section 2 regards the research’s theoretical background; Section 3 describes the methodology; Section 4 presents the results and analysis in the simulation and real domain; with the conclusions in Section 5.

2. Theoretical Background

2.1. PROSAIL

PROSAIL is one of the widely used and standard RTMs in remote sensing applications; it combines radiative transfer models of the leaf (PROSPECT) and canopy (SAIL) [24]. PROSAIL defines directional reflectance as a function of leaf property, canopy parameters, sun-surface-sensor observation geometry, and soil factors (${\rho }_{soil}$) [26]. The leaf parameters include the internal structure parameter of the leaf mesophyll (N), chlorophyll a + b concentration (${\mathrm{C}}_{ab}$), leaf equivalent water thickness (Cw), dry matter content ($C_m$), and brown pigments ($C_{bp}$). The canopy parameters comprise leaf area index (LAI), leaf inclination distribution function (LIDF), and average leaf inclination angle (ALIA). The geometry parameters also contain view zenith angle (VZA), sun zenith angle (SZA), and relative azimuth angle between sun and sensor (RAA) [50] Using PROSAIL in forward mode can simulate directional reflectance and provide training and testing datasets for an ML algorithm [51]. This method is particularly useful when satellite sensors cannot provide sufficient observations for BRDF modeling due to limited viewing geometry, cloud contamination, or poor atmospheric conditions [52]. PROSAIL can simulate only BOA directional reflectance; thus, requiring an additional step to couple the PROSAIL with an atmosphere RTM like the 6S code to produce the top of atmosphere (TOA) directional reflectance and to represent satellite data properly.

2.2. Atmospheric RTM: The 6S Code

The Second Simulation of a Satellite Signal in the Solar Spectrum (6S) was introduced by Vermote [53] to describe the scattering and absorption effects of atmospheric gases and aerosols on solar radiation propagation downward and upward through the atmosphere under cloudless conditions [54]. The multiple interactions of solar radiation with atmospheric layers are described by a method called the Successive Orders of Scattering (SOS). In 2005, a vector version of 6S (6SV) that accounts for radiation polarization was proposed [55]. The inverse mode of 6S can be used in retrieving atmospheric parameters, such as aerosol optical thickness (AOT), columnar water vapor (CWV), and O3 column concentration or atmospheric correction of satellite images [54]. However, through the forward application of 6S, it is possible to simulate TOA directional reflectance for various input conditions and use them as a training and testing database [56]. To account for the anisotropic effect of the surface in 6S calculation, the PROSAIL can be coupled with 6S to accurately describe the TOA directional reflectance. Figure 1 depicts a schematic overview of coupling PROSAIL and 6S and their parameters. According to Figure 1, after calculating the surface directional reflectance by PROSAIL (combination of PROSPECT and SAIL), the outputs are used as ground reflectance, which is one of the pre-requirements for running 6S.

2.3. One-Dimensional Convolutional Neural Networks (1D-CNNs)

In recent years, one-dimensional convolutional neural networks (1D-CNNs) have received more attention in remote sensing due to their low computational cost, small training dataset requirement, and compatibility in dealing with applications mainly based on one-dimensional spectral or time-series data from one-pixel imagery or field measurement data [57]. In contrast to two-dimensional (2D) CNNs, which are appropriate for image-based applications and employ convolutional filters over the two dimensions, the 1D-CNNs extract the spectral or temporal features by applying filters along one dimension of the data. However, as 2D-CNNs, the 1D-CNNs model includes the five following layers: (1) an input layer that receives the input data in a tensor format and then feeds it to the network; (2) hidden layers consisting of several learnable convolution filters to convolve with input data and extract hidden information and high-level features; (3) pooling layer (if necessary) to downscale the dimensionality of features derived by the convolutional filters while highlighting the most dominant features and reducing the noise. The 1D-CNN also includes (4) a flattening and fully connected layer that transfers the results of the convolutional layer into a one-dimensional vector form to (5) the output layer, which contains the same number of target variables [58]. Such architecture is used in this work, as shown in the methodology. Additionally, most CNN-based models use a particular activation function called Rectified Linear Units (ReLU), which outperforms previous activation functions such as Sigmoid and Tanh. Regularization techniques, like the dropout technique, that accidentally assign zero weights to many neurons and drop them during the training process can be used in 1D-CNN to significantly mitigate the overfitting problem and improve the generalization accuracy of the network [59].

3. Material and Methodology

The proposed methodology in this article is divided into three phases. The first phase involves simulating a sufficient dataset of directional reflectance values for a variety of sun-sensor-surface and atmospheric configurations using coupled leaf canopy (PROSAIL) and atmospheric RTMs (6S). The second phase consists of designing and training-validating a 1D-CNN architecture for BRDF estimation using the simulated dataset from the previous phase. The trained and validated network was evaluated using the k-fold cross-validation method and independently generated testing datasets unknown to the network. In the third phase, the final designed network was transferred from the simulation domain to the real domain to evaluate the performance of the network for real satellite data. In this way, the general steps of real OLCI data collection, sampling and network tuning and retraining were taken to evaluate the network performance in the real domain. Figure 2 depicts the flowchart of the proposed methodology in both simulation and real domain.

3.1. PROSAIL and 6S Parameterization to Generate Synthetic OLCI Data

Before designing the 1D-CNN in the second phase of the study, the training and validation datasets in the simulation domain were generated by running the coupled PROSAIL and 6S models in the forward mode using [60]. In this study, we considered only PROSAIL and 6S input parameters with the higher importance according to global or local sensitivity analysis found in the literature [61]. To ensure generalizable results, the incorporated input variables were drawn randomly from a uniform distribution extending to the minimum and maximum range of each parameter, as listed in

Table 1. Description and ranges of the PROSAIL+6S model key input parameters used to simulate BRDF.

Parameters	Symbol	Description	Unit	Range	Distribution
Sun-sensor geometry	SZA	Solar zenith angle	Deg	10–50	Uniform
	VZA	Viewing zenith angle	Deg	0–60	Uniform
	RAA	Relative azimuth angle	Deg	−185–185	Uniform
Leaf properties	$C_{ab}$	Chlorophyll a, b content	μg ${\mathrm{cm}}^{-2}$	3–30	Uniform
	N	Leaf structure parameters	unitless	1–3	Uniform
	$C_{ar}$	Leaf Carotenoids content	μg ${\mathrm{cm}}^{-2}$	10	Fixed value
	$C_w$	Equivalent water thickness	cm	0.002	Fixed value
	$C_m$	Dry matter content	g ${\mathrm{cm}}^{-2}$	0.002	Fixed value
	$C_{bp}$	Brown pigments	unitless	0	Fixed value
Canopy architecture	LAI	Leaf area index	${\mathrm{m}}^2$ ${\mathrm{m}}^{-2}$	0.1–6	Uniform
	$LID{F}_a$	Average leaf slope	Deg	15–60	Uniform
	$LID{F}_b$	Leaf inclination distribution	Deg	0	Fixed value
	Hspot	Hot spot parameter	unitless	0.01–0.9	Uniform
	${\rho }_{soil}$	Dry/Wet soil factor	unitless	0.5	Fixed value
Atmospheric	AOT	Aerosol Optical Thickness	unitless	0.1–0.5	Uniform

[44]. From the geometrical parameters, the VZA ranged between 0–60° to cover the viewing geometry of wide-angle remote sensing sensors likes Sentinel 3-OLCI, which has a maximum VZA of 55° [62]. The RAA range was set from −185° to +185°, allowing for the demonstration of TOA reflectance anisotropy at all possible azimuth angles and even slightly more than possible values (−180° to +180°) to avoid discontinuity on the borders. The SZA value also was taken between 10° and 50° which cover the seasonally varying SZA in the study area.

To determine the range of canopy structure parameters, LAI varied between 0.1–6, obtained from the literature [63,64,65]. The range of hot-spot, related to the ratio of mean leaf size over canopy height, varies between 0.01 (small leaves, tall canopies) and 0.9 (large leaves, short canopies) [66]. The leaf inclination distribution function (LIDF) consists of two input components, $LID{F}_a$ (average leaf angle) varied from 15° to 60° and $LID{F}_b$ (distribution bimodality) parameter was ignored [51].

From the leaf parameters, only $C_{ab}$ and leaf structure parameter (N) was considered as float variables for simulation. The $C_{ab}$ was ranged from 3 to 30. The N, which is linked to leaf characteristics like as thickness and intercellular properties, varied between 1 to 3, following the literature [67,68]. Due to the close relationship between $C_{ab}$ and $C_{ar}$ ($R^2$= 0.86), reported in [69], the $C_{ar}$ was assumed a constant value of 10 μg ${\text{cm}}^{-2}$. The remaining three leaf parameters, due to their negligible influence on directional reflectance in the visible and near infrared were fixed ($C_w=0.002$, $C_m=0.002$, and $C_{bp}=0$) according to values in different studies [70]. The background soil reflectance was assumed to be Lambertian, and a typical soil related to dry/wet factor value of (${\rho }_{soil}$ = 0.5) was set [71].

To set atmospheric conditions, only the AOT, which is the most influential optical atmospheric parameter in visible and near-infrared (NIR) spectral bands, was ranged between 0.1 to 0.5 [43]. A continental aerosol mixture and a standard predefined atmospheric model (US Standard1962) were selected. Water vapor was kept to the default value in the 6SV code (1.42 g/cm²).

Therefore, nine critical input variables consisting of three solar-sensor geometry angles (SZA, VZA, RAA), three canopy properties (LAI, LIDF, Hspot), two leaf traits (N, C_ab), and one atmospheric parameter (AOT) were selected for simulations. According to the finding of [51,55,57], a total number of 100,000 samples can help to reach a sufficient accuracy for canopy variables estimation. In addition, it covers a diverse and representative leaf-canopy-atmosphere-geometry configurations for such a training dataset for BRDF modeling. Thus, in this study, a combination of 100,000 input parameters was used to run the coupled PROSAIL+6S in forward mode to simulate the corresponding BRDF. Note that the generation of these simulations took about 90 h.

Since PROSAIL and 6S simulate directional reflectance between 400 nm and 2500 nm with 1 nm interval, the output should be resampled to a multispectral version dataset using the desired sensor Spectral Response Function (SRF) to match real-sensor observations [72]. This study uses the SRF of the OLCI as shown in Figure 3 for the resampling process.

The simulated spectral directional reflectance was resampled to the four OLCI bands, including Oa04, Oa06, Oa08, and Oa17 using the following equation:

$${\rho }_S\left(i\right)=\frac{{{\displaystyle \sum }}_{Ai}^{Bi}{w}_i\left(\lambda \right)\cdot E_0\left(\lambda \right)\cdot \rho \left(\lambda \right)}{{{\displaystyle \sum }}_{Ai}^{Bi}{w}_i\left(\lambda \right)\cdot E_0\left(\lambda \right)},$$

(1)

where ${\rho }_S\left(i\right)$ denotes the resampled directional reflectance value with respect to the specific spectral band i, $\rho (\lambda$) is the TOA directional reflectance, $E_0\left(\lambda \right)$ is the TOA spectral solar radiance, and $w_i\left(\lambda \right)$ indicates the weight of the band SRF of OLCI where A and B are respectively the minimum and maximum wavelengths of OLCI in each i band.

3.2. Real Sentinel 3-OLCI Data for Validation and Application

While OLCI synthetic data was used to develop the proposed 1D-DNN network in the simulation domain, a collection of real observations acquired from the OLCI sensor was used to ensure that the proposed network also performs in the real domain. The OLCI has a spatial resolution of about 300 m at the nadir and a swath width of 1270 km. It covers 21 spectral bands with wavelengths ranging from 390 nm to 1040 nm [73]. In this study, as given in

Table 2. The selected spectral bands of OLCI for this study.

Spectral Bands	Spectral Range (nm)	Center (nm)	Width (nm)	Spatial Resolution (m)
Oa04	438–448	490	10	300
Oa06	555–565	560	10	300
Oa08	660–670	665	10	300
Oa017	856–876	865	20	300

, only four OLCI bands, including Oa04, Oa06, Oa08, and Oa17 were used to represent respectively the blue, green, red, and NIR bands to cover the visible and NIR spectrum.

The minimum temporal resolution of OLCI at the equator over land and ocean is about 1 and 2 days, respectively, while at high altitude parts of the Earth, due to the rising overlap in adjacent swaths, it can reach less than 0.6 days over land and 1 day over the ocean with two-coplanar S3 satellites [73]. Therefore, a fixed ground point on the surface can be viewed at many different angles during the orbital repeat cycle especially at high latitudes, where there is a significant amount of overlap between images acquired from adjacent orbital paths. In addition, the OLCI has a push-broom imaging system with a field of view of 68.5° and the maximum VZA is limited to a 55°, allowing the acquisition of off-nadir observations, which is essential for BRDF modeling [62].

3.2.1. Collection of OLCI Data for Application in the Real Domain

The proposed method was applied to the OLCI real data to evaluate the performance of the developed 1D-CNN in the real domain. To achieve this, at first, a two-week collection of OLCI observations were collected from an area located in Canada with high latitude. This helps to ensure that a significant amount of overlap between images acquired from adjacent orbital paths is found and more observations from different angles can be collected for each pixel. Figure 4 shows the collected OLCI data and the study area, located in the province of Quebec, between the 46.04° and 47.8°N and 73.8° and 77.16°W that covers the part of the Boreal shield Ecozone.

All acquired OLCI images by both S3A and S3B satellites from 1st to 16th July 2019 over the study area were downloaded using Copernicus Open Access Application Programing Interface Hub (API) with a python library, namely, Sentinelsat (https://pypi.org/project/sentinelsat/) (accessed on 9 December 2022), that enables to query and download the Sentinel data.

3.2.2. Sampling of OLCI Data

The collected time series of OLC data in the real domain include different land covers with a large number of pixels recorded from different angles. For the aim of efficient processing and following the simulation domain experiments to achieve comprehensive coverage of the viewing and illumination angles, a random sampling method was developed based on the histogram of variables to choose an optimal subset of observations in the vegetation-covered areas, particularly forests and wetlands. This sampled dataset was used for post-fine-tuning and re-training the pre-designed network in the real domain. Therefore, under the assumption that the surface and atmosphere do not vary over two weeks, a subset of n pixels was selected from the image domain in both spatial and time dimensions according to the following sampling rules inspired by [74]:

• To determine the optimal n pixels, the average RMSE between the measured and the estimated reflectance by 1D-CNN in four bands was considered.
• The normalized Difference Vegetation Index (NDVI) of the selected pixels should be greater than 0.5 to cover the dense vegetation areas.
• Geometrical parameters, including VZA, SZA, and RAA values were divided into 10 bins using the histogram analysis of each variable, and then for each bin, a defined number of pixels were chosen randomly.

3.3. 1D-CNN Developing

Building a robust, accurate, and well-generalized CNN model is dependent on several factors, including preprocessing the training dataset, building the network architecture, and finally, tuning hyper-parameters such as learning rate, batch size, etc. [75]. The following sections describe these critical factors.

3.3.1. Data Preprocessing for the 1D-CNN

Differences in units and ranges of the input variables of PROSAIL+6S can prevent the network from properly learning from all features as required, resulting in poor performances and high generalization errors. First, scaling all dataset variables to specific ranges and orders of magnitude can mitigate these issues, lead to good performance, and accelerate network training and convergence. The ranges of all input variables (of uniform distributions) were rescaled between [0, 1] using Min-Max normalization [76], Second, to avoid the network from being influenced by random outliers generated by the coupled PROSAIL+6S, each output variable was analyzed independently to identify and remove the outliers before running the network with training and validation data. Third, instead of feeding entire datasets to the network in one iteration (one batch), they were divided into small batches and then introduced to the network in several iterations to increase the network’s speed and ensure efficient gradient descent optimization. Additionally, before training the network, the batched dataset was shuffled to enhance the training convergence and improve the network’s generalization ability. Finally, the shape of the simulated data, which was a matrix of n (number of samples) by m (number of variables), had to be matched to the 1D-CNN, which accepts the input data as a three-dimensional tensor (batch size, channels, length).

3.3.2. 1D-CNN Architecture Design

Finding the optimal CNN architecture highly depends on the problem to solve. This study used a trial-and-error approach inspired by previously developed 1D-CNNs [47] to establish an efficient and well-generalized architecture for BRDF characterization. As a result, various 1D-CNN architectures were designed and evaluated with a variety of depths, widths, kernel sizes, and a number of hidden, convolutional and pooling layers. Most conducted studies to develop a CNN-based model have proposed expanding the network feature space (channels) while decreasing the height and width of initial inputs by going deeper [77]. In light of this, we started with a network with a one-layer convolutional layer neural network as basic architecture, and then different convolutional layers were added to find a suitable architecture. Therefore, designing the architecture began with a shallow network and gradually expanded its depth to find more complex patterns and features. During the architecture design, it is recognized that increasing the network depth leads to more feature extraction, but the number of model parameters (weights and biases) also increases, resulting in more challenging in training and generalizing the network. Additionally, the trials indicated that increasing the number of filters improved the accuracy of the 1D-CNN model but diminished the model generalization performance. The architecture of final designed network started by reshaping the input data into a tensor and feeding it to a convolutional layer (Figure 5). The following two convolutional layers contained 32 and 64 filters to derive additional features from the input data. Next, two layers were also 1D-convolutional layers, involving 128 and 256 filters with the same kernel size, stride, and padding values of respectively 2, 1, and 1 (See

Table 3. Details of the proposed (1D-CNN) architecture.

Layer Type	Layer Configuration				Output Size		Learnable Parameters
	Filters size	Kernel	Stride	Pad	Channel	Length
Reshaped input layer	-	-	-	-	3	3	0
Conv-1D	32	2	1	1	32	4	(2 × 1 × 3 + 1) × 32 = 224
Conv-1D	64	2	1	1	64	3	(2 × 1 × 32 + 1) × 64 = 4160
Conv-1D	128	2	1	1	128	4	(2 × 1 × 64 + 1) × 128 = 16,512
Conv-1D	256	2	1	1	256	5	(2 × 1 × 128 + 1) × 256 = 65,792
Max-pooling		2	2	0	256	2	0
Flatten	-				512	1	0
Dense	512 neurons				-		512 × 128 + 128 = 65,664
Dense	128 neurons				-
Output layer	4 neurons				-		0
Total parameters	-				-		152,352

for more details). A polling layer was added to the network to down sample the data while preserving the most prominent features. Moreover, one flatten layer was also used to reshape the multidimensional output of previous layers into a one-dimensional array and prepare them for the next two fully connected dense layers, with respectively 512 and 128 neurons. Finally, the output layer was composed of four neurons corresponding to the four bands of the OLCI. The non-linear ReLU (Rectified Learned Unit) function: f(x) = max (0, x), (where x is the input to a neuron) was used as an activation function in all hidden layers due to its outstanding performances. Figure 5 illustrates the proposed 1D-CNN architecture. Table 3 provides additional information regarding the layer type and configuration, such as filters, kernels, stride, and pad size used in the proposed 1D-CNN architecture. The last column of Table 3 lists the number of learnable parameters for each layer.

3.3.3. 1D-CNN Hyper-Parameter Tuning

The network learning process is controlled by manually tuning a set of hyper-parameters, such as batch size, epoch number, type of loss and optimization function, and learning rate that defines step size by adjusting network parameters to minimize the loss function [78].

The network performance was monitored in terms of training time as well as training and testing accuracy, using different values for each of the hyper-parameters through different scenarios.

Table 4. Summary of the hyper-parameters of the proposed 1D.

Hyper-Parameters	Tested Values	Selected Values
Learning rate	0.01, 0.001, 0.0001	0.001
Epochs number	100, 250, 500	100
Batch size	16, 32, 64, 128, 256	32
Optimizer	SGD, ADAM, AdaMax	AdaMax
Loss function	MAE(L1Loss), MSE	MSE

provides a summary of tested and selected hyperparameters. For example, in one scenario, the learning rate was dropped by magnitudes of 0.01, 0.001 and 0.0001, and the batch size values were 16, 32, 64, 128 and 256. This hyperparameter diversity helped to monitor the network performance in different conditions. For instance, if the network was training too slowly or converged too fast, the batch size or the learning rate was changed.

All scenarios were implemented on a Windows PC with Intel Core i7-3820 CPU @ 3.60 GHz and 32 GB RAM. The PC GPU was NVIDIA GeForce GTX 1080.

In addition, the optimization model and loss function had a direct effect on how the network weights and biases were adjusted. Several different optimizers, including Stochastic Gradient Descent (SGD) and Adam (Adaptive moment estimation), demonstrated a good training performance on a few occasions, but the network encountered a generalization problem. Ultimately, an improved optimization method, known as AdaMax (Adaptive moment estimation with Maximum), demonstrated satisfactory performance in terms of accuracy, mitigating overfitting and generalization problems when compared to the other optimization models tested [79].

Furthermore, the mean squared error (MSE) function was used to calculate the network’s loss over 100 epochs. For the real domain, the same hyper-parameters were used, except for the batch size that was changed from 32 to 16 to fit with the real data size.

3.4. Performance Evaluation

The performance of the proposed 1D-CNN was examined independently using a set of testing samples. In this way, three criteria, including root-mean-squared error (RMSE), and correlation coefficient (R²), according to the following equations, were used to evaluate the degree of correlation and error between a certain number (N) of estimated directional reflectance $\left({\widehat{y}}_i\right)$ and simulated directional reflectance $(y_i$):

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left({\widehat{y}}_i-y_i\right)}^2}{N}},$$

(2)

$$R=\frac{\left(N{{\displaystyle \sum }}_{i=1}^N\left({\hat{y}}_i{y}_i\right)\right)-\left({{\displaystyle \sum }}_{i=1}^N{\hat{y}}_i\right)\left({{\displaystyle \sum }}_{i=1}^N{y}_i\right)}{\sqrt{\left(N{{\displaystyle \sum }}_{i=1}^N{\hat{y}}_i{ }^2-{\left({{\displaystyle \sum }}_{i=1}^N{\hat{y}}_i\right)}^2\right)\left(N{{\displaystyle \sum }}_{i=1}^N{y}_i^2-{\left({{\displaystyle \sum }}_{i=1}^N{y}_i\right)}^2\right)}},$$

(3)

4. Results

4.1. Network Training and Testing in the Simulation Domain

The entire simulated dataset was divided into two subsets: training and validation, in a proportion of respectively 80% and 20% of the total simulated dataset. To ensure an unbiased dataset splitting, the dataset was randomly shuffled, then divided into training and validation data. The testing dataset was generated using an independent procedure next presented in the network testing Section 4.2. As it was discussed in the previous section, different scenarios, including varying network architecture and tuning network hyper-parameters with various values, were evaluated to build an optimal network for BRDF modeling. As part of these experiments,

Table 5. Network performance for different learning rates and batch sizes (the best learning rate and batch size are highlighted in green color).

Learning Rate
Batch Size		0.01			0.001			0.0001
		Training Accuracy	Testing Accuracy	Time (Minutes)	Training Accuracy	Testing Accuracy	Time (Minutes)	Training Accuracy	Testing Accuracy	Time (Minutes)
	16	0.997	0.560	83.62	0.997	0.942	83.05	0.998	0.925	80.53
	32	0.998	0.670	38.27	0.998	0.987	37.92	0.988	0.811	38.08
	64	0.998	0.790	20.01	0.998	0.876	20.22	0.998	0.853	20.28
	128	0.998	0.910	12.37	0.999	0.925	11.47	0.997	0.864	11.51
	256	0.999	0.871	7.63	0.998	0.718	7.65	0.997	0.861	7.17

gives information about the outperformed network performance in terms of training, testing average accuracy, and training time for values of the two most important hyper-parameters: batch size and learning rate. As can be seen from Table 5, the training accuracy is similar during most of the experiments with a value of 0.99 or higher. Testing accuracy varies from very low (0.560) to very high (0.987), with a more stable performance for the learning rate of 0.001, except for the large batch size of 256 which led to decreasing testing accuracy. The latter also degraded when the batch size was increased above 256 (results not shown). The training time ranges between approximately 7 min, for the large batch size of 256, and more than 80 min, for the small batch size of 16. In addition, a high value of learning rate with low a value of batch size required more time to converge and created a noisy learning curve and an unstable training procedure (curves not shown). Finally, the highest testing accuracy was obtained for a batch size of 32 and a learning rate of 0.001 (highlighted with green color) during 100 epochs and 37.92 min. This model has a more generalization ability.

Figure 6 illustrates the variation of the loss function (MSE) and accuracy (${\mathrm{R}}^2 score$) during training and validating the proposed 1D-CNN model with a learning rate of α = 0.001 and a batch size of 32. The shape and evolution of these curves provide information about the model performance and, in general, about the network under- or over-fitting issues. As can be seen, both values of the loss function using training and validation datasets decrease sharply, and the accuracy curve comes near 1, indicating that the model rapidly learned and converged. The training and validation loss curves decrease rapidly until about 40 epochs, then flatten out and gradually decrease. We noticed that the network did not exhibit an under-fitting problem since the training loss is not significantly higher than the validation loss. Furthermore, the validation loss was not considerably higher than the training loss, ensuring the proposed 1D-CNN model did not suffer from overfitting; this is a sign of proper hyper-parameter selection and network learning.

Figure 7 displays the scatter plots of the simulated BRDF values by the 1D-CNN and the simulated BRDF values by PROSAIL+6S RTMs for four bands of the OLCI, including Oa04, Oa06, Oa08, and Oa017. The least square regression lines (Red) show a high agreement between the simulated and the estimated BRDF values along all bands as the coefficients of determination $\left({\mathrm{R}}^2\right)$ are over 0.99 and the RMSE values are about 0.003.

4.2. Simulation Testing Dataset for Network Evaluation

The performance and generalization ability of the proposed 1D-CNN model were evaluated using independently generated test data hidden from the network during the training and validation processes. To create the testing dataset, the values in

Table 6. Input variables values for the testing dataset.

N	$\mathbf{Cab} (\mathsf{\mu }\mathbf{g} {\text{cm}}^{-2})$	LIDFA (°)	Hspot	$\mathbf{LAI} ({\mathbf{m}}^2$${\mathbf{m}}^{-2})$	AOT	SZA (°)	VZA (°)	RAA(°)
1.5	10.5	30.5	0.3	3.5	0.15	30.5	0–50°	−180–180°

were selected for input variables, after ensuring that these values were not included in the training nor validation dataset.

The generated directional reflectance test dataset is visualized in polar coordinates to illustrate the visual anisotropy of the surface reflectance and to ease the interpretation, analysis, and comparison of the simulated and predicted BRDF patterns at any azimuth and zenith angles. Figure 8 and Figure 9 illustrate the polar representations of the simulated and the predicted BRDF using the physical RTM and 1D-CNN. Both figures demonstrate the BRDF variation in a polar space for the four OLCI bands: Oa04, Oa06, Oa08, and Oa17. Polar angles were used to represent relative azimuth angles (RAA), ranging from −180°, to 180°, with a 12° interval. Simultaneously, the VZA is denoted by a radius ranging from 0° to 50°, with a 1.667° interval.

The predicted BRDF pattern using the 1D-CNN model matches with the simulated BRDF pattern using RTM base model. Indeed, the estimated and simulated BRDF in four bands vary from low-directional reflectance values (dark blue) to high-directional reflectance values (vivid red). When the sensor direction and sun coincide, directional reflectance values in the backscattering direction are remarkably higher than in the forward and nadir directions (Figure 8 and Figure 9). This is the most characteristic and prominent feature of the BRDF, known as the hot-spot effect, which is characterized in the proposed 1D-CNN by a peak reflectance at a view zenith angle of 30° near the sun position at all the illustrated polar BRDF. These results indicate the outstanding prediction ability of the 1D-CNN for BRDF characterization.

The directional reflectance values estimated by the 1D-CNN are plotted against the corresponding simulated directional reflectance values by PROSAIL+6S (Figure 10). The results show the developed 1D-CNN BRDF model can approximate the PROSAIL+6S in a near-perfect performance for the test dataset also. The coefficient of determination $R^2$ between the predicted and simulated reflectance values for Oa04, Oa06, Oa08, and Oa17 bands of the OLCI is over 0.97, indicating a good agreement between the simulated directional reflectance by PROSAIL+6S and the predicted directional reflectance using 1D-CNN.

Table 7. Comparison of metrics in four bands of the OLCI for testing datasets.

Spectral Bands	${\mathbf{R}}^2$	RMSE
Oa04	0.996	0.0013
Oa06	0.977	0.0044
Oa08	0.979	0.0030
Oa017	0.997	0.0018

summarizes the R² and RMSE. Similarly, compared to the validation dataset, the 1D-CNN showed a significant accuracy achievement for the testing dataset and confirmed the general ability of the proposed 1D-CNN to perform BRDF estimation.

Analysis of Estimation Error

Figure 11 illustrates the polar plot of the Percentage Error (PE) between the simulated reflectance values by PROSAIL+6S and the estimated directional reflectance values by the proposed 1D-CNN model. The red region indicates that more PE has been raised, particularly in the direction of the hot-spot near the principal plane, where the sun and sensor were closer to each other.

Moreover, there was a tendency for higher error in several relative azimuth angles when the view zenith angle increased. Regarding the polar plots of the PE at four of the OLCI bands, the minimum PE error obtained was relatively small and respectively equaled 1.72%, 2.12 % and 2.44% for NIR (Oa017), green (Oa06) and red (Oa08) bands, respectively. In comparison, a more significant error of about 3% arose for the blue (Oa04) band. The result indicated that the proposed 1D-CNN model was capable of characterizing BRDF patterns with high accuracy in the majority of azimuth and view zenith angles, except the hot-spot area, which slightly reduced the model’s predictive ability.

4.3. Comparison of the BRDF Shape in Principal and Cross-Principal Plans

The BRDF polar plots enabled us to compare the directional reflectance pattern at any azimuth and zenith angle, particularly in the Principal Plane (PP) and the Cross-Principal Plan (CPP) which are frequently used directions in BRDF analysis. Indeed, we could effectively analyze the proposed 1D-CNN’s performance compared to the RTMs-based model. The BRDF pattern as a function of VZA along the PP and CPP at four bands of OLCI appears in Figure 12 and Figure 13. The blue and red curves represent the estimated and simulated BRDF pattern using respectively the 1D-CNN and the physical RTM. Generally, both simulated and estimated curves form a similar curve influenced by the hot-spot effect in PP and bell-shaped in CPP. As shown, the 1D-CNN model appears to match the simulated BRDF with a slight difference in the hot-spot region and when the VZA increases. However, the results not only demonstrated a high degree of consistency between the simulated and predicted BRDF patterns but also confirmed the excellent ability of the developed 1D-CNN model to predict the BRDF behavior.

4.4. Application to Real Data

After downloading the required observations, a series of pre-processing steps were applied to build a data cube of high-quality, cloud-free reflectance-based observations. In the first step, the downloaded OLCI Level-1 data products in Full Resolution (300 m) were converted from radiometrically calibrated TOA radiances (W/m²·sr·μm) to TOA reflectance data and then projected the data into the geographic coordinate system (WGS 1984). In the second step, a range of filtering was applied based on the quality flag layer to mask the low-quality observations. In addition, all contaminated pixels by cloud and cloud shadow affecting quantitative analyses in the BRDF modeling were masked out using the quality flags (QF) band associated with the OLCI images. The QF is a 32-bit integer: bits from 0 to 20 indicate saturation per band and the rest of the bits present other quality information. The bit 27 (quality flags bright) related to clouds and cloud shadows was used to mask the contaminated pixels [80]. Figure 14 shows the collected OLCI data at the Oa17 band, from 1st July to 16th July 2019 before and after cloud and cloud shadow masking.

4.5. 1D-CNN Evaluation in the Real Domain

To evaluate the performance of the pre-developed 1D-CNN BRDF modeling in the real domain, the relationship between the estimated and the measured directional reflectance was calculated in the terms of $R^2$, and RMSE at four bands, including Oa04 (Blue), Oa06 (Green), Oa08 (Red), and Oa17 (NIR). The scatterplots and boxplots between the measured and the estimated directional reflectance across all four bands are shown in Figure 15. When comparing the measured reflectance with the estimated reflectance values at the Oa04 band (Blue), the $R^2$ shows a strong agreement ($R^2$ = 0.87, RMSE = 0.007). However, at Oa17 (NIR) band, there is a slightly weak relationship between the measured and estimated reflectance values (Oa17 band; $R^2$ = 0.75, RMSE = 0.007), showing a remarkably lower agreement around 13% than the Blue. In addition, the correlation at the Oa06 band (Green; R² = 0.79, RMSE = 0.01) and Oa08 band (Red band; $R^2$ = 0.83, RMSE = 0.015) is also lower around 9% than Blue. This can be because of the lower signal-to-noise ratio of OLCI images in the Green, Red, and NIR bands compared to the blue band or other scaling issues which make reconstructing the BRDF accurately more challenging. This finding indicates that our proposed network works better with less noisy data than with highly noisy data. This can be because of this reason that our network is developed in a simulation domain that is not affected by real noises.

5. Conclusions

This study presented an emulation-based approach to approximate the RTMs-based TOA BRDF model using a 1D-CNN. The RTMs (PROSAIL and 6SV) were used to simulate a large number of input-output values of geometry-surface-atmosphere parameters and corresponding TOA BRDF for training, validating and testing the 1D-CNN. The developed 1D-CNN reached a high accuracy and computational speed. The proposed approach also had an acceptable performance in the real domain, and based on a collection of about two weeks of OLCI observations. However, the proposed method is mainly developed for the vegetation area and does not apply to surfaces like water which has a different BRDF pattern. In addition, the model is applicable only for the defined ranges of inputs in the simulation domain and needs retraining in the case of changing these ranges for other areas of interest. Finding an appropriate 1D-CNN architecture and the optimal hyperparameters were based on a trial-and-error approach, which requires more time and effort to set. But once the network was trained and tested, it can work with acceptable accuracy and speed. Furthermore, the coupled surface–atmosphere RTMs for simulating the TOA BRDF could be furtherly improved by involving topographic effects to simulate heterogeneous rugged areas viewed from an off-nadir direction. Working at TOA level helps to integrate multi-sensor datasets (MODIS, VIIRS and OLCI) in future improvements of the method with a minimum pre-processing effort (no atmospheric nor topographic corrections required).More importantly, this effort proved that emulation-based approach can be an effective solution for establishing the relationship between the remote sensing observations and surface and atmospheric parameters, demonstrating the DL model’s ability to and emulate the complex environmental and physical relationships in remote sensing applications [47,48]. The deep learning-based approach could be applied to emulate even complex combinations of physical and/or empirical models, encountered in some remote sensing applications (ex.: hydrology, sediment transportation, climate change effects, life cycle). Deep learning can establish nonlinear and complex connections between input and output variables with advantageous accuracy and computation cost.

6. Future Works

Our future works will focus particularly on the application of the proposed methodology to the real domain (actual EO data) with multi-source, multi-temporal satellite data that are available for BRDF modeling through a variety of cloud computing platforms (CCPs). In the real domain, different CCPs, like Google Earth Engine (GEE), can provide various remote sensing data ranging from coarse to fine spatial resolution acquired from multi-source sensors. Additional works will investigate how to apply DNN based emulators to the atmospheric correction of satellite images using atmospheric-based RTMs such as 6SV or retrieving the vegetation traits based on TOA reflectance data (RTMs inversion mode). BRDF will be considered in these processes since the model has already been developed. An integrated radiometric correction including atmospheric, topographic and BRDF effects on satellite acquisitions is also an avenue to look at.