Accuracy Assessment of Atmospheric Correction of KMSS-2 Meteor-M #2.2 Data over Northern Eurasia

Abstract

Atmospheric correction of satellite remote sensing data is a prerequisite for a large variety of applications, including time series analysis and quantitative assessment of the Earth’s vegetation cover. It was earlier reported that an atmospherically corrected KMSS-M (Meteor-M #2) dataset was produced for Russia and neighboring countries. The methodology adopted for atmospheric correction was based on localized histogram matching of target KMSS-M and MODIS reference gap-free and date-matching imagery. In this paper, we further advanced the methodology and quantitatively assessed Level-2 surface reflectance analysis-ready datasets, operatively produced for KMSS-2 instruments over continental scales. Quantitative assessment was based on accuracy, precision, and uncertainty (APU) metrics produced for red and near-infrared bands of the KMSS-2 instrument based on a reference derived from a MODIS MOD09 reconstructed surface reflectance. We compared error distributions at 5%, 20%, and 50% levels of cloudiness and indicated that the cloudiness factor has little impact on the robustness of the atmospheric correction regardless of the band. Finally, the spatial and temporal gradients of accuracy metrics were investigated over northern Eurasia and across different seasons. It was found that for the vast majority of observations, accuracy falls within the −0.010–0.035 range, while precision and uncertainty were below 0.06 for any band. With the successful launch of the most recent Meteor-M #2.3 with a new KMSS-2 instrument onboard, the efficiency and interoperability of the constellation are expected to increase.

1. Introduction

The Meteor-3M series of polar-orbiting satellites are designed to support the Russian Hydrometeorological and Environment Monitoring Service, research institutes, and other organizations with operative Earth-observation data. The satellites are being launched to monitor objects and phenomena in the atmosphere, hydrosphere, and cryosphere for inspecting helio-geophysical conditions in the near-Earth space environment, global climate monitoring, emergency situations tracking, and ecological monitoring [1]. At the moment, there are three operating Meteor-3M satellites, namely, Meteor-M #2, Meteor-M #2.2, and Meteor-M #2.3, and at least one more satellite is scheduled for launch.

The Meteor-M #2 meteorological satellite, which carries the KMSS-M instrument (an acronym for «Multispectral Satellite Imaging System», see Abbreviations), was launched on 8 July 2014 and placed on a sun-synchronous orbit at an altitude of about 830 km [2]. The KMSS-M includes two optical cameras, designated MSU-201 and MSU-202. The cameras were designed for measurements in three spectral bands—green (0.535–0.575 µm), red (0.630–0.680 µm), and near-infrared (0.760–0.900 µm) with a spatial resolution of 60 m and a swath width of about 960 km with both cameras [2] (

Table 1. Key specifications of KMSS series instruments onboard Meteor-M satellites.

Characteristic	KMSS-M (Meteor-M #2)	KMSS-2 (Meteor-M #2.2)
Technique	Pushbroom	Pushbroom
Altitude	832 km	820 km
Spectral bands	Green: 0.535–0.575 μm Red: 0.630–0.680 μm NIR: 0.760–0.900 μm	Green: 0.520–0.590 μm Red: 0.640–0.690 μm NIR: 0.785–0.900 μm
Spatial resolution at nadir	60 m	55 m
Revisiting rate	5 days	5 days
Radiometric resolution	8 bit	10 bit
Total swath	960 km	1020 km

). Meteor-M #2.2 was launched on 5 July 2019 with the KMSS-2 instrument onboard [1]. The most recent Meteor-M #2.3 was launched on 27 June 2023; it carries the KMSS-2 instrument and operates in flight test mode at the moment.

Note that there is a difference in viewing geometry for KMSS-M and KMSS-2 instruments. Specifically, in KMSS-M, three CCD linear sensors of corresponding spectral bands are spaced along the focal plane in the direction of the satellite’s orbital motion, which leads to interband parallax, pronounced on multispectral images of elevated objects such as clouds, mountain peaks, as well as to a time lag between recordings of reflectance in different spectral bands. In the domain of remote vegetation monitoring, along-path multi-angular measurements are attractive because they facilitate and improve biophysical properties retrieval [3,4,5]. Unlike KMSS-M, KMSS-2 features a classical nadir-viewing geometry (Figure 1).

KMSS data are advantageous for operational monitoring, change detection, and vegetation state assessment [6,7]. The Meteor-M\KMSS satellite system is characterized by a demanded combination of spatial resolution and revisit rate, providing that Russia is completely covered with observations every 3–5 days. In tandem, two Meteor-M satellites provide full coverage of Russia on a nearly daily basis. Until recently, much of the potential of KMSS data could not be effectively used due to a lack of atmospheric correction and poor compatibility with standard methods, widely used for optical remote sensing and automated workflows for preprocessing and thematic analysis.

Automated workflows for KMSS data processing are aimed at imagery georeferencing, clouds, and shadows detection, as well as atmospheric correction to provide seasonal and multiyear time series of surface reflectance and vegetation indices over northern Eurasia, including the grain belt of Russia [8]. KMSS Level-2 data prepared with the abovementioned workflow were shown to be qualitatively compatible with similar products of other sensors in this domain [8], indicating their interoperability with standard algorithms and developed workflows for time series analysis, satellite-based land use land cover (LULC) mapping and vegetation parameters retrieval [9,10,11,12,13,14,15,16]. Specifically, atmospherically corrected and gap-free time series of KMSS imagery were used for cropland thematic mapping over the Southern Federal District of Russia and to identify crop parcels using earlier developed routines [17].

Recently, the technology was performing in an operative mode for the first time to support timely and wide-scale remote vegetation monitoring with KMSS data, including continental-scale LULC mapping, crop types mapping, and crop state assessment. However, operatively-derived KMSS datasets were not publicly released yet, and neither the accuracy of produced atmospherically corrected KMSS datasets were analyzed, as it was done for other satellite systems such as TERRA/MODIS, Landsat-8/OLI, and Sentinel-2/MSI [18,19]. To assess the impact of various parameters used in the atmospheric correction, including choice of reference datasets, environmental variables, and workflow constraints, sensitivity analysis is performed [20,21,22].

It is worth noting that though the workflows are developed both for KMSS-M and KMSS-2 data processing, here we focus on operatively derived datasets, applicable for the KMSS-2 instrument only since the Meteor-M #2 mission is likely to be over in 2023 and only historical datasets for the KMSS-M instrument will be produced further on. Thus, the objective of this study was to (1) provide a description of KMSS-M and KMSS-2 data processing aimed at atmospheric correction and (2) perform an accuracy assessment for KMSS-2 atmospheric correction performed over northern Eurasia. In this study, we adopted commonly used APU metrics [18,23,24] to carry out a quantitative and comprehensive analysis of the datasets derived.

2. Materials and Methods

2.1. Satellite Data Acquisition

The Scientific Research Center (SRC) “Planeta” [25,26] acquires, processes, and distributes remote sensing data from various satellite systems, including data derived from the Russian Earth observation constellation: Kanopus-V, Resurs-P, Meteor-M, Electro-L, and Arktika-M series. The SRC Planeta relies on a network of three regional centers (European: Moscow, Obninsk, and Dolgoprudny; Siberian: Novosibirsk; and Far Eastern: Khabarovsk), as well as more than 70 stationary and mobile autonomous platforms in Russia, Antarctica, and maritime downlink stations. This network provides users with full and regular coverage of remote sensing data over Russia and Eurasia. To secure prompt handling of datasets from various centers, a joint data handling system was earlier developed in cooperation with the Space Research Institute of the Russian Academy of Sciences (IKI RAS) [25]. It is designed to provide operational access to distributed data archives, as well as their processing and joint analysis. The joint data-handling system was used for the acquisition of raw Meteor-M-series data with the following processing and storing outputs in the archives of the Center for Collective Use (CCU) [27,28,29] in the Space Research Institute.

Reference MODIS Terra\Aqua (MOD09 collection 6 product) granulated daily scenes data in sinusoidal projection were acquired through dedicated distribution networks (i.e., LP DAAC [30]).

2.2. Methodology

This section describes the methods used both for KMSS and auxiliary satellite data processing to facilitate KMSS atmospheric correction and its validation. Ancillary MODIS data were processed in order to generate reference imagery, i.e., reconstructed gap-free daily surface reflectance in NIR and Red bands. Processing of raw KMSS data included image georeferencing, cloud masking, and atmospheric correction. The final subsection describes the framework built for discrepancy analysis of atmospherically corrected KMSS imagery using APU metrics and MODIS reference datasets over large scales.

2.2.1. MODIS Data Processing

Methodology for georeferencing, cloud and shadow screening, as well as atmospheric correction of KMSS imagery, required the use of coarse-resolution reference imagery produced based on the MODIS MOD09 surface reflectance product in the Space Research Institute of RAS [31]. Time series of MODIS surface reflectance acquisitions were processed using the LOWESS algorithm [32], i.e., weighted sliding-window local regression method for pixel-wise reconstruction of dense seasonal or multiannual gap-free time series of daily surface reflectance or vegetation index. This approach is applicable to a time series of atmospherically corrected acquisitions from various satellite systems with a sufficient revisiting rate as it conveniently bypasses both cloud identification and temporal compositing stages of preprocessing, which may be sensor-specific. The implementation of the abovementioned method ingests a time series of observations with a viewing angle below 45 degrees, accompanied by respective values of their weight factor, i.e., a parameter sensitive to the presence and degree of interfering factors affecting observation quality. For MODIS daily surface reflectance reconstruction, the Normalized Difference Snow Index (NDSI) was used as a weight factor. For this aim, the MOD09GA product was used to discard high viewing angles and to produce NDSI at 500 m resolution. Then, seasonal time series of the target band’s surface reflectance at 250 m resolution were processed with LOWESS using the corresponding time series of NDSI values as weight factors.

In a previous study, it was shown that the mean absolute error (MAE) of the reconstruction was below 5% for the Red and NIR bands [31]. To assess MAE, the dedicated experiment was carried out when cloud- and shadow-free areas were manually selected on a yearlong time series of MODIS scenes over the European part of Russia. Afterward, each of these input scenes was successively excluded before the reconstruction of the seasonal time series, and the excluded image was compared to the correspondent LOWESS-reconstructed image over these selected areas.

It was also indicated that LOWESS-reconstructed MODIS reference imagery features higher geolocation accuracy than the original MOD09 datasets. Namely, using Sentinel-2 (MSI) imagery of decameter resolution as a reference, it was shown that the average geolocation error of LOWESS-reconstructed MODIS was below 20 m [33], making these MODIS images suitable for sub-pixel KMSS imagery georeferencing.

In the last step, MODIS-reconstructed data are projected into an equirectangular (geographical) projection and a KMSS-based grid with spatial resolution, fourfold of KMSS.

2.2.2. KMSS Data Processing

A detailed description of components of KMSS data processing workflow was given in several studies [6,8,31,33]. In this section, we focus on key elements of the workflow and their performance.

The general processing chain, aimed at atmospheric correction of KMSS imagery, consists of three units following each other: georeferencing, cloud and shadow identification, and localized moving window histogram matching (Figure 2).

Initially, the workflow implements rescaling and spatial re-arrangement of the source raw KMSS data from swaths to a one-degree granular grid in geographical projection following the introduction of new naming nomenclature for derived products to facilitate further spatiotemporal analysis, cataloging, and distribution. Within the nomenclature, the gridded product list includes product denomination, acquisition and processing date and time, sensor, and granule ID.

The calibration of raw KMSS data usually takes place on various downlink stations using calibration files distributed earlier by the sensor manufacturer. Ingestion of new calibration files into the workflow is not automatic and thus may pass with delay, while such an update is not always supported on some downlink stations. This factor results in the coexistence of various calibrations, ruling out the interoperability with standard approaches to atmospheric correction, like DOS or FLAASH. KMSS itself lacks onboard calibration hardware, except for the availability of dark current adjustment. This explains some complexity in methodological approaches and preference for spatial and multisource analysis rather than multispectral data for some operations, like georeferencing and cloud identification.

The rescaling provides that input data are linearly recomputed and coerced into a fixed range of values before georeferencing starts. KMSS imagery georeferencing issues originate from several sources: errors in the identification of satellite orientation (and thus in viewing geometry), varying distortions of optical systems due to external factors, and secondary shifts due to parallax bands alignment at the Earth’s surface level with digital elevation model. Accounting for this, the georeferencing step is performed in two tiers: first, evaluation and compensation of systematic offset, and second, identification of smaller displacements for the regular spatial grid of nodes with the following bilinear interpolation between nodes. The secondary objective of this step is to provide a reference for a single-band cloud detection algorithm in the absence of stable radiometric calibration.

Both tiers of georeferencing are based on the identification of the optimal position of the KMSS image or its part relative to the reference image where the value of Pearson’s correlation factor (PCF) between two sources is at maximum. The minor differences between tiers are the extent, search area, spatial step size, and gridding. Specifically, at the first tier, the procedure compensates for the systematic shift of the KMSS image using statistics collected from the target and adjacent granules on a MODIS (fourfold sized) pixel grid, while at the second tier, the procedure operates with only a squared image block around the specific node of the original image on a sliding KMSS pixel grid. When calculating PCF at the second tier, the portions of the KMSS images are sequentially coarsened within a sliding window of 4 by 4 pixels, and thus derived degraded image is compared with the MODIS reference. While searching for the optimal position of the image KMSS node, blocks and grids are moved with a spatial step equal to the original KMSS pixel size. Since the approach employs spatial correlation and is based on the matching of terrain surface patterns between a cloud-free reference image and the original KMSS image, it is applicable when a fraction of a clear surface is not marginal for the KMSS image patch, whether it is a granule or just a node block. At this step, no cloud mask is present, and in order to arrange and qualify the nodes, a gradient analysis of PCF together with node-returned shifts is performed. Gradient analysis starts from the nodes with PCF higher than the threshold—these nodes are considered as cloud-free (qualified)—and iteratively attaches adjacent nodes if both the relative shift and relative PCF decrease value for them are small, and thus attached nodes attain «qualified» status as well. At the end of the georeferencing, qualified nodes provide vectors of compensation shifts, which are interpolated throughout the node grid to produce continuous coverage, while unqualified nodes are considered «mostly clouded» and are used for cloud detection. Such an implemented georeferencing step provides necessary baseline information for single-band cloud identification.

The resultant mean absolute error of geolocation for all qualified nodes processed over the Russian grain belt does not exceed 40 meters for KMSS imagery acquired during the season of 2020, as was shown earlier [33] using Sentinel-2 (MSI) imagery as a reference. It was also shown that geolocation misalignments between source KMSS and reference Sentinel-2 imagery are not uncommon to exceed 3000 meters.

Since KMSS data are considered as not calibrated, there is no use in relying on multiband tests and indices for cloud detection; this methodology focuses on spatial statistics with the assistance of a MODIS-based reference. Moreover, in such cases, spectral bands could only be treated independently, which may affect the accuracy of cloud detection even though a high-quality reference is involved.

For cloud detection, the methodology uses a node status: if a node is qualified, its image block represents mostly a clear surface, while unqualified nodes and their blocks correspond to clouded areas. This information is further used for (i) evaluation of KMSS –MODIS covariance, building a linear regression and per-band recalibration coefficients for cloud-free nodes, and (ii) node-adaptive threshold value determination for cloud detection in a considered band. Eventually, cloud detection is performed based on joint analysis of uncalibrated band reflectance, node threshold value, and the degree of (positive) deviation from the KMSS-MODIS trend line [6,8]. Specifically, the KMSS pixel is masked as a «cloud» if three conditions are met: (1) its node is «unqualified» or adjacent to such, (2) pixel value surpasses the KMSS-MODIS trend line by more than two standard deviations, and (3) pixel value is higher than this node threshold value. Note that the conditions above are met for snowed pixels as well because of the high brightness of snow on remote sensing images in both bands.

Cloud shadow identification is based on Sun illumination geometry reconstructed using geographic coordinates and image acquisition time. The area of potential cloud shadows is determined using illumination geometry as well as the degree of (negative) deviation from the KMSS-MODIS trend line within a potentially shadowed area [6,8]. That is, a pixel is attributed to «shadow» if two conditions are met: (1) it is inside a casted area, defined by cloud mask and illumination geometry, and (2) its value does not reach the KMSS-MODIS trend line by more than one standard deviation.

Finally, atmospheric correction is performed, which is based on a histogram-matching approach, where the KMSS surface reflectance histogram is stretched to match the reference (Figure 3). For better performance and to enhance method flexibility, it utilizes a sliding window of varying size linked to the nodes grid. The sliding window collects localized cloud-free distribution of image pixel band values for associated data sources: cloud-masked KMSS images and the MODIS reference matched by band, date, and territory, which is masked with a KMSS cloud mask. For each focal sliding window centered in the node, after reaching minimum statistical requirements during the size variation, linear transformation coefficients are calculated using the least squares method and applied to the associated KMSS distribution within the window.

The conversion applied to a KMSS pixel value yields a proxy estimate for bottom-of-atmosphere surface reflectance. As a result, each pixel of the KMSS image obtains as many estimates as many focal windows were included. Moreover, each focal window’s estimate is characterized by the agreement level between matched histograms, which is further used as a weight factor for the estimate. The resultant surface reflectance value SR is calculated as the weighted average of n estimates sr_i, with each window’s weight factor w_i. That is,

$$SR=\frac{{\sum }_{i=1}^n{w}_{i}s{r}_i}{{\sum }_{i=1}^n{w}_i},$$

(1)

where the window’s weight factor is inversely proportional to the mean absolute pairwise difference between MODIS and KMSS coerced histograms for each k-th reflectance stratum from the stratum range {$k\in [0\dots N]$} corresponding to reflectance range from 0 to 1:

$$w_i=\frac{1}{{\sum }_{k=1}^N\left|A_k^M-A_k^K\right|},$$

(2)

where A^M and A^K correspond to MODIS and KMSS-matched histogram values (red and blue lines in Figure 3d), respectively. Both KMSS and MODIS histograms are equal in volume, as the statistics are gathered from the same area using a single cloud and shadow mask, so each histogram’s cumulative value may be taken as 1 for simplicity (numerator in Formula (2)). Therefore, both A^M and A^K correspond to each histogram’s bins normalized on histogram cumulated value.

Formula (2) suggests that focal windows with higher discrepancies between two sources yield lower weight and, thus, lower impact on the atmospheric correction result. This enhances the robustness of the approach in case of erroneous or unreliable cloud and shadow identification.

Since reference MODIS imagery is based on the LOWESS-reconstructed time series of a large number of satellite observations with relatively low viewing angles, this atmospheric correction of KMSS images suppresses radiometric effects arising due to differences in observation and illumination angles across a wide KMSS swath width. In addition, this implementation allows the processing of each spectral band independently, which provides additional flexibility.

2.2.3. Accuracy Assessment Framework

The statistical metrics accuracy (A), precision (P), and uncertainty (U) were estimated as below:

$$A=\frac{1}{n}\sum _{i=1}^n{∆\rho }_{\lambda },$$

(3)

$$P=\sqrt{\frac{1}{n-1}\sum _{i=1}^n{({∆\rho }_{\lambda }-A)}^2},$$

(4)

$$U=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left({∆\rho }_{\lambda }\right)}^2},$$

(5)

where n corresponds to the number of samples in a stratum, λ identifies the spectral band (either Red or NIR), and ${∆\rho }_{\lambda }$ equals to the signed difference between KMSS and MODIS reflectance at wavelength λ and is calculated as follows:

$${∆\rho }_{\lambda }={\rho }_{\lambda }^{SRKMSS}-{\rho }_{\lambda }^{SRMODIS}.$$

(6)

To effectively perform accuracy assessment over continental scales, a specific framework was developed (Figure 4). To produce APU metrics in an explicit way, we started with the calculation of the pixel-wise difference between the target KMSS and the associated MODIS reference image. Toward this end, the following steps were made. First, the KMSS image was resampled to the MODIS pixel grid and resolution when KMSS pixels under the MODIS pixel’s footprint were averaged. Second, for each KMSS observation per granule, the MODIS reference pixel was compared with the coarsened KMSS pixel, and thus, the produced residual was attributed to the reference value, forming a two-dimensional density distribution. The distribution was stored in the form of a table and provided the number of cases when a specific residual was found for a given reference surface reflectance stratum. Each table, derived for a specific date and granule, was then stored in a separate file associated with the necessary metadata information.

2.2.4. Samples Density Map

The data sampling map was generated from a KMSS one-degree granule grid (Figure 5). The samples were taken from the cloud-free area of each observation of each granule involved in the analysis. The sample size (dot sizes in Figure 5) is linked to the upper left corner of the respected granule.

To assess the impact of cloudiness on atmospheric correction and thus the sensitivity to cloud identification, we performed the analysis for three different cloudiness levels: 5%, 20%, and 50% according to cloud masks, derived using the cloud masking algorithm. The dot size in Figure 5 depicts the number of observations acquired for a given cloudiness level, whereas the dot color represents the level of cloudiness under consideration. This results in sample size heterogeneity throughout the study area, which is associated with varying mean area cloudiness within the target season as well as with the impact of other hindering factors, like snowcaps in elevated areas, e.g., the northern Far East and Siberian Plateau.

Nevertheless, we believe that the high sampling density and wide extent of the area of analysis ensured the high representativeness of the results.

3. Results

3.1. MODIS Gap-Free Daily Coverage

A time series of gap-free MODIS daily surface reflectance was produced for the snow-free season of 2022 to provide a reference for the KMSS processing workflow. The time series contains reconstructed daily observations, which are not affected by hindering factors, in the Red and NIR MODIS bands. Processing results were in the original MOD09 sinusoidal projection and covered the area of 20 MODIS tiles, corresponding to the grain belt of Russia and its surroundings. The example of daily NDVI coverage produced from reconstructed daily NIR and Red is depicted below (Figure 6).

3.2. KMSS Data Processing

The computational capacities of CCU were used to process all KMSS-2 acquisitions obtained during the season of 2022 over the Russian grain belt. Approximately 2190 total swathes (by both MSU-221 and MSU-222 cameras) acquired over Russia were converted into 150,063 granule dates (which roughly equals 68 dates per granule). For each granule, the referencing, cloud, shadow masking, and atmospheric correction were performed. The processing of a single NIR-band KMSS-2 image and the outputs are depicted on Figure 7.

3.3. APU Database

The APU framework was used to generate the database for accurate assessment of the area of interest. Processing all MODIS-KMSS pairs with regular one-degree sampling over the grain belt of Russia resulted in a file database with a total uncompressed storage volume of more than 4 TB, and more than 260,000 files were produced for both bands and all cloudiness levels. Database entities featured a specific structure of their names, which included granule ID, band, acquisition time, and cloudiness level. This naming provides fast data handling during its aggregation through successive tables merging. The resultant APU metrics provided below were derived from the database while merging the entities passed through the predefined filters on band and cloudiness levels.

3.4. Accuracy Assessment

The accuracy metrics were calculated using all acquisitions obtained over the grain belt during the snow-free part of the 2022 season using the APU database for three very different levels of cloudiness. To illustrate the impact of cloudiness on accuracy, precision, and uncertainty for a given band, corresponding APU metrics were compared on one plot (Figure 8). Quantity distributions, i.e., the number of observations for a given reference surface reflectance stratum for three levels of cloudiness, are depicted as dotted curves with levels of grey.

The accuracy curves for different levels of cloudiness follow the same trends and almost coincide over most dense portions of the distributions. The metrics demonstrate notable peaks around low reference values, which indicate that dark objects are the most challenging cases to account for with a given approach. For the Red band, these objects represent dark, dense vegetation, like coniferous canopies, while in NIR, they most often correspond to open dark soils, like chernozems. This overestimation is explained by the inability of the histogram-matching approach to deal with atmospheric scattering.

One may observe the notable distribution of an accuracy metric when part of the curve lies in the positive domain while the other part crosses into the negative domain. This is another consequence of adopting the histogram matching approach, which performs worse for the edges of the distribution than for its center.

To understand the general impact of cloudiness and, thus, the quality of cloud identification on the accuracy over a continental-scale area extent, we calculated weighted average accuracy metrics using observation distribution as a weight, thus producing an integrated estimate of metrics for three levels of cloudiness (

Table 2. Aggregated statistics for all sampled acquisitions in NIR and RED bands at three levels of cloudiness according to routinely derived cloud masks.

Cloudiness Level	Accuracy MW 1 RED/NIR	Precision MW RED/NIR	Uncertainty MW RED/NIR
5%	0.0070/0.0044	0.0286/0.0339	0.0296/0.0348
20%	0.0065/0.0038	0.0287/0.0370	0.0295/0.0379
50%	0.0059/0.0034	0.0277/0.0384	0.0285/0.0394

¹ MW stands for mean weighted values; bold indicates the lowest values among levels of cloudiness.

). One may see that mean weighted accuracy metric errors do not grow with cloudiness level; instead, the metrics slightly decrease for both bands because extended distributions of higher cloudiness generally spread on strata with lower absolute errors. Meanwhile, mean weighted precision and uncertainty demonstrate multidirectional variations against cloudiness level.

The variations of mean weighted accuracy with the cloudiness level are notably small and never exceed 0.002, while other metrics’ variations do not exceed 0.005 for any spectral band. The comprehensive analysis of the cloud identification accuracy is beyond the scope of this study, but evidently, the cloud factor has a low impact on the robustness of the atmospheric correction due to its inherent mechanism aimed at reducing the influence of disturbing factors.

Given the above, the detailed analysis of derived metrics could be performed no longer considering the cloudiness level parameter, fixing its value. Thus, in order to assess temporal and spatial trends of APU metrics over extended regions of northern Eurasia, we performed corresponding stratifications of the source dataset at a cloudiness level of 50%, hereby involving the analysis of all the data available. Considering temporal stratification, we divided the target year 2022 into three temporal strata—winter–spring (months January till May), summer (June till August), and autumn–winter (September till December) periods in order to account for seasonal variations. In order to consider the spatial heterogeneity of northern Eurasia, we stratified the region into a regular grid of rectangular cells (2° in latitudinal direction and 8° in longitudinal direction) to facilitate visual interpretation of the result over the longitudinally elongated continent. For each non-empty cell, the mean weighted value of APU metrics was produced within the abovementioned temporal strata (Figure 9 for the NIR KMSS-2 band and Figure 10 for the Red KMSS-2 band).

The lowest absolute value and the highest spatial homogeneity for accuracy both for NIR and Red bands correspond to summer months, the period which is also characterized by the most dense sampling throughout northern Eurasia. Precision and uncertainty metrics demonstrate no significant spatial gradient in the summer period either. In particular, accuracy falls within the −0.005–0.006 range for both bands, while precision and uncertainty are below 0.055 for the NIR band and below 0.025 for the Red band in the summer period. The winter–spring period corresponds to both the highest values and spatial heterogeneity of APU metrics.

The spatial variability of APU metrics from the maps above is given in

Table 3. A range between the 5th (left value) and 95th (right value) percentiles of mean weighted APU values, that were observed over northern Eurasia for each of the three seasons. All values are multiplied by 10². The highest and the lowest values for each metric and band throughout the seasons are in bold.

Season	MW 1 Accuracy, 10−2		MW Precision, 10−2		WM Uncertainty, 10−2
	Red	NIR	Red	NIR	Red	NIR
January–May	0.03–3.42	−0.13–1.94	0.82–6.27	1.88–5.60	0.84–6.98	1.90–5.98
June–August	−0.01–0.22	−0.44–0.53	0.70–2.17	2.41–5.17	0.71–2.25	2.45–5.35
September–December	−0.03–1.63	−0.12–0.83	0.78–4.55	1.57–4.29	0.80–4.88	1.58–4.54

¹ MW stands for mean weighted values.

in view of the interpercentile range between the 5th and 95th percentiles. It confirms that the lowest variability corresponds to the summer season, while the highest APU values correspond mostly to the winter–spring period.

4. Discussion

In this study, we analyzed and assessed the effectiveness of the localized histogram matching approach for atmospheric correction of operational KMSS-2 data over continental scales. For the first time, KMSS-2 atmospherically corrected datasets were produced in an operative mode and quantitatively assessed over the grain belt of Russia using MODIS MOD09 reconstructed daily gap-free surface reflectance imagery as a reference.

It was noted that high georeferencing and radiometric characteristics of LOWESS-reconstructed reference imagery secured both adequate atmospheric correction and accuracy assessment, resulting in good performance. It was also shown that atmospheric correction performs well in an operational mode, provided by the advancement of the previously developed methodology. Specifically, it was found that for the vast majority of observations, accuracy falls within the −0.010–0.035 range, while precision and uncertainty were below 0.06 for any band. The accuracy metrics graphs indicate that dark objects are the most challenging to account for because atmospheric scattering is not considered within this approach.

To assess atmospheric correction errors at different cloudiness levels, we compared error distributions and their mean weighted values at 5%, 20%, and 50% levels of cloudiness, and we found that their variations do not exceed 0.002 for an accuracy metric and do not exceed 0.005 for other metrics for both bands; this robustness is attributed to the effect of multiple rated estimators with following weighted averaging, implemented within an atmospheric correction algorithm.

Considering the spatial and temporal heterogeneity of the results, it was found that minimum discrepancy of metrics values is observed in the summer period when accuracy falls within −0.005–0.006 and both precision and uncertainty are below 0.055, while the highest heterogeneity and the maximum APU values are reached in the winter–spring period when accuracy reaches 0.035, and both precision and uncertainty reach 0.07 for the Red band (0.02 and 0.06 for the NIR band, accordingly). The latter may be attributed to the snow factor, as it persists on KMSS-2 imagery for the most part of the winter–spring and autumn–winter periods. Given that no significant spatial gradient of APU values in the summer period is observed, metrics variability is connected with disturbing factors rather than with vegetation cover properties. During changes in snow cover, the MODIS-KMSS histogram adjustment performs worst: unmasked snow under vegetation canopies on KMSS-2 imagery changes contrasts and deforms reflectance histograms, while reconstructed reference MODIS imagery suppresses contaminating factors and portrays snowless cover—these factors hamper accurate superimposing of the histograms. Snow-specific masks are regarded as a way to a solution to this issue. Also, within the abovementioned seasons, low sun zenith angles, high latitudes, and rapid illumination changes within several hours of difference for Terra, Aqua, and Meteor-M satellite flybys increase non-linear histogram discrepancy between source and reference imagery.

This type of atmospheric correction does not account for multi-angular observations because KMSS imagery is matched, though in a spatially adjusted manner, to a result of time series reconstruction, where this factor is suppressed. It is worth noting that the resultant surface reflectance is not equal to surface reflectance derived using atmospheric transfer models and a solution for the bidirectional reflectance distribution function (BRDF). Thereby, the produced KMSS dataset is designed for LULC mapping and LULC change mapping and cannot be used as a dataset for the analysis of surface-reflectance properties.

5. Conclusions

The study presents the methodology and a quantitative assessment of atmospheric correction for KMSS imagery with a focus on northern Eurasia. We believe that the results provided within the study benefit the interoperability of newly produced KMSS datasets and give insights into the performance of the method that is applicable to more recent and higher-resolution imagery.

Future studies include further development of methodology towards better georeferencing involvement of convolution neural networks for accurate cloud identification, as accounting for multi-angular observations in atmospheric correction. Planned processing of multiyear KMSS datasets will provide opportunities for a long-term assessment of the methodology’s performance over multiple years.

The developed processing chain relies on rather sophisticated approaches, which eventually pay off in a good overall performance. Nevertheless, further advancement of these methodologies may be reached in a more straightforward way if new KMSS-class devices onboard scheduled Meteor-M missions and ground segments both will be able to process and provide KMSS data of higher quality. With the successful launch of the most recent Meteor-M #2.3 with a new KMSS-2 instrument onboard on 27 June 2023, the efficiency and interoperability of the constellation are expected to increase.

Data preparation and processing for this study were performed on the capacities of CCU “IKI-Monitoring” [27].