Combining Gaussian Process Regression with Poisson Blending for Seamless Cloud Removal from Optical Remote Sensing Imagery for Cropland Monitoring

: Constructing optical image time series for cropland monitoring requires a cloud removal method that accurately restores cloud regions and eliminates discontinuity around cloud boundaries. This paper describes a two-stage hybrid machine learning-based cloud removal method that combines Gaussian process regression (GPR)-based predictions with image blending for seamless optical image reconstruction. GPR is employed in the ﬁrst stage to generate initial prediction results by quantifying temporal relationships between multi-temporal images. GPR predictive uncertainty is particularly combined with prediction values to utilize uncertainty-weighted predictions as the input for the next stage. In the second stage, Poisson blending is applied to eliminate discontinuity in GPR-based predictions. The beneﬁts of this method are illustrated through cloud removal experiments using Sentinel-2 images with synthetic cloud masks over two cropland sites. The proposed method was able to maintain the structural features and quality of the underlying reﬂectance in cloud regions and outperformed two existing hybrid cloud removal methods for all spectral bands. Furthermore, it demonstrated the best performance in predicting several vegetation indices in cloud regions. These experimental results indicate the beneﬁts of the proposed cloud removal method for reconstructing cloud-contaminated optical imagery.


Introduction
The growing population and climate change present significant challenges to food security [1].The increasing climate variability has a significant impact on food production, emphasizing growing needs for sustainable agricultural management to meet food security needs [2].Consistent cropland monitoring and thematic information extraction are essential for sustainable agricultural management [3].In this context, remote sensing imagery is regarded as an important source of information owing to its ability to provide periodic thematic information in croplands.The representative thematic information derived from remote sensing imagery includes crop type maps [4,5] and crop yield prediction information [6][7][8].Remote sensing-based crop monitoring often requires image time series to extract information on the growth cycles of crops of interest.However, the acquisition of optical remote sensing imagery is greatly affected by atmospheric conditions, making it challenging to collect multi-temporal cloud-free optical images.The presence of clouds and cloud shadows in optical imagery greatly reduces the usability of data for multi-temporal cropland monitoring.When optical imagery contains clouds and cloud shadows, cloudcontaminated regions are first detected and then masked out in the imagery.Representative cloud detection tools include Sen2Cor [9], MACCS-ATCOR Joint Algorithm (MAJA) [10], and AgroShadow [11].Even though cloud regions are accurately detected, these regions are excluded from further analysis.Thus, image reconstruction or gap-filling in cloud by applying Poisson blending.The main contribution of this approach is to provide a novel pipeline for cloud removal that includes a modification procedure for GPR-based initial predictions that is specific to Poisson blending-based discontinuity elimination.As initial prediction results are directly fed into the Poisson blending procedure, uncertainty in GPR-based prediction affects the gradient computation for discontinuity elimination.In this work, uncertainty-weighted predictions are used as input for the gradient computation to alleviate the impact of the prediction uncertainty on Poisson blending.The uncertaintyweighted scheme in particular is proposed because GPR provides the prediction uncertainty, and the gradient information (not the original image values) is considered within Poisson blending.The potential of the proposed approach was demonstrated via cloud removal experiments using Sentinel-2 images of two agricultural sites.

Methods
In this study, cloud removal is undertaken to reconstruct cloud-contaminated remote sensing imagery from the prediction date (T P ) using the cloud-free imagery from the reference date (T R ).The two-stage hybrid cloud removal method proposed in this study combines GPR-based predictions with Poisson blending-based discontinuity elimination (Figure 1).
Agronomy 2023, 13, x FOR PEER REVIEW 3 of 16 by applying Poisson blending.The main contribution of this approach is to provide a novel pipeline for cloud removal that includes a modification procedure for GPR-based initial predictions that is specific to Poisson blending-based discontinuity elimination.As initial prediction results are directly fed into the Poisson blending procedure, uncertainty in GPR-based prediction affects the gradient computation for discontinuity elimination.
In this work, uncertainty-weighted predictions are used as input for the gradient computation to alleviate the impact of the prediction uncertainty on Poisson blending.The uncertainty-weighted scheme in particular is proposed because GPR provides the prediction uncertainty, and the gradient information (not the original image values) is considered within Poisson blending.The potential of the proposed approach was demonstrated via cloud removal experiments using Sentinel-2 images of two agricultural sites.

Methods
In this study, cloud removal is undertaken to reconstruct cloud-contaminated remote sensing imagery from the prediction date (T ) using the cloud-free imagery from the reference date (T ).The two-stage hybrid cloud removal method proposed in this study combines GPRbased predictions with Poisson blending-based discontinuity elimination (Figure 1).

Initial Prediction Using Gaussian Process Regression
GPR first quantifies the relationships between T and T using training samples extracted from non-cloud regions both in reference and prediction images.The reflectance values in the cloud masks of the target imagery at T for each spectral band are then predicted using the quantified relationships and reflectance values from the cloud masks in the reference imagery.The initial prediction result based on GPR at T is an image in which cloud masks are replaced with GPR-based prediction results, while non-cloud regions are retained.
GPR learns a probability distribution over functions via a stochastic Gaussian process (GP) [39,40].In GPR, given a set of training samples,  = {( ,  ),  = 1, ⋯ , }, the output  is formulated as the sum of some unknown latent function  at an input  , and the Gaussian noise with zero mean and variance  :

Initial Prediction Using Gaussian Process Regression
GPR first quantifies the relationships between T R and T P using training samples extracted from non-cloud regions both in reference and prediction images.The reflectance values in the cloud masks of the target imagery at T P for each spectral band are then predicted using the quantified relationships and reflectance values from the cloud masks in the reference imagery.The initial prediction result based on GPR at T P is an image in which cloud masks are replaced with GPR-based prediction results, while non-cloud regions are retained.
GPR learns a probability distribution over functions via a stochastic Gaussian process (GP) [39,40].In GPR, given a set of training samples, D = {(x n , y n ), n = 1, • • • , N}, the output y n is formulated as the sum of some unknown latent function g at an input x n , and the Gaussian noise with zero mean and variance σ 2 n : Unlike conventional parametric regression models, the latent function g is regarded as a random variable following a particular distribution.In GPR, g(x) is assumed to be distributed as a GP, which is defined as a collection of random variables, any finite number of which have a joint Gaussian distribution [39] (p.13).A GP is completely specified by its mean function µ(x) and covariance function k(x, x ): By combining the zero mean GP prior and the Gaussian likelihood computed from training samples, the posterior distribution over the unknown output for the new observation x * can be analytically computed within a Bayesian framework.After computing , the predictive mean (µ(x * )) and the predictive variance (σ 2 (x * )) at the new observation x * are finally obtained, as follows: where The attractive advantage of GPR over other regression models is its ability to provide prediction uncertainty estimates (i.e., the predictive variance in Equation ( 4)) together with prediction values.
The covariance function in Equation ( 2) is specified by the kernel function measuring the similarity between inputs of a function.The commonly used kernel function is the square exponential kernel, also called the radial basis function (RBF) kernel: where σ 2 g is the signal variance or the vertical length scale, l is the horizontal length scale, and σ 2 n is the noise variance.δ xx is 1 if x = x and is 0 otherwise.The optimal values of the three hyperparameters are usually determined through marginal likelihood maximization.
In this study, an uncertainty-weighted scheme based on the following two strategies is presented to reduce the impact of the prediction uncertainty on discontinuity elimination.
(1) The first strategy is inverse uncertainty weighting.Larger weights are assigned predictions with smaller prediction variances, while smaller weights are given to predictions with larger prediction variances.(2) The second strategy is mean bias correction.When inverse uncertainty weights are normalized, the resulting weighted predictions will likely have much smaller values than the original GPR predictions.Consequently, gradient computation in Poisson blending tends to return smaller gradient values, yielding smoothed blending results.
To avoid smoothing effects, the ratio of the mean values from the initial and weighted predictions was empirically considered another weighting factor to preserve the first momentum of the predictions.By considering the empirical mean ratio-based correction term as another weighting factor, the final prediction has the same mean value as the initial predictions but different variations.
Based on the above two strategies, for initial predictions based on GPR {(x k , y * (x k )), k = 1, • • • , K} within a specific cloud mask, the final weighted predictions (y * ω ) for a specific cloud region are formulated as follows: where σ is the predictive standard deviation of GPR.y 0 = 1/K∑ K k=1 y * (x k ) and y 1 = 1/K∑ K k=1 σ(x k ) −p y * (x k ) are the mean values of initial predictions and inverse uncertainty weighted predictions, respectively.These two values are computed for each cloud region.
Within Equation (6), a weighting power value (p) that controls how proportional the weights are to the inverse of the uncertainty needs to be determined.The larger weight power is likely to generate weighted values that are quite different from the original values.Based on our preliminary tests using different power values, a power value of 0.5 was empirically selected in this work.Thus, the weight assigned to each prediction is proportional to the inverse of the square root of the standard deviation.Hereafter, the final weighted predictions in Equation ( 6) are referred to as the GPR prediction and are used as input for Poisson blending.

Discontinuity Elimination Using Poisson Blending
Poisson blending is a guided interpolation framework [41] and has been applied to the seamless cloning of different images.The key principle of Poisson blending is to seamlessly compose some parts of a target image with a source image through image editing in the gradient domain, not in the original value domain.
Within the process of cloud removal, the GPR prediction image and the cloudcontaminated image correspond to the source and target images, respectively.Given the cloud mask in the target image (Ω) and its boundary (∂Ω), let f Ω be an unknown image reflectance function, defined over the interior of Ω (i.e., target reflectance values within the cloud mask) and let f NC be a known image reflectance function defined outside of Ω in the target image (i.e., reflectance values in non-cloud regions).Poisson blending aims to find the unknown function f Ω , which reconstructs reflectance values within Ω without any discontinuity around ∂Ω (Figure 2).To this end, image editing proceeds such that (1) the source image has the same reflectance value on ∂Ω as the target image and (2) the gradient for each pixel over the interior Ω of the source image should preserve the original gradient of the source image.By imposing these constraints, the reconstructed cloud region features a variation in reflectance corresponding to the GPR prediction results and is also consistent with the reflectance of the non-cloud region in the target image around the cloud boundary.where  is the predictive standard deviation of GPR. = 1/ ∑  * ( ) and  = 1/ ∑ ( )  * ( ) are the mean values of initial predictions and inverse uncertainty weighted predictions, respectively.These two values are computed for each cloud region.Within Equation (6), a weighting power value () that controls how proportional the weights are to the inverse of the uncertainty needs to be determined.The larger weight power is likely to generate weighted values that are quite different from the original values.Based on our preliminary tests using different power values, a power value of 0.5 was empirically selected in this work.Thus, the weight assigned to each prediction is proportional to the inverse of the square root of the standard deviation.Hereafter, the final weighted predictions in Equation ( 6) are referred to as the GPR prediction and are used as input for Poisson blending.

Discontinuity Elimination using Poisson Blending
Poisson blending is a guided interpolation framework [41] and has been applied to the seamless cloning of different images.The key principle of Poisson blending is to seamlessly compose some parts of a target image with a source image through image editing in the gradient domain, not in the original value domain.
Within the process of cloud removal, the GPR prediction image and the cloud-contaminated image correspond to the source and target images, respectively.Given the cloud mask in the target image (Ω) and its boundary (∂Ω), let  be an unknown image reflectance function, defined over the interior of Ω (i.e., target reflectance values within the cloud mask) and let  be a known image reflectance function defined outside of Ω in the target image (i.e., reflectance values in non-cloud regions).Poisson blending aims to find the unknown function  , which reconstructs reflectance values within Ω without any discontinuity around ∂Ω (Figure 2).To this end, image editing proceeds such that (1) the source image has the same reflectance value on ∂Ω as the target image and (2) the gradient for each pixel over the interior Ω of the source image should preserve the original gradient of the source image.By imposing these constraints, the reconstructed cloud region features a variation in reflectance corresponding to the GPR prediction results and is also consistent with the reflectance of the non-cloud region in the target image around the cloud boundary.Under the guidance of vector field v defined over the interior of Ω in the source image, Poisson blending is formulated using the following the minimization problem [41]: where ∇ is the gradient operator.
Solving Equation ( 7) is equivalent to finding the solution of the following Poisson equation with Dirichlet boundary conditions [41]: Under the guidance of vector field v defined over the interior of Ω in the source image, Poisson blending is formulated using the following the minimization problem [41]: where ∇ is the gradient operator.Solving Equation ( 7) is equivalent to finding the solution of the following Poisson equation with Dirichlet boundary conditions [41]: where ∆ and div are the Laplacian operator and the divergence of a vector field, respectively.
Agronomy 2023, 13, 2789 6 of 16 A set of four connected neighbors for the pixel in the target image is usually considered to approximate the differential operators in the discrete image domain.Based on this finite difference approach, a numerical solution of Equation ( 8) is then obtained by solving a sparse linear system [42].

Study Area and Data
Cloud removal experiments were conducted on images of two agricultural sites in Korea, Gimje (Site 1) and Hapcheon (Site 2), as shown in Figure 3.The two sites are major rice and garlic/onion production regions in Korea, respectively.Thus, cloud-free image time series are required for crop monitoring.The transplanting and harvesting periods for rice grown at Site 1 are in May and October, respectively.Garlic and onions at Site 2 have the highest vegetative vitality from late April to early May.Harvesting begins in late May.The total areas of the two sites are 2500 ha and 676 ha, respectively.where ∆ and div are the Laplacian operator and the divergence of a vector field, respectively.
A set of four connected neighbors for the pixel in the target image is usually considered to approximate the differential operators in the discrete image domain.Based on this finite difference approach, a numerical solution of Equation ( 8) is then obtained by solving a sparse linear system [42].

Study Area and Data
Cloud removal experiments were conducted on images of two agricultural sites in Korea, Gimje (Site 1) and Hapcheon (Site 2), as shown in Figure 3.The two sites are major rice and garlic/onion production regions in Korea, respectively.Thus, cloud-free image time series are required for crop monitoring.The transplanting and harvesting periods for rice grown at Site 1 are in May and October, respectively.Garlic and onions at Site 2 have the highest vegetative vitality from late April to early May.Harvesting begins in late May.The total areas of the two sites are 2500 ha and 676 ha, respectively.In this study, Sentinel-2 images were used as input for cloud removal experiments (Figure 3).The Sentinel-2 images have often been utilized for crop monitoring because the combined constellation of Sentinel-2A and -2B satellites provides images every five days.Out of the twelve spectral bands of the level-2A bottom-of-atmosphere (BOA) products [43], the reflectance values of four spectral bands, including green, red, red-edge, and near infrared (NIR) bands, were considered for cloud removal experiments because they have frequently been used to calculate several vegetation indices for crop monitoring (Table 1).The red-edge 1 band with a spatial resolution of 20 m was resampled to 10 m using bilinear resampling to match the spatial resolution of all the spectral bands to 10 m.Among the cloud-free Sentinel-2 images acquired in 2021 for each site, 20 August and 14 April were experimentally selected as the prediction dates for Site 1 and Site 2, respectively, because vegetative vitality was high at both sites on those dates.Based on preliminary tests, reference dates showing a high correlation with the prediction dates were selected for the two sites (26 July and 10 March, respectively).In this study, Sentinel-2 images were used as input for cloud removal experiments (Figure 3).The Sentinel-2 images have often been utilized for crop monitoring because the combined constellation of Sentinel-2A and -2B satellites provides images every five days.Out of the twelve spectral bands of the level-2A bottom-of-atmosphere (BOA) products [43], the reflectance values of four spectral bands, including green, red, red-edge, and near infrared (NIR) bands, were considered for cloud removal experiments because they have frequently been used to calculate several vegetation indices for crop monitoring (Table 1).The red-edge 1 band with a spatial resolution of 20 m was resampled to 10 m using bilinear resampling to match the spatial resolution of all the spectral bands to 10 m.Among the cloud-free Sentinel-2 images acquired in 2021 for each site, 20 August and 14 April were experimentally selected as the prediction dates for Site 1 and Site 2, respectively, because vegetative vitality was high at both sites on those dates.Based on preliminary tests, reference dates showing a high correlation with the prediction dates were selected for the two sites (26 July and 10 March, respectively).

Experimental Design
The synthetic cloud masks shown in Figure 3 were first generated in cloud-free Sentinel-2 imagery at T P to evaluate quantitative prediction performance within the cloud removal experiment.The cloud masks in this study included a large amount of clouds composed of various types of land cover.The cloud type was assumed to be thick clouds that block signals from the land surface, and the cloud masks include both thick clouds and shadows.When clouds are widely distributed over the study area, they may be located at the edge of the image.In this case, f NC is not available when interpolating reflectance values around edge pixels.To solve this limitation, the GPR prediction values of the edge pixels were assumed to be f NC when clouds are located at the edge of the image.Based on this assumption, GPR prediction values were first assigned to the edge of the cloud masks that actually belonged to clouds.Poisson blending was then applied by regarding GPR prediction values as non-cloud values at T P .In considering this case, some cloud masks were particularly placed at the edges of the image, as shown in Figure 3. Finally, the cloud masks, consisting of 1 for the cloud region and 0 for the non-cloud region, were prepared and utilized as input for the cloud removal experiment.The actual reflectance values for each spectral band of the cloud region were used to compute quantitative evaluation metrics.The cloud masks occupy approximately 26% of the study area at both sites.
Based on our previous studies [27,38] and computational efficiency, 1% of the noncloud pixels in the study area, 2500 and 676 for Site 1 and Site 2, respectively, were randomly extracted and then used as training samples for GPR model training.
The prediction performance of the proposed cloud removal method was compared with that of two existing cloud removal methods, including modified NSPI (MNSPI) [29] and geostatistical NSPI (GNSPI) [44].The two methods were selected because they are hybrid methods, like the proposed method, and their source code is publicly available [45].MNSPI reconstructs persistent missing regions through the weighted combination of similar pixels based on spectral and spatial distances from auxiliary spatial information [29].In GNSPI, missing information is predicted by combining regression-based temporal information with kriging of residuals [44].GNSPI utilizes multi-temporal images to extract spectrally similar neighboring pixels.For a fair comparison with MNSPI and the proposed method using a single reference image, this study utilized a single reference image to implement GNSPI.The common parameters for MNSPI and GNSPI are the number of land cover types and the size of the moving window for extracting neighboring pixels, which were set to 7 and 5, respectively.In addition, the GPR prediction was compared with the final prediction to investigate the effect of Poisson blending-based discontinuity elimination.
All cloud removal methods were applied to each of the four spectral bands of the Sentinel-2 image, and prediction results for each spectral band were compared quantitatively and qualitatively.The four spectral bands are often used to calculated the vegetation index, which is the essential information source for crop monitoring.Thus, to highlight the importance of cloud removal for crop monitoring, the reflectance values of the four spectral bands were utilized to calculate the vegetation index.The predictive performance of different cloud removal methods was then evaluated by comparing the accuracy of the calculated vegetation index.This study considered three vegetation indices that can be calculated from the four spectral bands: (1) normalized difference vegetation index (NDVI), (2) normalized difference red-edge (NDRE), and (3) normalized difference water index (NDWI): where ρ denotes the reflectance of a specific spectral band.Two evaluation metrics, the relative root mean square error (rRMSE) and structural similarity index measure (SSIM), were used to measure predictive performance.Since each spectral band has different ranges of spectral reflectance, rRMSE was considered for a relative comparison.SSIMs measure spatial similarity by comparing structural information between actual and predicted reflectance values in cloud regions [46].Smaller rRMSE values indicate higher prediction accuracy.On the other hand, the closer the SSIM value is to one, the better the structural similarity.
Given actual reflectance values (y(x m )) and predicted values ( ŷ(x m )) for a specific spectral band in cloud regions consisting of M pixels (x m , m = 1, • • • , M), rRMSE and SSIM are calculated as follows: where µ y and σ 2 y are the mean and variance values for the actual reflectance values, respectively.µ ŷ are σ 2 ŷ are the mean and variance values for the predicted reflectance values, respectively.COV denotes the covariance between actual and predicted reflectance values.C 1 and C 2 are two constraint constants.
The Scikit-learn library [47] and the Python code of Poisson image editing [48] were utilized to implement GPR and Poisson blending, respectively.All data processing, including weighted predictions and evaluation metrics computation, was implemented using Python programming.

Reflectance Prediction
Table 2 lists the quantitative evaluation results for each spectral band.The proposed method achieved the best predictive performance at both sites.MNSPI showed a worse predictive performance for all spectral bands at both sites, except GNSPI which was worse for the red band at Site 1. Significant improvements in rRMSE compared with the worst method were found in the red-edge and NIR bands at both sites.The improvements in the rRMSE of the proposed method over that of MNSPI were 19.27% and 13.34%, respectively, for the red-edge and NIR bands at Site 1.The proposed method also increased the rRMSE by 11.19% and 12.54% for the red-edge and NIR bands, respectively, at Site 2, compared to MNSPI.The superiority of the proposed method was also shown in the SSIM, except for the green band at Site 2. The maximum relative improvement in SSIM was obtained for the red-edge and NIR bands.These quantitative comparison results indicate that the proposed method can effectively capture the structural characteristics and overall quality of the actual reflectance values in cloud regions.The predictive performance of GPR was slightly worse than that of GNSPI, except for the green band at Site 1 and the NIR band at Site 2. However, GPR outperformed MNSPI for all spectral bands.When comparing the GPR prediction with the final prediction of the proposed method, Poisson blending increased the rRMSE and SSIM for all spectral bands at both sites.Poisson blending interpolates from the actual value of the non-cloud regions to the inside cloud regions to maintain the spectral pattern.Thus, the proposed method yielded a decrease in rRMSE and an increase in SSIM.Since the GPR predictions achieved better predictive performance in a relative sense, the final prediction obtained after applying Poisson blending demonstrated the best prediction accuracy, confirming the potential of the proposed hybrid method for use in cloud removal.
Figure 4 shows the zoomed subareas around the cloud masks.Compared to the GPR predictions, the final prediction results clearly showed natural variations in reflectance around the cloud outline upon visual inspection, demonstrating the necessity of discontinuity elimination via Poisson blending after generating the GPR prediction.
of the actual reflectance values in cloud regions.The predictive performance of GPR was slightly worse than that of GNSPI, except for the green band at Site 1 and the NIR band at Site 2. However, GPR outperformed MNSPI for all spectral bands.When comparing the GPR prediction with the final prediction of the proposed method, Poisson blending increased the rRMSE and SSIM for all spectral bands at both sites.Poisson blending interpolates from the actual value of the non-cloud regions to the inside cloud regions to maintain the spectral pattern.Thus, the proposed method yielded a decrease in rRMSE and an increase in SSIM.Since the GPR predictions achieved better predictive performance in a relative sense, the final prediction obtained after applying Poisson blending demonstrated the best prediction accuracy, confirming the potential of the proposed hybrid method for use in cloud removal.
Figure 4 shows the zoomed subareas around the cloud masks.Compared to the GPR predictions, the final prediction results clearly showed natural variations in reflectance around the cloud outline upon visual inspection, demonstrating the necessity of discontinuity elimination via Poisson blending after generating the GPR prediction.The effect of seam removal via Poisson blending is clearly revealed upon comparison with MNSPI and GNSPI at both sites (Figures 5 and 6).The MNSPI prediction shows noticeable differences in spectral patterns between cloud and non-cloud regions around cloud outlines at Site 1 (Figure 5).Due to the low reflectance in the NIR band in most cloud regions, the image appears relatively dark, and there are also outliers at some cloud The effect of seam removal via Poisson blending is clearly revealed upon comparison with MNSPI and GNSPI at both sites (Figures 5 and 6).The MNSPI prediction shows noticeable differences in spectral patterns between cloud and non-cloud regions around cloud outlines at Site 1 (Figure 5).Due to the low reflectance in the NIR band in most cloud regions, the image appears relatively dark, and there are also outliers at some cloud boundaries.Discontinuities near cloud outlines are somewhat less pronounced in the GNSPI prediction than in the MNSPI prediction, but spectral distortions still exist inside the cloud region.In addition, the spectral reflectance of some pixels near the cloud boundary is low, resulting in irregular spectral patterns within the crop parcel.In contrast, the proposed method restored the reflectance of the cloud region so that it was similar to the actual reflectance.Spectral discontinuity around the cloud outline is alleviated, and the reflectance inside the crop parcel is distributed evenly and continuously.
boundaries.Discontinuities near cloud outlines are somewhat less pronounced in the GNSPI prediction than in the MNSPI prediction, but spectral distortions still exist inside the cloud region.In addition, the spectral reflectance of some pixels near the cloud boundary is low, resulting in irregular spectral patterns within the crop parcel.In contrast, the proposed method restored the reflectance of the cloud region so that it was similar to the actual reflectance.Spectral discontinuity around the cloud outline is alleviated, and the reflectance inside the crop parcel is distributed evenly and continuously.Similar to Site 1, the reconstruction results of the proposed method were also visually better than those of the two existing methods at Site 2 (Figure 6).As shown in subarea A, MNSPI prediction showed, in relative terms, the most prominent discontinuity, and GNSPI also showed blurred boundaries in some crop parcels.The two existing methods failed to capture the structural characteristics of the relatively smaller crop parcels than at Site 1. boundaries.Discontinuities near cloud outlines are somewhat less pronounced in the GNSPI prediction than in the MNSPI prediction, but spectral distortions still exist inside the cloud region.In addition, the spectral reflectance of some pixels near the cloud boundary is low, resulting in irregular spectral patterns within the crop parcel.In contrast, the proposed method restored the reflectance of the cloud region so that it was similar to the actual reflectance.Spectral discontinuity around the cloud outline is alleviated, and the reflectance inside the crop parcel is distributed evenly and continuously.Similar to Site 1, the reconstruction results of the proposed method were also visually better than those of the two existing methods at Site 2 (Figure 6).As shown in subarea A, MNSPI prediction showed, in relative terms, the most prominent discontinuity, and GNSPI also showed blurred boundaries in some crop parcels.The two existing methods failed to capture the structural characteristics of the relatively smaller crop parcels than at Site 1. Similar to Site 1, the reconstruction results of the proposed method were also visually better than those of the two existing methods at Site 2 (Figure 6).As shown in subarea A, MNSPI prediction showed, in relative terms, the most prominent discontinuity, and GNSPI also showed blurred boundaries in some crop parcels.The two existing methods failed to capture the structural characteristics of the relatively smaller crop parcels than at Site 1.

Vegetation Indices Prediction
The quantitative evaluation results of the vegetation index predictions are summarized in Table 3.The mean of the actual NDWI values at both sites was negative, so the absolute mean value was used to calculate the rRMSE of the NDWI predictions.As expected, the proposed method yielded the best predictive performance for the three vegetation indices with the lowest rRMSE and the largest SSIM.Like the reflectance predictions, MNSPI showed the worst prediction accuracy and structural similarity at both sites, except for the NDVI prediction of GNSPI at Site 1.The worst rRMSE of GNSPI in NDVI prediction at Site 1 is attributed to its worse accuracy in red band prediction.The predictive performance of the NDWI prediction regarding GPR was better than that of MNSPI and GNSPI at both sites.The maximum increase in the rRMSE of the proposed method, compared with the worst predictor, was observed in its NDRE prediction (13.25% at Site 1 and 10.69% at Site 2) followed by the NDWI prediction.The reflectance values predicted in cloud regions are directly fed into the calculation of the vegetation index.Thus, the quality of the computed vegetation index is subject to errors in the reflectance predictions.As two spectral bands are utilized to calculate the vegetation index, the errors of each spectral band together affect predictive performance.Although the proposed method yielded the smallest rRMSE in NIR band prediction, the relatively large errors in red band prediction resulted in a smaller improvement in the rRMSE of the NDVI prediction.Meanwhile, the smallest rRMSE in the prediction of the red-edge and NIR bands of the proposed method led to the most significant improvement in the rRMSE of the NDRE prediction.For SSIM, the proposed method showed the best structural similarity for all three vegetation indices at both sites.In particular, the GPR prediction achieved the second-largest structural similarity.

Contribution of the Study
The superiority of our method in this study was attributed to both GPR-based predictions with high quality and seam removal via Poisson blending.It should be noted that the purpose of Poisson blending employed in the second stage is not to improve prediction accuracy.Instead, Poisson blending is used to remove seams and enhance naturalness in reflectance around cloud outlines.Thus, a substantial improvement in prediction accuracy cannot always be expected when applying Poisson blending alone.Therefore, initial predictions with reasonably high accuracy should be used as input for Poisson blending.Compared to random forest regression applied in previous studies [26,49], GPR showed a better prediction performance in cloud removal [38].In addition, the potential of GPR was demonstrated in many existing gap-filling studies [35,36].Based on previous studies, GPR was selected in this study to generate initial predictions for cloud removal.Consequently, although the temporal relationship between the reference and prediction dates was purely modeled and used for prediction in GPR, the performance of the GPR prediction was compatible with or better than that of MNSPI and GNSPI, which take a spatio-temporal hybrid approach.One of the advantages of GPR over other machine learning-based regression models is the availability of the predictive uncertainty information (i.e., predictive variance) along with predictions [39].To date, however, few studies have applied predictive uncertainty information for interpretation or value-added product generation.In this study, predictive standard deviation was combined with predictive value to give less weight to uncertain predictions.Since the gradient, rather than the reflectance value itself, is modeled within the process of Poisson blending, weighted prediction based on two strategies was utilized as the input for Poisson blending.
In their application of Poisson blending to cloud removal, existing hybrid approaches first blend two images acquired on different dates by substituting the non-cloud pixels in the reference image for the cloud region in the prediction image [30,31].Taking such a cut-and-paste approach, Poisson blending may fail to alleviate the discontinuity around cloud outlines when significant changes occur between reference and prediction dates.To mitigate this limitation, the proposed method replaced cloud regions with GPR predictions based on quantitative relationships between reference and prediction dates and predictive uncertainty.By enhancing naturalness through Poisson blending, prediction results could exhibit more significant structural coherence and maintain the overall spectral patterns in the cloud regions, as confirmed by the smallest rRMSE and the highest SSIM in the comparisons at the two cropland sites.The benefits of our method can be further verified through a comparison of the temporal differences in NDVI between reference and prediction dates.The two sites exhibited different temporal relationships between the reference and prediction dates.At Site 1, the average NDVI values for the reference and prediction dates were similar, with values of 0.72 and 0.79, respectively.In contrast, at Site 2, the average NDVI values for the reference and prediction dates were significantly different (0.27 vs. 0.51), indicating a substantial variation in vegetation vitality between the two dates.The predictive performance of the proposed method at Site 2 outperformed that at Site 1, as shown in Table 2.This result confirms the superiority of our method over conventional cut-and-paste approaches, particularly when there are spectral differences between the two dates.

Future Research Directions
While our method has demonstrated superior predictive performance, additional refinements and potential applications should be included in future work.
In this study, the relationship between the reflectance values of the reference and prediction dates in GPR was modeled independently for each spectral band.The spectral bands considered in this study have a reasonable correlation.Thus, it is worthwhile to evaluate the potential of multi-output Gaussian process regression, which can capture quantitative relationships between correlated multiple inputs [50].
Predictive uncertainty in GPR is epistemic uncertainty that indicates how confident the model is with respect to its predictions [51].Thus, predictive uncertainty may not always be related directly to prediction errors.As shown in Equation ( 4), predictive variance is dependent on covariance but independent of data values, like kriging variance in geostatistics [52].As the covariance depends solely on the location of data points, the predictive variance of any data sets separated by the same distance is the same, regardless of the data values.Prediction errors usually differ in the prediction of actual high and low values.Thus, the predictive uncertainty in GPR provides only indirect information on its measure of prediction accuracy, which is the limitation of the uncertainty-weighted approach in this study.An alternative approach is to quantify the relationship between the error and the predictive uncertainty in the non-cloud region and then convert the GPR variance into error-related values since errors can be calculated in the non-cloud region.
Recently, several studies have applied cloud removal to high spatial resolution satellite images, such as PlanetScope and WordView images, using deep learning and objectbased information [53,54].Our cloud removal method can be applied to cloud removal from high spatial resolution imagery without the need for further modification since GPR-based regression modeling and Poisson blending are independent of the spatial resolution of the input imagery.However, it is worth noting that high spatial resolution imagery often exhibits more spectral variations and irregular shapes in surface objects compared to low spatial resolution imagery.These differences in spatial structures and land-cover types between cloud and non-cloud regions may affect the prediction results of our method.Therefore, future work should include extensive experiments using high spatial resolution imagery.
The dense image time series generated via cloud removal can act as essential information sources for environmental monitoring, as well as cropland monitoring.For example, the characterization and classification of forest types often requires a series of multi-spectral remote sensing images to capture phenological changes in the forest, like crop monitoring and classification.However, the occurrence of cloud-contaminated images may degrade predictive performance [55].Thus, optical images reconstructed by cloud removal can be effectively employed for the multi-temporal analysis of forests.
With regard to image time series construction, cloud removal can be combined with other approaches.For example, cloud removal may be employed as a preprocessing step in spatio-temporal image fusion for generating fine spatial and temporal resolution imagery [56].Spatio-temporal image fusion requires both fine temporal resolution but coarse spatial resolution imagery and coarse temporal resolution but fine spatial resolution imagery.Images contaminated by clouds cannot be used as input for spatio-temporal image fusion.Recently, Tange et al. [57] utilized gap-filled images as input to generate all-sky MODIS land surface temperature.Similarly, the cloud removal method proposed in this work can be applied to reconstruct cloud images and combined with a spatio-temporal image fusion method specifically designed for small-scale cropland monitoring [58].The benefits of our cloud removal method for the spatio-temporal fusion of vegetation indices will be evaluated in future work.

Conclusions
The presence of cloud contamination in optical remote sensing imagery poses a significant challenge to cropland monitoring using remote sensing image time series.To address limitations in the usability of cloud-free images, this study presented a two-stage hybrid cloud removal method.Our cloud removal method is a hybrid approach leveraging the advantages of temporal (GPR) and spatial (Poisson blending) methods.The method begins with the generation of regression-based predictions, followed by discontinuity elimination through Poisson blending.Since the quality of input imagery greatly affects the effectiveness of Poisson blending, GPR was employed as the primary regression model due to its superior predictive capability in cloud removal and the availability of predictive uncertainty information.Cloud removal experiments using Sentinel-2 images at two cropland sites demonstrated the effectiveness of our method.The results from our method exhibited reduced spectral discontinuities around cloud outlines and fewer spectral distortions within cloud regions.In comparative evaluations with two existing hybrid methods, our method outperformed them both in terms of rRMSE and SSIM.Moreover, our method was also able to reconstruct vegetation index values within cloud regions, achieving the highest accuracy and structural similarity.These findings confirm that our method can be effectively applied to the construction of cloud-free optical image time series for crop monitoring and thematic mapping.

Figure 1 .
Figure 1.Schematic diagram of the processing flow of the two-stage cloud removal method proposed in this study (T : reference date; T : prediction date).

Figure 1 .
Figure 1.Schematic diagram of the processing flow of the two-stage cloud removal method proposed in this study (T R : reference date; T P : prediction date).

my 2023 ,
13, x FOR PEER REVIEW 5 of 16

Figure 2 .
Figure 2. Poisson blending procedures for seam elimination around the cloud boundary.Ω and ∂Ω represent the interior and boundary of the cloud mask, respectively.Arrows within the cloud mask in the initial prediction result denote the gradient field.

Figure 2 .
Figure 2. Poisson blending procedures for seam elimination around the cloud boundary.Ω and ∂Ω represent the interior and boundary of the cloud mask, respectively.Arrows within the cloud mask in the initial prediction result denote the gradient field.

Figure 3 .
Figure 3. False color composites of Sentinel-2 images in two sites (near infrared-red-green bands as R-G-B): (a) Sentinel-2 imagery on 20 August 2021 at the Gimje site (Site 1); (b) Sentinel-2 imagery on 14 April 2021 at the Hapcheon site (Site 2).White polygons shown in the Sentinel-2 imagery represent synthetic cloud masks.

Figure 3 .
Figure 3. False color composites of Sentinel-2 images in two sites (near infrared-red-green bands as R-G-B): (a) Sentinel-2 imagery on 20 August 2021 at the Gimje site (Site 1); (b) Sentinel-2 imagery on 14 April 2021 at the Hapcheon site (Site 2).White polygons shown in the Sentinel-2 imagery represent synthetic cloud masks.

Figure 4 .
Figure 4. Visual comparisons between GPR and final predictions in subareas at the two study sites.

Figure 4 .
Figure 4. Visual comparisons between GPR and final predictions in subareas at the two study sites.

Figure 5 .
Figure 5. Visual comparison of three prediction results with a false color composite of the actual imagery (NIR-red-green bands as R-G-B) in three subareas at Site 1.The white line marked in the image represents the synthetic cloud mask, the yellow boxes represent the three subareas, and the black dots indicate very low reflectance values.

Figure 6 .
Figure 6.Visual comparison of three prediction results with a false color composite of the actual imagery (NIR-red-green bands as R-G-B) in three subareas at Site 2. The white line marked in the image represents the synthetic cloud mask, the yellow boxes represent the three subareas, and the black dots indicate very low reflectance values.

Figure 5 .
Figure 5. Visual comparison of three prediction results with a false color composite of the actual imagery (NIR-red-green bands as R-G-B) in three subareas at Site 1.The white line marked in the image represents the synthetic cloud mask, the yellow boxes represent the three subareas, and the black dots indicate very low reflectance values.

Figure 5 .
Figure 5. Visual comparison of three prediction results with a false color composite of the actual imagery (NIR-red-green bands as R-G-B) in three subareas at Site 1.The white line marked in the image represents the synthetic cloud mask, the yellow boxes represent the three subareas, and the black dots indicate very low reflectance values.

Figure 6 .
Figure 6.Visual comparison of three prediction results with a false color composite of the actual imagery (NIR-red-green bands as R-G-B) in three subareas at Site 2. The white line marked in the image represents the synthetic cloud mask, the yellow boxes represent the three subareas, and the black dots indicate very low reflectance values.

Figure 6 .
Figure 6.Visual comparison of three prediction results with a false color composite of the actual imagery (NIR-red-green bands as R-G-B) in three subareas at Site 2. The white line marked in the image represents the synthetic cloud mask, the yellow boxes represent the three subareas, and the black dots indicate very low reflectance values.

Table 1 .
Summary of Sentinel-2 images used for the experiments (NIR: near infrared).

Table 2 .
Band-wise accuracy statistics of different cloud removal methods at the two study sites (rRMSE: relative root mean square error; SSIM: structural similarity index measure; MNSPI: modified neighborhood similar pixel interpolator; GNSPI: geostatistical neighborhood similar pixel interpolator; GPR: Gaussian process regression; NIR: near infrared).The best metric is shown in bold.

Table 2 .
Band-wise accuracy statistics of different cloud removal methods at the two study sites (rRMSE: relative root mean square error; SSIM: structural similarity index measure; MNSPI: modified neighborhood similar pixel interpolator; GNSPI: geostatistical neighborhood similar pixel interpolator; GPR: Gaussian process regression; NIR: near infrared).The best metric is shown in bold.

Table 3 .
Accuracy statistics of different cloud removal methods for three vegetation indices prediction at the two study sites (VI: vegetation index; NDVI: normalized difference vegetation index; NDRE: normalized difference red-edge; NDWI: normalized difference water index).The best metric is shown in bold.