Atmospheric correction for HY-1C CZI images using neural network in western Pacific region

ABSTRACT With a spatial resolution of 50 m, a revisit time of three days, and a swath of 950 km, the coastal zone imager (CZI) offers great potential in monitoring coastal zone dynamics. Accurate atmospheric correction (AC) is needed to exploit the potential of quantitative ocean color inversion. However, due to the band setting of CZI, the AC over coastal waters in the western Pacific region with complex optical properties cannot be realized easily. This research introduces a novel neural network (NN) AC algorithm for CZI data over coastal waters. Total 100,000 match-ups of HY-1 C CZI-observed reflectance at the top-of-atmosphere and Operational Land Imager (OLI)-retrieved high-quality remote sensing reflectance (R rs) at the CZI bands are built to train the NN model. These reflectance data are obtained from the standard AC algorithm in the SeaDAS. Results indicate that the distributions of the CZI retrieved R rs were consistent with the quasi-synchronous OLI data, but the spatial information from the CZI is more detailed. Then, the accuracy of the CZI data for AC is evaluated using the multi-source in-situ data. Results further show that the NN-AC can successfully retrieve R rs for CZI and the coefficients of determination in the blue, green, red, and near-infrared bands were 0.70, 0.77, 0.76, and 0.67, respectively. The NN algorithm does not depend on shortwave-infrared bands and runs very fast once properly trained.


Introduction
As a new generation ocean color sensor entailing four spectral bands with central wavelengths of 460, 560, 650, and 825 nm, the Coastal Zone Imager (CZI) onboard the HY-1C satellite launched in 2018 has provided continuous observations of China's coastal waters. The corresponding details of the spectral response function (SRF) are shown in Figure 1. The CZI camera is a new wide-band remote sensor with a swath width of approximately 950 km and a spatial resolution of 50 m. The revisit period is 3 d, and the overpass time is approximately 10:30 local time. This sensor has been designed to support the research of coastal zone areas, which are the regions closely related to human activities and where strong landsea interaction occurs.
The traditional Atmospheric Correction (AC) methods are not suitable for CZI because of its band setting. The standard AC algorithms (Gordon 1997;Gordon and Wang 1994) in Sea-viewing Wide Field-of-view Sensor Data Analysis System (SeaDAS) use two nearinfrared (NIR) bands to estimate the contribution of aerosols at visible bands. This method is sufficient in the open ocean areas but is challenging to operate in coastal areas where the black ocean assumption often fails (Gordon and Wang 1994;Zhang et al. 2007). Consequently, an alternative AC algorithm was constructed for AC in coastal waters with non-zero NIR water-leaving radiances (L w ) by using Shortwave Infrared (SWIR) bands. In this scheme, the black ocean assumption may still be tenable owing to the strong absorption of water in SWIR, and the L w in the visible bands can be obtained by extrapolation (Wang and Shi 2007). Several modified algorithms based on this principle were subsequently developed, and improvements in accuracy, at least for selected validation areas, were achieved (Jiang and Wang 2014;Wang, Shi, and Jiang 2012;Wang and Shi 2007). However, as CZI only comprises the four bands mentioned above, the lack of SWIR renders the AC for CZI a daunting task.
The Neural Network (NN) method is a powerful tool for prediction, recognition, function approximation, and pattern classification (Haykin and Network 2004). The feasibility of applying the NN method to AC has been investigated in several previous studies (Brajard et al. 2012;Brajard, Moulin, and Thiria 2008;Chen et al. 2014;Fan et al. 2017;Hamre et al. 2013;Jin and Stamnes 1994;Li et al. 2020;Schroeder et al. 2007). Brajard, Moulin, and Thiria (2008) used a trained NN to replace the radiative transfer model in the forward simulation to obtain the Top-Of-Atmosphere (TOA) radiances. Fan et al. (2017) employed a coupled atmosphere-ocean radiative transfer model to simulate the radiance at the TOA (L TOA ) and the remote sensing reflectance (R rs ) just above the surface. They also trained two NNs to directly derive the aerosol optical depth and R rs from the L TOA , respectively (Fan et al. 2017). Although the studies demonstrated that the NN method could be successfully used for AC, it could hardly be applied to CZI images. For example, only a few in-situ data could be synchronized with CZI, hence the lack of representativeness of the training samples. An alternative solution is to use an auxiliary sensor, but its setting needs to be like that of the CZI in the visible bands. The auxiliary sensor should also possess NIR and SWIR bands to be able to produce R rs with the standard AC methods. The Operational Land Imager (OLI) onboard Landsat-8 generates high-quality aquatic science products (Franz et al. 2015;Pahlevan et al. 2017), that can be applied as an auxiliary sensor for CZI owing to the following reasons: (1) OLI has a wider spectrum coverage than CZI, including NIR and SWIR bands. (2) The CZI camera visible bands (central wavelength: 460, 560, 650, and 825 nm; bandwidth: 420-500, 520-600, 610-690, and 760-890 nm) roughly correspond to those of 3,4,and 5 (central wavelength: 482,561,655 and 865 nm; after band conversion, as shown in Figure 2. (3) The relatively high Signal-To-Noise Ratio (SNR) and the similar spatial resolution from OLI will be beneficial for the AC of CZI. (4) The training strategy, as to be described in the following paragraph, can help the model to further consider the differences (e.g. bandwidths) between two sensors.
This study proposes a new AC algorithm for CZI based on the NN, a scheme that differs from those in previous studies in terms of NN structure and training process ( Figure 2). First, we calculate the TOA reflectance (ρ TOA ) from the CZI data and R rs from the OLI  data by using the standard AC method (hereafter called as SWIR-AC) embedded in SeaDAS, which is one of the best AC methods for OLI (Xu et al. 2020). Then, as shown in Figure 1, we apply the band conversion of R rs from OLI to CZI based on the band SRF to build the simulated dataset for model training (Pinto et al. 2016). Based on the dataset, we take the ρ TOA of the four bands and observation geometries (Solar Zenith Angle (SOLZ), Sensor Zenith Angle (SENZ), and Relative Azimuth Angle (RAA)) as the inputs of NN and R rs as the outputs to train the network. Finally, we validate the NN method on a testing dataset and in-situ data. The detailed validation results are detailed in the succeeding sections.

Simulated dataset
For the model training and validation, a simulated dataset was built by selecting the synchronous remote sensing images of OLI Level 1 and CZI Level 1B for 2019 from https://earthexplorer.usgs.gov/ and https:// osdds.nsoas.org.cn/#/, respectively. The corresponding locations are marked in red boxes in Figure 3(a). A total of three match-up pairs of the images were obtained with the cloud coverage < 3%, and the time difference < 3 hours.
The top-of-atmosphere reflectance (ρ TOA ) from the top-of-atmosphere radiance (L TOA ) of CZI is calculated using where θ 0 and F 0 are the SOLZ and extraterrestrial solar irradiance, respectively. The observation geometries (cosine of the SOLZ, cosine of the SENZ, and cosine of the RAA) and ρ TOA at four bands were used as the inputs of NN. As for the outputs of NN, R rs for four bands were obtained from the corresponding OLI images via the SWIR-AC embedded in SeaDAS (ver. 7.5). The spectral response difference (Figure 1) between the two sensors was compensated by the spectral band adjustment factor β λ ð Þ, as calculated by Equation (2), using the relative response functions and field measurements R field rs (N= 83) obtained for the South China Sea in September, 2018.
where SRF is the spectral response function, and the β values at the 460, 560, 650, and 825 nm bands are 0.886211, 0.9995, 1.020512, and 1.240233, respectively. The spectral band conversion of R rs from OLI to CZI was accomplished using Equations (2) and (3). The OLI images were resampled to the spatial resolution of 50 m to ensure pixel homogenization. Quality control methods for Level-3 global products were used to ensure the quality of the simulated data. All pixels were screened using the 12-flags bit  (25). In the list, the numbers in the parentheses represent the bit position of the 32-bit flag value. Then, the following criteria were used to extract the high-quality R rs samples: (1) The percentages of pixels with valid R rs (460 nm) values in the 3 × 3 box are checked.
If the valid percentage is larger than 50%, then the data in the box are further checked by the following quality assessment, otherwise it is discarded.
(2) The means and standard deviations (SDs) of the valid R rs (460 nm) values within the box are calculated. Pixels with R rs (460 nm) values beyond the range of mean±1.5SD in the box are discarded.
(3) The means and SDs of valid R rs (460 nm) values are recalculated for the remaining valid pixels, and the coefficients of variation (CVs) are determined (i.e. SD divided by mean) to check for spatial heterogeneity. If the CV was is than 0.15, then the box is adopted in the next step; otherwise, it is discarded.
Standards (1) to (3) can help to ensure the uniformity of the spatial range of data and avoided the noise caused by instruments or stray clouds. In total, we achieved 2,268,327 match-up pairs. Among them, 100,000 pairs were randomly selected as the simulated dataset. Then, the selected pairs were divided into three independent groups: a training subset (70,000 data points), a testing subset (15,000 data points), and a validation subset (15,000 data points). The spectral shapes of the simulated dataset are shown in Figures 3(b,c), and the maximum, mean, and minimum values are shown in Table 1. In view of ensuring the representativeness of the training dataset, a simple empirical criterion was applied to classify the different types of water: R rs (650) < 0.0005 [sr −1 ] for clear water, R rs (650) > 0.012 [sr −1 ] for turbid water, and inbetween values of these two ranges for moderately turbid water. The rates of the three types of water in the training dataset were 79%, 17%, and 4% for clear water, turbid water, and moderately turbid water, respectively. The percentages indicate that the training dataset comprised various water types.

Neural network training
Previous studies have demonstrated that NN can approximate nonlinear functions (Chen et al. 2014;Liu et al. 2021;D'Alimonte, Zibordi, and Berthon 2004;D'Alimonte and Zibordi 2003;Géron 2019;Cybenko 1989;Hornik, Stinchcombe, and White 1989;Pinkus 1999). Therefore, the NN is suitable for solving our inverse problem of deriving R rs from ρ TOA . This section describes the NN model designed for the CZI sensor.
In the construction of our NN, an important issue was to find the optimum number of hidden layers and neurons. The choice generally depends on many variables, but this approach hinders the identification of the best solution in many cases. We used a practical approach to determine the number of hidden layers and neurons by utilizing many NN configurations with different numbers of hidden neurons (In particular, 10, 11, 13, and 15 neurons were selected, and we did not add more hidden layers due to the high-accuracy objective for the testing dataset). We recorded the process and outcomes of these NNs as a means of determining the optimal performance and computing time configuration for the training process. Finally, we used a single-hidden layer configuration with 11 neurons, in which 0.99 was obtained for the coefficient of determination (R 2 ). Our NN architecture is shown in Figure 4. The input layer entails seven elements, including the cosines of SOLZ, SENZ, and RAA, and the ρ TOA at the four wavelengths of 460, 560, 650 and 825 nm. The output layer entails four elements corresponding to the R rs values of the four wavelengths mentioned above. The input and output parameters were normalized using the mapminmax function embedded in the MATLAB NN toolbox, as follows: where x and y are the values of the training data before and after normalization, and x min and x max are the minimum and maximum values of x. A hyperbolic tangent sigmoid function was utilized as the activation function of the neurons, and a linear function was employed to transfer the hidden layer neurons to the output. Moreover, we adopted the Levenberg-Marquardt backpropagation algorithm for the training process of the NN. In the training process, the simulated dataset with 100,000 match-up pairs, was divided into model training, validation, and testing datasets based on the approximated random selection of 70%, 15%, and 15%, respectively. The training dataset was used to optimize the weights and biases of the NN according to the error. The validation dataset was used to measure the network generalization and to halt the training when the generalization shows no apparent improvement. The testing dataset was not involved in the training process, because it would be used for the unbiased evaluation of the performance of the trained network. The initial weights and biases of the network were randomly obtained with values from 0 to 1 (Nguyen and Widrow 1990). Then, the training algorithm was used to iteratively update the weights and biases in accordance with the training error between the simulated R rs values and the NN outputs. In each iteration, the validation dataset was used to monitor the performance of the current NN by computing the mean square error (MSE). If the MSE increases by more than ten iterations, then the training process is completed; otherwise, the training stops after 1,000 iterations. The relatively uncomplicated model can also help to prevent over-fitting.

In-situ data
To test the performance of the proposed AC method, we collected in-situR rs data from different sources, including (1) station observation data in western Pacific region, such as the Socheongcho, ARIAKE_TOWER, MuPing, and DongTou stations, and (2) ship-borne measurements data in the Bohai Sea. The red, green, and yellow dots in Figure 3 represent the specific locations of the different observation sources. The following aspects were considered in the match-up determination: time difference of < 3 hours, valid pixels of >50% in the 3 × 3 box after quality control, and CV of < 0.15 (Section 2.1). After removing some outliers and data without CZI matches, there are 65 match-up pairs left. The mean R rs at the 460, 560, 650, and 825 nm were 0.0133, 0.0172, 0.0065, and 0.0009 sr −1 , and the corresponding CVs were 0.8744, 0.6718, 0.9697, and 1.080. The detailed in-situ data information is discussed in the succeeding sections. Figure 4. Architecture of the neural network that models the atmospheric correction for the HY-1 C CZI data. The SOLZ, SENZ and RAA is cosine of the solar zenith angle, sensor zenith angle, relative azimuth angle, respectively.

Station observation data
Aerosol Robotic Network-Ocean Color (AERONET-OC), equipped with additional deployments of autonomous radiometer systems on offshore platforms, provides consistent and accurate measurements to support satellite ocean color validations (Zibordi et al. 2009). Two AERONET-OC stations (Socheongcho and ARIAKE_TOWER stations) were selected to obtain Level 2.0 data from NASA (https:// aeronet.gsfc.nasa.gov/new_web/data.html). The two stations represent different aquatic optical properties. The Socheongcho station in the Yellow Sea is located far away from land with clear surrounding water, while the ARIAKE_TOWER station in the Sea of Japan is located close to land with surrounding turbid water. As these stations have long-term observations, they offer a superiority alternative approach of validating the time consistency of the NN algorithm's performance. All stations use the same instrument, including calibration and post-processing procedure, indicating that discrepancies introduced by the differences in the measurement methods can be avoided. The normalized water-leaving radiance (L wnÀ f=Q), which was corrected for non-isotropic effects, was used in this research. The L wnÀ f=Q data at the wavelengths corresponding to the CZI were subsequently converted into R rs by using the same method employed in the literature (Pahlevan et al. 2017). Due to the lack of corresponding bands for the CZI, a linear interpolation method was used to realize the spectral band matching between satellite data and measured data (Li et al. 2020). 443, 490, 551, and 667 nm were used for interpolation to obtain R rs (460), R rs (560), and R rs (650). Due to the lack of the AERONET-OC band corresponding to 825 nm, only three visible bands were compared.
MuPing and DongTou stations are in the Yellow Sea and the East China Sea and they were designed by the National Satellite Ocean Application Service (NSOAS) to obtain long-term in-situ values of apparent optical properties, inherent optical properties, atmospheric optics, hydrology, meteorology, and main components concentration of ocean color. Both stations can provide high-frequency in-situR rs (i.e. nearly half an hour per observation) as measured by a CE-318 spectroradiometer. The data processing scheme follows the NASA ocean optics protocol (Mueller et al. 2003). 442, 490, 620, 667, 779, and 865 nm were used for linear interpolation to obtain R rs (460), R rs (650), and R rs (825).

Ship-borne measurements data
The ship-borne measurements data used in this study were measured by an autonomous ship-borne abovewater hyperspectral radiometer (350-2500 nm with ~1 nm increments) during a cruise on 29 October 2018 in Bohai Sea. The data processing scheme followed the NASA ocean optics protocol (Mueller et al. 2003). The total distance of the ship was 5.4 km and 9 sample points were measured on the way within four hours starting at 10:00 am local time. The values of the sampling points were collected in a short time, indicating the unique superiority for validation of AC methods.

Accuracy assessment
The accuracy of the R rs prediction was assessed by comparing the retrieved R rs values with the true values obtained from the simulated dataset and in-situ data. The comparison was quantified using the mean Absolute Percentage Deviation (APD) and the Root-Mean-Square Error (RMSE) between the retrieved R rs and true R rs :

Cross-sensor agreement
We extracted the concurrent OLI and CZI TOA reflectance values within three hours to analyze the correlations for all four bands of CZI via linear regression. As shown in Figure 1, differences exist between the spectral coverage and center wavelengths of OLI and CZI in all four bands. The agreement and difference over CZI and OLI images, were visualized and quantified by randomly selecting and comparing 100,000 ρ TOA pairs of CZI and OLI As shown in Figure 5, although the data have some discretization, most of them were centered around the 1:1 line for all four bands, indicating a strong agreement between OLI and CZI. The positive linear trends can be used to describe the relationships between the CZI and OLI data. The slope values were extremely close to 1, and the R 2 values are higher than 0.79. The results depict a generally strong agreement between the two sensors.

Performance of the new NN algorithm with simulated data
The performance of our NN algorithm was initially evaluated with the simulated dataset. The training and testing datasets (70% and 15% of the simulated dataset) were used to evaluate the trained NN. The performance assessment was based on statistical parameters, including R 2 , APD, and RMSE. Figure 6 shows a comparison of the scatter plots of the R rs values retrieved by the NN-AC and the simulated R rs values obtained from the training and testing datasets. Figure 6(a) shows the performance of the NN on the training dataset. The retrieved R rs (λ) values from the OLI images are consistent with the simulated values. The R 2 values are larger than 0.98 at the blue, green, and red bands, and larger than 0.93 at the NIR band. The APD values of the training datasets are lower than 15% at the visible bands, and the RMSE values do not exceed 0.0006 sr −1 for all bands. More importantly, the results obtained from statistical parameters on the testing datasets are extremely close to those of the training datasets, as shown in Figure 6(b).
The comparison of the retrieved R rs values and the true values on the training dataset indicate that the newly built NN-AC method can be used to simulate the process of radiative transfer. The homologous results on the testing dataset indicate that the NN-AC can accurately learn the training dataset.
To examine the impact of uncertainty, we added a 3% random error to the ρ TOA of the testing dataset. The results are shown in Figure 7. An increase in APD from 8.65%, 4.14%, 14.95%, and 13.59% to 18.45%, 31.53%, 35.66%, and 45.25% can be observed for the four bands, respectively. Regardless of the high scattering, the scatter points are almost within the range of 1.5 SDs at the visible bands. In the NIR band, the APD increases by 30% due to the random error; this trend may be caused by the relatively low signal and SNR values. Overall, the NN-AC can resist the influence of random error to a certain extent. However, when the signal values are extremely low, the instability of the method is apparent.

Evaluation of R rs retrieval with CZI data
The performance of the NN-AC was also evaluated by applying the algorithm to an HY-1 C CZI image in an optically complex coastal area. Taking the Yellow Sea as an example, as shown by the blue box in Figure 3. This region has relatively clean water and turbid water. Figures 8(a,b) show the contrast between the NN-AC and SWIR-AC for the turbid water and clear water areas, respectively. The retrieved R rs values by NN-AC were extremely close to the results of SWIR-AC for two types of water with R 2 higher than 0.97 and RMSE lower than 0.0085. The results indicate that the NN-AC method is suitable for water bodies of varying degrees of turbidity. Figures 9(a-h) show the highly similar spatial distributions of the R rs values retrieved from NN-AC and SWIR-AC. Many black points can be observed in Figures 9(e-h) presenting the AC failure of SWIR-AC; this problem does not exist in the NN-AC results. The details are magnified in Figures 9(a,b,e,f) for easy viewing. The comparison indicates that the NN-AC method can achieve the AC of CZI, and its spatial distribution of R rs is smoother than that of SWIR-AC. Figures 9(i-l) show the APD spatial distributions at the four bands for the NN-AC and SWIR-AC, while Figure 9(i) shows the mean APD histogram. The APDs at the first three bands are stable at approximately 20%. The relatively small error for the high-turbidity water indicates the NN-AC method can adapt to coastal waters. Notably, the APD in the NIR band is close to 40%, which may be caused by the low signal in this band; this finding indicates that the NN-AC method is unstable in clean water with low signal. The uncertainty caused by the band conversion may also lead to a relatively high APD, as depicted by spectral distance of approximately 40 nm in the central wavelength of the NIR band between the OLI and CZI data. The CZI's bandwidth is nearly twice than that of OLI.

Validation of the NN algorithm with in-situ data
The accuracy of the established algorithm was further validated using the multi-source in-situ data. The detailed information about the match-up selection has been previously presented in Section 2.3. Figure   and 0.69, respectively, whereas those in turbid water (ARIAKE_TOWER and Ship-borne measurements) are 0.69, 0.70, 0.59, and 0.45, respectively. The results indicate that the method proposed in this study can perform better in clear water and achieve competitive performance in turbid water because of the good consistency with the in-situ data, especially at the visible bands. More importantly, the results show that the  NN-AC algorithm can accomplish the AC for CZI in coastal waters regardless of the spatial or temporal difference.
The accuracy of the NN algorithm is limited by the original algorithm. This limitation suggests that the issue may be solved by improving the AC algorithm used for OLI images. However, a much higher precision training dataset is difficult to obtain, as complex algorithms are required to handle the AC of water with complex optical properties. As a probable solution, instead of training a single gigantic NN for handling different ocean conditions, the NN-AC can utilize a targeted dataset representing different types of atmospheric and marine conditions for model fine-tuning. Spatial location and data statistical analytic technique may be introduced when selecting the NN that can best match the marine conditions based on the geographical coordinates and the satellite measurements of ρ TOA to yield the optimum retrieval results.

Conclusions
The implementation of AC of CZI images over coastal waters is challenging due to the lack of NIR and SWIR bands. In this study, we proposed a new AC algorithm based on the NN for CZI data to solve the problem. Different from the approach of traditional AC algorithms that use a small volume of in-situ data, we generated numerous high-quality simulated matchups from CZI and OLI images for the NN training.
The new NN algorithm was validated using the simulated dataset and the results showed that it could accurately and effectively learn the training dataset. Furthermore, the algorithm attained a similar spatial distribution pattern and even a higher spatial continuity of R rs when evaluated with the retrieved R rs from the OLI image by using SeaDAS. Then, we further validated the algorithm by multi-source in-situ data, and the results showed that the NN algorithm could accomplish the goal of AC for the CZI sensor in coastal waters regardless of the spatial or temporal difference. The NN algorithm runs very fast once the network has been properly trained, and is therefore, suitable for operational use.
Meanwhile, there are some deficiencies in this study. Our method can sufficiently simulate the results of SWIR-AC, but some errors still exist between retrieved R rs and in-situ data. If the quality of the training data can be further improved in the future, then the performance of this algorithm can also be enhanced. Jianqiang Liu received the MS degree in marine meteorology from the National Marine Environmental Forecasting Center, State Oceanic Administration, Beijing, China, in 1989. He is the Chief Designer of the ground application system for the Chinese-France ocean satellite and new generational ocean color satellites and was a Deputy Chief Designer of the ground application system for the HY-1 and HY-2 satellites. He is one of the founders of satellite ocean remote sensing in China and plays an important role in the development of Chinese Ocean Satellite and manned space flight. His research interests include ocean color and sea wave satellite data processing, applications of ocean remote sensing data and Antarctic research.
Guangping Xia is currently pursuing her MS degree in Photogrammetry and Remote Sensing at LIESMARS, Wuhan University. She is interested in ocean color remote sensing.
Tong Yue is currently pursuing her BS degree in spatial information and digital technology at Wuhan University. She is interested in ocean color remote sensing.
Ruqing Tong is currently pursuing her MS degree in Photogrammetry and Remote Sensing at LIESMARS, Wuhan University. She is interested in ocean color remote sensing.
Liqiao Tian is a professor in LIESMARS, Wuhan University. He received the MS and PhD degrees in Cartography and GIS from Wuhan University. His current research interest is the application of remote sensing.
Kohei Arai received the BS, MS, and PhD degrees from Nihon University in 1972University in , 1974University in , and 1978. He moved to Saga University as a Professor in the Department of Information Science in April 1990. Since 1998, he has been appointed as an Adjunct Professor of the University of Arizona, Tucson. He is Science Calibration Team Leader for the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) project. His major concerns are atmospheric optics, vicarious calibration, nonlinear optimization theory, inverse problem solving, and pattern recognition/understanding.

Linyu Wang received the MS and BS degrees in Land
Resource Management from China University of Geosciences (Wuhan), in 2014 and 2017. He is interested in the application of remote sensing.