Evaluation of an Extended PICS (EPICS) for Calibration and Stability Monitoring of Optical Satellite Sensors

...........................................................................................x

Vuppula [16] presented a new technique known as the "PICS Normalization Process" (PNP) that combines OLI observations of multiple PICS into a single time series dataset with greater temporal resolution. She used the OLI data from six PICS (Libya-1, Libya-4, Sudan-1, Niger-1, Niger-2, and Egypt-1) and normalized them to Libya-4. The temporal resolution was increased to approximately 4 to 5 days compared to the 16 day OLI revisit time. Unfortunately, this combination method could not guarantee the generation of a dataset with daily or nearly daily acquisitions.
In their paper, Shrestha et al. [17] wanted to identify "optimal" regions that were common across all image bands. They used an unsupervised classification technique to generate a set of 19 classes or "clusters" of spectrally similar OLI image pixels of North used to estimate the overall temporal and spatial variation in radiometric measurements over a major Cluster 13 sub-region (~40% of total area covered by Cluster 13 pixels).
In addition, OLI data from Clusters 1, 3, 4, 16, and 19 were analyzed to evaluate sensors across their wider operating dynamic range.
The paper is presented as follows: Section 1 provides a brief review of the topic and previous research performed in this area. Section 2 describes the dataset, mask generation process and application approach in greater detail. Section 3 presents produced results and critical analysis of the results when this process is applied to a number of different sensors. Finally, Section 5 offers some concluding remarks.

SDSU Processed Google Earth Engine (GEE) Derived Data and Mosaic of North Africa
In order to develop an EPICS, Shrestha et al. [17] used GEE [18], which utilizes Google's worldwide processing and storage resources to archive and process freely available image data [19] from Landsat and other satellite sensors, to produce the analysis dataset in their initial stage of work. Inside the GEE environment, OLI was chosen as it has an established calibration accuracy within 3% [10,20]. Temporally filtered and down sampled statistics image datasets (collection of 25 band Image containing the pixel wise temporal mean, temporal standard deviation, temporal uncertainty and the number of valid pixels used to generate those statistics) for the VNIR, SWIR, and Cirrus multispectral bands were retrieved from GEE as 1° latitude by 1° longitude georeferenced chip files. The chips were locally mosaicked into a continental scale, 300m spatial resolution image of North Africa that covered an area between 36°N to 15°N latitude and 18°W to 35°E longitude. A detailed explanation on the development of the analysis dataset can be found in [17]. This work will take advantage of the already developed EPICS to evaluate its usability for calibration and stability monitoring.

Classification Map of North Africa
As mentioned in the introduction section, Shrestha's unsupervised classification algorithm ran on the cloud-screened and temporally filtered mosaic image (mentioned in Section 2.1) identified 19 distinct clusters of spectrally similar surface cover. Figure   6 1 shows a map of the identified clusters and their distribution throughout North Africa.
The bright green pixels represent Cluster 13.

Figure 1.
North Africa cluster map.

Cluster 13 as EPICS Candidate Cluster
Shrestha's original Cluster 13 was found to have an estimated spatial uncertainty (i.e. the ratio of spatial standard deviation to spatial mean of the mosaic filtered pixels) within 3% in all but the Coastal/Aerosol and Blue bands, which had estimated spatial uncertainties within 5%. Table 1 provides the estimated mean TOA reflectance, spatial uncertainty and average temporal uncertainty of Cluster 13 calculated from the 300m resolution mosaic filtered cluster pixels. Although in the original algorithm the maximum temporal filtering threshold was set at 5%, the unsupervised classification algorithm by Shrestha resulted in less than 3.35% temporal uncertainty across all bands for Cluster 13. The detailed process for determining the uncertainty values will be found in [17]. In addition, Cluster 13 contains a large number of pixels which are aggregated into contiguous regions. Due to its greater degree of pixel aggregation, larger pixel counts, and lower temporal and spatial variability in TOA reflectance measurements, EPICS analyses described in this paper focused on this particular cluster.

Creation of Cluster 13 Zone-Specific Masks
To aid and speed up the retrieval of Cluster 13 TOA reflectances from native resolution satellite images, raster masks were created using the boundaries defined in 8 the KML vector file which matched up with the sensor's specific spatial resolution. The masks where generated at UTM-zone size, which was to improve efficiency of processing. This section describes the mask creation procedure in greater detail.
First, the KML polygon lat/lon coordinates were converted to binary masks whose pixels were registered to the corresponding UTM map projection coordinates; the resulting geo-referenced masks possessed a spatial resolution matching, in this case, the 30 meter resolution Landsat images. Images potentially crossing multiple UTM zone boundaries were accounted for by oversizing the mask dimensions by 1.5°, which also resulted in more efficient processing. This procedure is applicable for any sensor as long as the masks are generated to match the sensor's spatial resolution. Figure 2 shows two Cluster 13 pixel masks generated with respect to UTM zones 29 and 34, which are shown in grey. The red outlines indicate the cluster region boundaries, the white pixels inside the mask correspond to valid Cluster 13 pixels, and the blue parallelograms represent the footprint of Landsat WRS2 path/rows (considered in additional detail and 34 (shaded region on right).
in Section 3). In total, seven masks were needed to represent the entire Cluster 13 region across North Africa, with each mask covering an area of approximately 2,500,000 km 2 .

Application of Cluster 13 Zone-Specific Masks
The generated Cluster 13 zone specific masks from previous section and quality control/assessment masks (which flags each pixel as clear, clouds, water/snow/ice, fill, etc. and are generally provided with satellite scenes), were applied to each image, and all resultant good pixels were collected. The spatial mean and spatial standard deviation of TOA reflectances were calculated using all the good Cluster 13 pixels of all the scenes available in the local archive (which is a local cache of all good images from Landsat and Sentinel mission, containing less than 5% cloud cover). The process was This shows how the UTM zone mask is applied on the satellite imagery. In this case, the Landsat-8 scene is ready to be masked out by the Cluster 13 pixel mask and QA band. Figure 3b shows a filtered TOA reflectance image (magnification applied for better visual representation) for the Coastal/Aerosol band image where the gray level pixels correspond to valid Cluster 13 TOA reflectance pixels.

Additional Data Filtering
Original Cluster 13 [17] was generated from the filtered 300 m spatial resolution mosaic, with the summary statistics in all input data constrained to a maximum of 5% temporal uncertainty. Table 1 shows that the resulting spatial uncertainties of Cluster 13 across all bands are all within 5%. Consequently, the spatial uncertainty in the native resolution image statistics after application of the georeferenced zonal pixel masks was also expected to be within 5%. So, if an individual data point in a particular band of the Cluster 13 time series dataset exhibited an estimated spatial uncertainty above the 5% threshold, or the data point appeared to deviate significantly from the overall TOA reflectance trend, the corresponding source image was considered "suspect" with respect to clouds/shadows or other conditions not identified in the Quality Assessment (QA) band filtering. The image was then inspected visually; if a previously unidentified cloud/shadow or other artifact was observed, the entire scene's statistics were excluded from further processing.

Development of Cluster-Based EPICS BRDF Model
The Cluster 13 time series trends from the different sensors exhibit variability in the TOA measurements, especially in the longer wavelength bands. Several factors can affect this variability including seasonal atmospheric aerosol/water vapor changes. The most significant contributor to this seasonal variation is BRDF effects [22] due mainly to solar position change. The angular dependencies were normalized using the procedure described below.
Image products for the sensors analyzed in this work include information of perpixel values of the solar and sensor zenith and azimuth angles (For Landsat-8 OLI, per pixel angle band information are in the metadata). Recall from Section 2.6, the georeferenced pixel masks used to generate the Cluster 13 reflectance statistics were also applied to the associated per-pixel angle images to calculate average values for corresponding sensor viewing and solar zenith and azimuth angles. This information was then used to develop a four-angle model which is based on Farhad's [23] procedure.
As the modeling and manipulation of data using computer software (MATLAB) favors the use of Cartesian co-ordinate system and the sensor view and solar angles in the angle band image are given with respect to a three-dimensional spherical coordinate system, the first step in his model generation was to project the angles into a twodimensional Cartesian coordinate space. Kaewmanee [24] extended Farhad's linear model with quadratic terms for both solar zenith and solar azimuth angles, including an additional potential interaction term between them. For the purposes of this analysis, Kaewmanee's approach has been further extended by including all possible interaction terms between the sensor view and solar angles and quadratic terms for the sensor view angles.

2
(1) where is the model predicted TOA reflectance, 0−14 are the model coefficients, and x1, y1, x2 y2 are the Cartesian coordinates representing the planar projections of the solar and sensor view angles originally given in spherical coordinates where SZA, SAA, VZA, and VAA are the solar zenith/azimuth and sensor viewing zenith and azimuth angles, respectively.
All terms were assumed to be required for effective characterization, as a cluster can contain pixels from widely separated regions possessing distinct, and generally unknown, BRDF characteristics. It has been found that, some terms in the above model Here, is the observed mean TOA reflectance from each scene and is the model predicted TOA reflectance. For this analysis, was set to the mean TOA reflectance of the respective time series.

Cluster 13 Imaging Frequency
The

Cluster Optimization
In order to minimize local processing and storage demands, one WRS-2 path/row pair was selected for each Landsat cycle day, such that the image contained the largest number of Cluster 13 pixels. Based on these criteria, nine paths were excluded from the initial analysis. In addition to the set of useable paths/rows on each cycle day, Table 2 shows the resulting single path/row for each cycle day, the corresponding geographic area (in km 2 ) covered by the path/row image, and the total intersecting Cluster 13 region pixel count. The individual path/rows represent the selected "optimized" path/row (shown in the purple boxes in Figure 2 covering approximately 40% of the total Cluster 13 area). Furthermore, Table 2 shows the additional path/row pairs that could be used if a cloud-free acquisition of the optimized path/row area was not available. In this paper, additional path/rows were not considered due to storage limitation and processing optimization.   without any further cloud screening and correction (Right).

Traditional PICS vs. EPICS
The key advantage of cluster-based calibration method is the ability to perform daily/near daily evaluation of a sensor's stability and calibration. In Figure 4, 1434 cloud-free OLI scenes of Cluster 13 (using optimized path/rows), acquired since launch to August 2018, were used to generate Cluster 13 TOA reflectance time series. The However, the spatial uncertainty of Cluster 13 is higher than the spatial uncertainty of Libya-4. This is expected as Libya-4 was chosen specifically to reduce variability by finding a "very" homogenous region, where as Cluster 13, allowed for more variation (5% spatial uncertainty criterion that was set in the Shrestha's classification [17]), in order to achieve greater spatial extent. Some of the extra spatial uncertainties are also due to the variation of solar and view geometry within the individual scenes in the current cluster-based analysis.

Cluster 13 Region Similarity
Recall that clustering algorithm mentioned in [17] ensured spectral similarity of Cluster 13 pixels to within a temporal and spatial uncertainty of 5% (mentioned in Section 2.3). An analysis was performed to determine whether images from the individual optimized path/rows exhibited similar spectral behavior within the initial uncertainty. For this analysis, cloud-free images of the optimized path/row regions were processed as described in Section 2.6 to determine the summary statistics (temporal 20 mean, temporal uncertainty, and average spatial uncertainty) for the valid Cluster 13 pixels. The TOA reflectance information from each path/row pair for each cycle day was normalized for BRDF effects as described in Section 2.8, then checked to ensure overall spectral similarity to within 5% uncertainty. Figure 6 a-g show the resulting plot of each path/row's temporal mean reflectances, temporal standard deviations and average spatial standard deviations for all the bands. For the purposes of this analysis, the "site" label on each plot's horizontal axis is a short-hand notation representing the "optimized" path/row as indicated in Table 2. The estimated temporal standard deviation and total standard deviation due to combined spatial and temporal uncertainty are represented by error bars with smaller and larger caps, respectively. For this analysis, the total uncertainty estimate assumes independence between the temporal and spatial uncertainties: where utotal is the propagated uncertainty due to temporal and spatial uncertainty, utemporal is the temporal uncertainty and uavg spatial is the corresponding mean spatial uncertainty. Figure 6 a-g reveal that most of the site's total standard deviation lies within ±5% of the overall Cluster 13 mean TOA reflectance. However, site 5 in Figure   6a,b and site 12 & 13 in Figure 6f reveal that the deviation of these site's mean reflectance from Cluster 13 mean is closer to 5%. It means that their contribution to Cluster 13 temporal variability is higher than other sites. This phenomenon suggests that if one of those sites is specifically selected for Cluster 13 behavior estimation, it could underestimate or overestimate the Cluster 13 mean behavior. These relative higher deviations of some site's mean TOA reflectance compared to Cluster 13 mean raises the question "How many random sites within Cluster 13 are required for 21 calibration and stability assessment of a sensor?" An answer to this question is presented in the next section.

Expected Behavior of Random Cluster 13 Location
As mentioned earlier, classification algorithm from [17] used to identify the various clusters assumes all data points within a given cluster exhibit the same general spectral behavior. Assuming this spectral similarity, any randomly chosen location within a cluster should, in principle, be able to serve as a source of a representative dataset for the entire cluster. An analysis to test this hypothesis was performed using image data acquired over the 16 sites (recall from Table 2 that "site" is used to indicate the optimized Cluster 13 WRS-2 path/row area imaged on a given cycle day), with the goal of determining the minimum number of sites required to achieve a specified uncertainty in the estimated mean reflectance.
Time series (TS) reflectance datasets, corrected for BRDF effects as described in Section 2.8, were generated for each site. The overall mean TOA reflectance was calculated using all possible distinct combinations of multiple sites (i.e. the overall mean TOA reflectance of the time series for two distinct sites, three distinct sites etc).      Looking at the 3-site combination in x-axis of Figure 9, the absolute differences for all bands tend to reach zero and the distribution standard deviations are all within ±0.01 reflectance unit. This case is mentioned in the previous paragraph which produces an uncertainty of ~1%. It suggests that prediction of the Cluster 13 TOA reflectance's within 1% uncertainty is possible using only three sites. Again, from Table 4, it is predictable that using single site will provide the highest uncertainty (the worst-case scenario) which is around 2% across coastal and blue bands. All these results suggest that the Cluster 13 mean TOA reflectance can be estimated within an uncertainty of maximum 2% using only one site and within 1% uncertainty using only three sites leaving a choice to trade-off between accuracy and number of sites selection.

Validation
As mentioned in the Introduction section, Shrestha's original clusters were generated from the GEE derived Landsat OLI dataset which-i) lacked BRDF correction for Sun and sensor geometry variations; and ii) required down sampling the original image data to 300 m resolution, so Cluster 13 temporal trending was evaluated using OLI, Landsat-7 ETM+, and Sentinel 2A/2B MSI image data at their native spatial resolutions, including BRDF correction as needed. Figure 10 shows the cluster-level temporal trend for all OLI bands after applying the BRDF correction described in Section 2.8. The corresponding statistics are presented in Table 5. Additionally, temporal mean values, temporal uncertainty values and temporal revisit intervals are also mentioned at the top portion of the figure. The temporal uncertainties of Cluster 13 lie within 3% across all the bands which suggests that Cluster 13 is as stable as traditional PICS sites while offering much larger ROI which also results in minimal infuse for any specific localized land change or variability. Some of the bands such as Green, NIR, and SWIR1 bands are extremely stable -less than 1.5 %. Again, some seasonality is still left in the SWIR2 channel producing relatively higher temporal uncertainty compared to its nearby longer wavelength bands which needs to be further investigated. Overall, temporal uncertainties from the optimized path/row pairs are less than the Cluster 13 temporal uncertainties predicted from [17] (also shown in Table 1) except for Coastal/Aerosol and Blue bands. This reduction in uncertainty values was expected because the BRDF effects in longer wavelengths were minimized in the trending showed in Figure 10 whereas corresponding values from [17] lacked this compensation. Sentinel 2B MSI image tiles were selected such that a significant portion of the optimized WRS2 path/rows were included within it, as listed in Table 6. Similarly, Landsat 7 ETM+ scenes were also selected by optimizing its intersection with Cluster 13 (i.e. same as OLI). But, only nine path/rows of Landsat 7 ETM+ were selected due to data storage limitations. The spectral differences between the Landsat OLI, Sentinel 2A MSI, Sentinel 2B MSI, and Landsat 7 ETM+ were compensated by applying spectral adjustment factor (SBAF) to Sentinel and Landsat 7 data sets. The SBAF correction process can be found in [25]. This compensation was done to ensure a better comparison between the Cluster 13 TOA reflectance measurements from different sensors.  Figure 11. Validation of Cluster 13 mean temporal TOA reflectance values using OLI, Sentinel 2A/2B MSI, and ETM+ Sensors.

Extension of Dynamic Range Using Lower Reflectance Clusters
The 19 clusters from Shrestha's analysis have their own distinct spectral signatures which provide a wide dynamic range for sensor calibration especially at longer wavelengths as shown in Figure 12.

Increase of Sensitivity to Detect Change in the Sensor
This section describes how the temporally rich dataset from Cluster 13 can help to detect "sensor change" quicker than traditional PICS. Due to the harsh environment in space, sensor response often decays with time and the reflectance trend starts to produce nonzero slope. Because of temporal variability present in the time series, the magnitude of the slope needs to exceed a certain minimum value, which is determined by computing the statistically significant minimum detectable trend.
The following equation in Weatherhead et al. [27,28] allows estimation of the required number of years (N) to detect a trend of magnitude at the 95% confidence level and with 50% probability: where, σ is month-to-month variability and ∅ is the 1-month lag autocorrelation in the time series. This equation can also be used to estimate minimum detectable trend when N (in years) is given as input. However, it was necessary to make sure that the units of σ and are the same in the equation.
Using this equation, Bhat et al. [11] [27,28], σ in the above equation was expressed as month to month variability (in percentage) by dividing the overall standard deviation of the monthly averaged TOA reflectance's by the overall mean values. Similarly, N is computed in terms of years by dividing the total number of months by 12. One-month lag autocorrelation coefficient calculation process is described as following.
Autocorrelation (also known as serial correlation) is a statistical method which is widely used in time series analysis to detect non randomness in a dataset by measuring the correlation of a signal with a delayed copy of itself. It is often expressed as a function (autocorrelation function) of the time lag between the two copies of signal.
Mathematically, the autocorrelation function measures the correlation between and + , where k = 0, 1, 2, ... K are time lags and is a stochastic process. According to [29], the autocorrelation for lag k is: where, is the estimate of the autocovariance and 0 is the sample variance defined In the above equations is the the sample mean of the time series defined as: and, N is the number of observations present. As both Cluster 13 and Libya-4 time series is converted to a monthly sampled time series, the autocorrelation at lag k=1 is the one-month lag autocorrelation.    Figure 10, it is visible that the temporal trend of green band after BRDF correction has very less temporal variation while the SWIR2 channel has the largest variations. As the Equation (8) used to calculate minimum detectable trend also depends on the variability of time series, Green band takes less amount of time and the SWIR2 channel shows less improvement using Cluster based method.

Summary and Conclusion
This paper focuses on the application of EPICS (Cluster) for stability monitoring of optical satellite sensors. EPICS provides a significant improvement of the temporal revisit period of calibration time series in contrast to the temporal revisit period offered 39 by traditional PICS with similar or 1%~2% higher temporal uncertainties depending on bands. One of the clusters, Cluster 13 (using limited regions and cloud-free scenes only), offers temporal revisit period of potentially as good as 1.4 days for Landsat 8 OLI in contrast to an average of every 18-20 days obtained from traditional PICS. By using all the regions of Cluster 13 and a nominal cloud consideration (~around 30% scene rejection due to cloud cover) a temporal revisit period of 0.33 day (~three cloud free collects everyday) can be obtained by a moderate resolution sensor. Furthermore, for the sensors having a wide field of view, Cluster 13 can offer even less than this revisit period. This improvement in the temporal revisit period resulted in better (depending on bands the increase in sensitivity as large as ~2 times) sensitivity of drift detection.
The temporal uncertainty of Cluster 13 was analyzed and validated using Landsat 8 OLI, Landsat 7 ETM+, Sentinel 2A MSI, and Sentinel 2B MSI. All sensors agree that Cluster 13 is temporally stable within 3%. However, spatial uncertainty of Cluster 13 is around 5% which is slightly more variation than that of traditional PICS. This extra spatial variation is pronounced due to its large spatial extent across the continent. Even though Cluster 13 with its hugely extended target size has spatial uncertainty of 5%, still the goal of matching traditional PICS temporal uncertainty is achieved. The near daily/daily calibration opportunity and increased sensitivity of change detection outweighs traditional PICS based methods, while also reducing the impact of a single location "change".
This paper also suggests that a single random location from Cluster 13 can be a