Virtual image-based cloud removal for Landsat images

ABSTRACT

The inevitable thick cloud contamination in Landsat images has severely limited the usability and applications of these images. Developing cloud removal algorithms has been a hot research topic in recent years. Many previous algorithms used one or multiple cloud-free image(s) in the same area acquired on other date(s) as reference image(s) to reconstruct missing pixel values. However, it remains challenging to determine the optimal reference image(s). In addition, abrupt land cover change can substantially degrade the reconstruction accuracies. To address these issues, we present a new cloud removal algorithm called Virtual Image-based Cloud Removal (VICR). For each cloud region, VICR reconstructs the missing surface reflectance by three steps: virtual image within cloud region construction based on time-series reference images, similar pixel selection using the newly proposed temporally weighted spectral distance (TWSD), and residual image estimation. By establishing two buffer zones around the cloud region, VICR allows automatic selection of the optimal set of time-series reference images. The effectiveness of VICR was validated at four testing sites with different landscapes (i.e. urban, croplands, and wetlands) and land change patterns (i.e. phenological change, abrupt change caused by flooding and tidal inundation), and the performances were compared with mNSPI (modified neighborhood similar pixel interpolator), WLR (weighted linear regression) and ARRC (AutoRegression to Remove Clouds). Experimental results showed that VICR outperformed the other algorithms and achieved higher Correlation Coefficients and lower Root Mean Square Errors in surface reflectance estimation at the four sites. The improvement is particularly noticeable at the sites with abrupt land change. By considering the difference in the contributions from the reference images, TWSD can select more reliable similar pixels to improve the prediction of abrupt change in surface reflectance. Moreover, VICR is more robust to different cloud sizes and to changing reference images. VICR is also computationally much faster than ARRC. The framework for time-series image cloud removal by VICR has great potential to be applied for large datasets processing.

KEYWORDS:

• Cloud removal
• Time-series images
• Landsat
• Similar pixels
• Gap filling
• Thick clouds

1. Introduction

Remote sensing imagery acquired by optical satellite sensors are important earth observation data sources and play irreplaceable roles in the studies of global climate change and natural resources management (Kokhanovsky 2013; Nijland, Reshitnyk, and Rubidge 2019; Wang et al. 2021; Chen et al. 2022). Among the optical satellite sensors, the Landsat series, including Thematic Mapper (TM), Enhanced Thematic Mapper (ETM+), Operational Land Imager (OLI) are possibly one of the most widely used owing to their medium spatial resolution, regular revisit cycle, and long-term data continuity (Huang et al. 2021; Zhu et al. 2020). However, cloud contaminations in Landsat images are inevitable because approximately one-third of global land area is obscured by clouds. For some areas, the probability of cloud cover is even higher (Shen et al. 2015; Cao et al. 2020). Clouds and their shadows lead to significant reductions in the availability of valid observations and restrict real-time monitoring of land surface change (Wang et al. 2021; Li et al. 2021). Compared to thin clouds, thick clouds result in missing image values as they completely block the electromagnetic energy reflected from the land surface. Reconstruction of the images covered by thick clouds and their shadows (hereafter we call cloud removal) has become a hot research topic for the purpose of wider application of time-series Landsat images (Shen et al. 2015; B. Chen et al. 2017; Luo et al. 2020; Zhang et al. 2020; Li et al. 2016).

The basic idea of cloud removal approaches is to predict unknown pixels in cloud/shadow areas with the assist of known information. According to the sources of auxiliary information, cloud removal methods can be divided into three categories: (1) single image-based methods, (2) multi-sensor image-based methods, and (3) single-sensor reference image(s)-based methods. Single image-based methods restore cloud pixels based on the assumption that the missing data and the surrounding cloudless data share the same statistical distribution or geometric structures. Typical methods include spatial interpolation (Siravenha et al. 2011; Li et al. 2014), and deep learning algorithms such as Generative Adversarial Networks (Xu et al. 2022; Zheng, Liu, and Wang 2021). These methods do not require additional data sources, while produce unreliable results for large clouds. The second category of methods uses auxiliary information from different sensors such as synthetic aperture radar (SAR) images that are less affected by cloud, or MODIS images that have short revisit cycle (Shen et al. 2019; Gao et al. 2021; Meraner et al. 2020). Due to radiometric inconsistencies and spatial resolution discrepancies, these methods have limited performances in areas with complex spatial patterns (Shen et al. 2019; Li et al. 2019). Single-sensor reference image(s)-based methods utilize cloud-free images acquired by the same sensor on other date(s) as reference image(s). The spectral and spatial relationships between the neighboring pixels and the cloud pixels (also called target pixels) are quantified based on the reference image(s) and used to fill the missing values (Chen et al. 2011; Zhu et al. 2012; Zeng, Shen, and Zhang 2013; Vuolo, Wai-Tim, and Atzberger 2017; Malek et al. 2018; Cao et al. 2020; Zhang et al. 2020; Hu et al. 2019). As images acquired by the same sensor (e.g. Landsat TM/ETM+/OLI) are not affected by radiometric or geometric inconsistencies, this type of methods has become the mainstream for thick cloud removal. Representative and commonly used methods include Neighborhood Similar Pixel Interpolator (NSPI) (Chen et al. 2011), modified neighborhood similar pixel interpolator (mNSPI) (Zhu et al. 2012), and weighted linear regression model (WLR) (Zeng, Shen, and Zhang 2013).

Most previous research relied on one cloudless image as the reference image (Zhu et al. 2012; Cheng et al. 2014; Zhang et al. 2018; Xu et al. 2022). The reconstruction accuracies depend on the similarity between the reference image and the cloud image. Generally, larger time interval between the reference image and the cloud image indicates weaker correlation between them (Lin et al. 2013; Du et al. 2019). Recent methods proposed utilization of multi-temporal reference images (B. Chen et al. 2017; Cao et al. 2020; Zhang et al. 2020; Zhou et al. 2022). For example, B. Chen et al. (2017) sorted multi-temporal reference images and utilized the most similar three reference images. Cao et al. (2020) proposed AutoRegression to Remove Cloud (ARRC) algorithm which uses multi-year (normally 3–4 years) Landsat surface reflectance time series to restore missing data. Zhang et al. (2020) proposed a spatio-temporal patch group deep learning framework based on multi-temporal images. Incorporation of multiple reference images effectively reduce the dependency on a single reference image. However, it is challenging to determine the optimal set of reference images because too much dependency on time series reference information might lead to adverse effect when abrupt land cover change occurred in the cloud image (Cao et al. 2020; Zhang et al. 2021).

Another key to the success of cloud removal algorithms is how to integrate the reference information to fill the missing values. Two types of strategies, i.e. pixel-based and patch-based strategies, have been proposed. Pixel-based strategy aims to predict values of each target pixel based on neighboring pixels. Similarity metrics between neighboring pixels and the target pixel are usually developed to select similar neighboring pixels (Zhu et al. 2012; Cheng et al. 2014; Chen et al. 2017; Cao et al. 2020; Padhee and Dutta 2019). However, it is difficult to determine whether the similar pixels are still similar to the target pixel at the time of cloud image acquisition. Abrupt land cover change can lead to obvious change in the target cloud pixel, thus the utilization of “similar” pixels during other time leads to reconstruction errors (Chen et al. 2018; Li et al. 2019). We expect that paying more attention to the reference images with greater similarities with the targe image might help to alleviate the errors in similar pixel selection. Patch-based strategy involves learning relationships between pairs of cloudless and the simulated cloud-contaminated image patches as training datasets. Deep learning models, such as spatio-temporal patch group deep learning framework have been developed (Zhang et al. 2020). However, to date the proposed deep learning models primarily aimed to restore textural and structure patterns in the target image rather than achieve high accuracy of the surface reflectance estimation of individual pixels (Zhang et al. 2020, 2021; Xu et al. 2022). For example, Xu et al. (2022) aimed to reconstruct Landsat Level-1 digital number images. The 16-bit images were down sampled to 8-bit to be better suitable for deep learning models. In addition, requirements for large number of training datasets and the complex parameter settings of deep learning restrict its wide applications for operational cloud removal.

In this study, we propose a new cloud removal algorithm for Landsat surface reflectance (SR) images that we call Virtual Image-based Cloud Removal (VICR). It combines virtual image construction and pixel-based prediction by using temporal, spectral, and spatial information in time-series reference images. The idea of virtual image is to produce a synthesis reference image within cloud region based on time-series reference images, which is closer to the image at the prediction time than any of the original reference images. Although the concept of virtual image was developed for spatio-temporal fusion of Landsat and MODIS data (Wang et al. 2020), it is introduced in cloud removal for the first time. In addition, VICR provides two functions different from previous cloud removal algorithms: (1) it allows automatic selection of the optimal set of temporal reference images; (2) it provides a new similarity index, namely temporal weighted spectral distance (TWSD), for similar pixel selection that aims to reduce prediction errors caused by temporal changes. We expect VICR to achieve more stable cloud removal performance in landscapes with complex spatial pattern and abrupt change than the previous cloud removal algorithms based on temporal reference image(s). In addition, we expect that VICR is computational efficient and can be practically applied for cloud removal for time-series images.

2. Virtual Image-based Cloud Removal (VICR) algorithm

VICR performs cloud removal for each cloud region separately. The basic idea of VICR is that Landsat SR image at band n within a cloud region can be reconstructed by the sum of a virtual image that is linearly transformed from time-series reference images and a residual image, which can be represented as the following equation:

(1) ${\widehat{L}}_{p\_cr}\left(n\right)=L_{v\_cr}\left(n\right)+{\gamma }_{cr}\left(n\right)$(1)

where ${\widehat{L}}_{p\_cr}\left(n\right)$ indicates the Landsat SR image within the cloud region at band n at the prediction time (thereafter called target time), $L_{v\_cr}\left(n\right)$ indicates the virtual SR image at band n within the cloud region, and ${\gamma }_{cr}\left(n\right)$ indicates the residual image within the cloud region.

VICR includes three main steps: (1) construct the virtual image $L_{v\mathrm{\_}cr}\left(n\right)$ in the cloud region by local linear transformation of time-series reference images; in this step, the optimal set of reference images are determined (Section 2.1); (2) for each pixel in the cloud region, select similar pixels in the neighborhood of the cloud region with TWSD (Section 2.2); (3) construct residual image ${\gamma }_{cr}\left(n\right)$ for cloud region based on similar pixels and then predict SR image ${\widehat{L}}_{p\_cr}\left(n\right)$ in the cloud region (Section 2.3) (Figure 1).

Figure 1. The flowchart of VICR.

2.1 Construct virtual image within cloud region

VICR constructs virtual image within cloud region for each band. For band n, $L_{v\mathrm{\_}cr}(n)$ is a linear transformation of time-series reference image in the same region. It can be expressed as:

(2) $L_{v\_cr}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_cr}(n)+b\left(n\right)$(2)

where $L_{t\_cr}(n)$ is the reference image in the cloud region acquired at time t. Note that $L_{t\_cr}(n)$ should be cloudless in the cloud region. T is the total number of reference images. $a_t\left(n\right)$ is the coefficient of the t^th reference image and $b\left(n\right)$ is a constant. Both parameters are to be estimated, and the optimal set of time-series reference images are to be determined.

2.1.1 Estimation of $a_t\left(n\right)$ and $b\left(n\right)$

To estimate $a_t\left(n\right)$ and $b\left(n\right)$, we assume that the image in the cloud region and its neighborhood image around the cloud region share the same transformation function. Therefore, we generate a buffer area outside the cloud region using dilation operation with a window size of 31 × 31 pixels (“imdilate” function in MATLAB). The buffer area has approximately 15 pixels width, and we call it “buffer zone 1” (buf1) (Figure 2). Next, $a_t\left(n\right)$ and $b\left(n\right)$ are estimated using the linear regression model fitted between time-series reference images within buf1 ($L_{t\_buf1}(n)$) and the image at the target time within buf1 ($L_{buf1}(n)$).

Figure 2. Establishment of buffer zone 1 (buf1), buffer zone 2 (buf2) around the cloud region.

(3) $L_{buf1}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf1}(n)+b\left(n\right)+{\gamma }_{buf1}\left(n\right)$(3)

where $a_t\left(n\right)$ and $b\left(n\right)$ are derived using the least squares method. ${\gamma }_{buf1}\left(n\right)$ represents the residual image of band $n$ within buf1.

2.1.2 Selection of the optimal set of time-series reference images

In the above process, it is clear that the number and the acquisition time of the reference images affect the accuracy of regression fitting. We expect that a set of reference images can be selected to minimize ${\gamma }_{buf1}\left(n\right)$ so that the virtual image estimated from Eq. 3(3) $L_{buf1}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf1}(n)+b\left(n\right)+{\gamma }_{buf1}\left(n\right)$(3) is as close as the real image at the target time. For this purpose, we build a new buffer area with the same window size outside buf1, which is called buffer zone 2 (buf2) (Figure 2). The virtual image at the target time within buf2 ($L_{v\_buf2}$) is obtained through the regression parameters $a_t\left(n\right)$, $b\left(n\right)$ and the reference images within buf2 ($L_{t\_buf2}(n)$) (Eq. 4(4) $L_{v\_buf2}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf2}(n)+b\left(n\right)$(4) ).

(4) $L_{v\_buf2}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf2}(n)+b\left(n\right)$(4)

We hope that $L_{v\_buf2}$ is very close to the real image at the target time $L_{buf2}$. We assume that the reference images that yield the most accurate estimation of $L_{v\_buf2}$ result in the highest fitting accuracy in buf1 as well as in the cloud region. We call these reference images the optimal set of reference images. In this study, we take an iterative strategy to select the optimal set of reference images.

Step 1. Select the image that is temporally closest to the target image and has no cloud/shadow in the cloud region as the first reference image (t₁ in Figure 3). In case that two images are selected and have the same temporal interval with the target image, the image acquired before the target date is preferred as the first reference image. Then, calculate coefficients $a_1\left(n\right)$ and $b\left(n\right)$ by Eq.3(3) $L_{buf1}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf1}(n)+b\left(n\right)+{\gamma }_{buf1}\left(n\right)$(3) .

Figure 3. An illustration of selection and sorting of reference images. The tick marks denote the image acquisition dates. The time interval between neighboring tick marks is 16 days. The red tick mark denotes the date when the target image was acquired. The images on the gray tick marks have cloud pixels in the cloud region, therefore they are not considered as reference images. The images on the blue tick marks are cloud free within the cloud region. The arrows show the sequence of image selection. It starts from the temporally closest image before target date (shown in green square arrowhead), and then visit the image on the other side of target date. In this example, the images on t₁, t₂, … , and t₇ are added as input in sequence.

Step 2. Calculate $L_{v\_buf2}(n)$ using the coefficients (Eq. 4(4) $L_{v\_buf2}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf2}(n)+b\left(n\right)$(4) ) and then calculate the root mean square error of $L_{v\_buf2}$ averaged over all bands ($\overline{RMSE}$) following Eq. 5(5) $\overline{RMSE}\text{=}\frac{\text{1}}{\text{N}}{\displaystyle {\sum }_{n=1}^N\sqrt{\frac{{\displaystyle {\sum }_{k=1}^K(L_{v\_buf2}(x_k,y_k,n)\text{-}{L}_{buf2}}(x_k,y_k,n){)}^2}{K}}}$(5) .

(5) $\overline{RMSE}\text{=}\frac{\text{1}}{\text{N}}{\displaystyle {\sum }_{n=1}^N\sqrt{\frac{{\displaystyle {\sum }_{k=1}^K(L_{v\_buf2}(x_k,y_k,n)\text{-}{L}_{buf2}}(x_k,y_k,n){)}^2}{K}}}$(5)

where $\left(x_k,y_k\right)$ is the position of the $k\mathrm{t}\mathrm{h}$ pixel and $K$ represents the total number of pixels in buf2. $N$ represents the total number of image bands.

Step 3. Select another image on the opposite side of the previous reference image date on the time axis (Figure 3). If the image is cloud free within the cloud region, then consider it as the second reference image and add it in the input of Eq. 3(3) $L_{buf1}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf1}(n)+b\left(n\right)+{\gamma }_{buf1}\left(n\right)$(3) (t₂ in Figure 3); otherwise, skip it and go to the other side of the image on the time axis. Calculate the virtual image within buf2 and its $\overline{RMSE}$. $R$epeat Step 3 until $\overline{RMSE}$ does not continue to decline. At this time, it is considered that the selected images are the optimal set of time-series reference images.

2.1.3 Generation of $L_{v\mathrm{\_}cr}(n)$ and ${\gamma }_{buf1}\left(n\right)$

Based on the optimal set of time-series reference images, we can generate the virtual image in the cloud region at target time $(L_{v\mathrm{\_}cr}(n))$ using Eq. 2(2) $L_{v\_cr}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_cr}(n)+b\left(n\right)$(2) and the residual image in buf1 (${\gamma }_{buf1}\left(n\right)$) using Eq. 3(3) $L_{buf1}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf1}(n)+b\left(n\right)+{\gamma }_{buf1}\left(n\right)$(3) .

2.2 Construct TWSD and select similar pixels

In sections 2.2 and 2.3, we aim to estimate residual image ${\gamma }_{cr}\left(n\right)$ for the cloud region. The residual of each pixel in the cloud region is estimated by weighted average of the residuals of its similar pixels in ${\gamma }_{buf1}\left(n\right)$. To capture the change in SR in temporal domain, we propose a new similarity index, i.e. temporal-weighted spectral distance (TWSD, Eq. 6(6) $\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}=\sqrt{\frac{{\sum }_{t=1}^T{\sum }_{n=1}^N\left|a_t\left(n\right)\right|\ast {\left[L_{t\mathrm{\_}cr}\left(x_i,y_i,n\right)-L_{t\mathrm{\_}buf1}\left(x_j,y_j,n\right)\right]}^2}{N}}$(6) ), in which both spectral and temporal distances between pixels are considered.

(6) $\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}=\sqrt{\frac{{\sum }_{t=1}^T{\sum }_{n=1}^N\left|a_t\left(n\right)\right|\ast {\left[L_{t\mathrm{\_}cr}\left(x_i,y_i,n\right)-L_{t\mathrm{\_}buf1}\left(x_j,y_j,n\right)\right]}^2}{N}}$(6)

where $\left(x_i,y_i\right)$ is the position of the $i$th pixel in the cloud region and $\left(x_j,y_j\right)$ is the position of the $j$th pixel in buf1. $\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}$ represents the distance between the $i$th pixel and the jth pixel in both spectral and temporal domain. Unlike previous spectral distance metrics (X. Zhu et al. 2012; B. Chen et al. 2017; Cao et al. 2020), we incorporate the transformation coefficient $a_t\left(n\right)$in the calculation of the distance. Larger $a_t\left(n\right)$ indicates that the corresponding reference image has a higher contribution to the target image, thus the neighboring pixels at time t should be given greater weights. For each pixel in the cloud region, we select the 20 pixels with the lowest $\mathrm{T}\mathrm{W}\mathrm{S}\mathrm{D}$ (i.e. M = 20) as the similar pixels of the target pixel. These similar pixels are used to estimate the residuals of the target pixel and then to obtain the residual image in the cloud region.

2.3 Estimate residual image and recover cloud image

For a target pixel $\left(x_i,y_i\right)$ in the cloud region, its residual is estimated using linear weighted sum approach considering the spatial and spectral-temporal distances to its similar pixels. The spatial distance $\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{ij}$ between the target pixel $i$ and its similar pixel $j$ is calculated as their Euclidean distance.

(7) $\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{ij}=\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}$(7)

We use $\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}$ to represent the spectral-temporal distance between pixel $i$ and $j$. Then, $\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{ij}\hspace{0.17em}$ and $\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}$ are normalized through:

(8) $D_{ij}=\frac{\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{ij}-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_i\right)}{max\left(\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_i\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_i\right)}+1$(8)

(9) $S{T}_{ij}=\frac{\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i\right)}{max\left(\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i\right)}+1$(9)

where $D_{ij}$ and $S{T}_{ij}$ represent the normalized spatial distance and spectral-temporal distance respectively. $\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_i$ and $\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i$ represent the arrays composed of the spatial distances [$\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{i1},\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{i2},\dots ,\mathrm{D}\mathrm{i}\mathrm{s}{\mathrm{t}}_{iM}$] and spectral-temporal distances [$\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{i1},\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{i2},\dots ,\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{iM}$] between the target pixel $i$ and the M similar pixels, respectively (M = 20). The $max\left(\cdot \right)$ and $\mathrm{m}\mathrm{i}\mathrm{n}\left(\cdot \right)$ are to find the maximum and minimum values of the array. The weight $w_{ij}$ of the similar pixel $j$ of the target pixel $i$ is calculated as:

(10) $w_{ij}=\frac{1/\left(D_{ij}S{T}_{ij}\right)}{{\sum }_{j=1}^{M}1/\left(D_{ij}S{T}_{ij}\right)}$(10)

Finally, for the cloud region we obtain the residual value of the target pixel i by weighted summation of the residuals of the similar pixels in buf1.

(11) ${\gamma }_{cr}\left(x_i,y_i,n\right)\text{=}{\displaystyle {\sum }_{j=1}^M{w}_{ij}}\times {\gamma }_{buf1}\left(x_j,y_j,n\right)$(11)

When the residual image of the cloud area ${\gamma }_{cr}\left(n\right)$ is constructed, it is added to the virtual image $L_{v\_cr}$ obtained in Eq. 1(1) ${\widehat{L}}_{p\_cr}\left(n\right)=L_{v\_cr}\left(n\right)+{\gamma }_{cr}\left(n\right)$(1) to predict the cloud free image in the cloud region ${\widehat{L}}_{p\_cr}$. Following the above procedures, we reconstruct each cloud region separately. Once the images in all cloud regions are reconstructed, they are combined to obtain the final cloud-free image scene.

3. Experimental data and evaluation

3.1 Testing sites

Four sites with different landscape heterogeneities and land change patterns were selected to evaluate the VICR method (Figure 4). The first site is at the Jianghan Plain in Hubei Province, China (112.8°E, 30.2°N). This site is an important agricultural production area in Central China. The crop types are dominated with rapeseed and paddy rice rotations (Fang et al. 2021) and the croplands are fragmented, presenting heterogeneous spatial patterns. The second site is in the Lower Gwydir Catchment, Australia (149.36°E, 29.08°S). Land use include grazing, dryland cropping, horticulture, and irrigated lands (Shi et al. 2014). A flood event occurred in mid-December 2004, resulting in an abrupt land cover change shown on the Landsat 5 TM image on December 12^th, 2004. The third site covers urban and suburban areas in Beijing (116.45°E, 39.95°N), China, presenting great spatial heterogeneity and vegetation phenological dynamics. The fourth site is at the coastal wetland in Yellow River Delta in Shandong Province, China (119.2°E, 37.75°N). As an area with intense land sea interactions, land surface reflectance changes frequently due to tidal inundation and sediment accumulation. In recent years, the landscape patterns in the Yellow River Delta have changed dramatically due to rapid expansion of an invasive plant species Spartina alterniflora (Z. Wang et al. 2021).

Figure 4. Actual images and cloud-simulated images in standard false color composition over (A) Jianghan Plain (Landsat 8 OLI), (B) Lower Gwydir Catchment (Landsat 5 TM), (C) Beijing (Landsat 8 OLI), and (D) Yellow River Delta (Landsat 8 OLI). Acquisition dates are in the parentheses in the format of MM/DD/YY. The image sizes for the testing sites A-C are 1200 × 1200 pixels, and that of testing site D is 1397 × 1234 pixels. The number of simulated cloud pixels in (A-D) are 93,312, 107,529, 89813, and 103,174, respectively.

3.2 Landsat data and preprocessing

To evaluate the VICR method, we selected a cloud-free Landsat imagery at each site and simulated a cloud region (hereafter we call target image) (Figure 4). The cloud region has around 90,000 ~ 110,000 pixels. Time-series Landsat images with cloud cover less than 80% acquired before and after the target date were considered as candidate reference images. All Landsat images were Collection 1 Tier 1 Level 2 surface reflectance products downloaded from Google Earth Engine. For each Landsat image, we used the quality assurance band to detect cloud and cloud shadow regions and mask them. We did not use Landsat 7 ETM+ images because the images had missing strips after the scan-line corrector failed in 2003.

3.3 Assessment metrics

The performance of VICR was evaluated by comparing the reconstructed images and the actual image in cloud region. Quantitative assessment was conducted using Correlation Coefficient (CC), Root Mean Square Error (RMSE) and Structural Similarity Index Measure (SSIM). CC represents the degree of agreement between the estimation and the real value. RMSE measures the statistical difference of pixel values between the predicted image and the actual image. SSIM evaluates the overall structure similarity between the two images (Zhou et al., 2004). The quantitative evaluation was carried out for six bands of Landsat OLI/TM sensors, i.e. blue, green, red, near infrared red (NIR), short-wave infrared 1 (SWIR1), and short-wave infrared 2 (SWIR2) bands.

4. Experiment design and results

Six experiments were designed to evaluate the performances of VICR method from multiple aspects. The specific experiment design and their purposes of each experiment are shown in Table 2. We selected mNSPI (modified neighborhood similar pixel interpolator), WLR (weighted linear regression), and ARRC (AutoRegression to Remove Clouds) for comparison because they all involved using both temporal and spatial auxiliary information in reference image(s) for cloud removal. The source codes are publicly available except ARRC (mNSPI: https://xiaolinzhu.weebly.com; WLR: http://sendimage.whu.edu.cn/send-resource-download). mNSPI and WLR are based on a single reference image and have been widely used (Zhu et al. 2012; Zeng, Shen, and Zhang 2013). mNSPI is provided in IDL code, and WLR is provided with an executable program. ARRC is a recently-proposed cloud removal method based on time series Landsat images (Cao et al. 2020). It combines predictions from a long-term component using autoregression of time series observations (3 years of Landsat images were used as recommended by ARRC) and a short-term component based on a single cloudless reference image. As the source code is not available, we rewrote ARRC on MATLAB. VICR was also coded on MATLAB. All processes in this study were carried out on a workstation (Dell-T8520, CPU: i9 -10,900× 3.70 GHz, RAM: 32GB). The code for VICR is available at https://github.com/wangzp-ludan/VICR.

Table 1. Reference images and their input orders for VICR in Jianghan Plain, Lower Gwydir Catchment, Beijing and Yellow River Delta.

Jianghan Plain		Lower Gwydir Catchment		Beijing		Yellow River Delta
Reference image	Order	Reference image	Order	Reference image	Order	Reference image	Order
11/14/14	11	08/06/04	11	03/04/17	11	05/25/17	12
12/16/14	10	08/22/04	10	04/21/17	8	07/12/17	10
01/17/15	8	09/23/04	8	05/07/17	6	09/30/17	9
03/22/15	6	10/09/04	6	05/23/17	4	01/20/18	6
07/28/15	3	10/25/04	4	07/10/17	2	03/09/18	5
Target image: 09/14/15		11/26/04**	1	Target image: 09/12/17		06/29/18**	2
10/16/15*	1	Target image: 12/12/04		09/28/17*	1	Target image: 07/31/18
11/01/15**	2	12/28/04*	2	10/30/17**	3	08/16/18*	1
02/05/16	4	01/13/05	3	11/15/17	5	09/17/18	3
05/11/16	5	01/29/05	5	12/01/17	7	10/03/18	4
04/09/16	7	02/14/05	7	12/17/17	9	03/12/19	7
08/31/16	9	03/02/05	9	01/02/18	10	04/13/19	8
12/05/16	12	04/03/05	12	02/03/18	12	05/31/19	11

Note: Image date in MM/DD/YY.

*represents the cloudless reference image for mNSPI, WLR and ARRC short-term component in Experiment III.

** represents the cloudless reference image for mNSPI, WLR and ARRC short-term component when the image * is not available in Experiment IV.

Table 2. Specific information on the six experiments.

Experiment design	Purpose
Experiment I: Assessment on the selection of optimal set of reference images	To verify the reliability of the automatic selection of optimal reference images.
Experiment II: Assessment on TWSD	To evaluate the effectiveness of the Temporal Weighted Spectral Distance (TWSD) metrics.
Experiment III: Quantitative evaluation and comparison with other cloud-removal methods	To compare the accuracies of VICR with ARRC, mNSPI and WLR algorithms to verify its advantage.
Experiment IV: Sensitivity to reference images	To evaluate the sensitivity of VICR accuracy to reference image compared to ARRC, mNSPI and WLR.
Experiment V: Sensitivity to different cloud sizes	To analyze the influence of different cloud sizes on the stability of the accuracies.
Experiment VI: Performance on real cloud-contaminated images	To use real cloud-contaminated images for cloud removal and compare the reconstruction effects with ARRC, mNSPI and WLR.

4.1 Experiment I: Assessment on the selection of optimal set of reference images

4.1.1 Experiment design

VICR uses the fitting parameters obtained from buf1 to predict virtual image in buf2 ($L_{v\_buf2}$) at the target time, and the RMSE between the $L_{v\_buf2}$ and the actual image in buf2 ($L_{buf2}$) is used to guide determination of the optimal set of time-series reference images. The basic assumption is, with changing reference images the prediction errors in the cloud region present the same variation patterns as the $\overline{RMSE}$s of $L_{v\_buf2}$, because buf2 and the cloud region are spatially close. To justify the assumption, we conducted VICR based on the input of 1 to 12 reference images. The reference images were selected and added to the reference image list following the procedures in Section 2.1 (Figure 3), and the order of the input images is shown in Table 1. The $\overline{RMSE}$ of $L_{v\_buf2}$, the $\overline{RMSE}$ of the virtual image in the cloud region ($L_{v\_cr}$), and the $\overline{RMSE}$ of the reconstructed image in the cloud region (${\widehat{L}}_{p\_cr}$) were calculated for four testing sites, respectively.

4.1.2 Experiment results

In general, the $\overline{RMSE}$s of $L_{v\_buf2}$, $L_{v\_cr}$ and ${\widehat{L}}_{p\_cr}$ showed consistent pattern with increasing number of reference images (Figure 5). They all presented decreasing trend with first a few reference images and then gradually stabilize. For Beijing, the $\overline{RMSE}$s of $L_{v\_buf2}$, $L_{v\_cr}$ and ${\widehat{L}}_{p\_cr}$ all reached the minimum when three reference images were used as input (Figure 5c). For Jianghan Plain, Lower Gwydir Catchment and Yellow River Delta, the $\overline{RMSE}$ of $L_{v\_buf2}$ stoped decreasing when the 7^th, 6^th and 6^th reference image was added, respectively (Figure 5a,b,). Thus, the first 7, 6, and 6 images were selected as the optimal set of reference images for the three sites, respectively. With these reference images, unfortunately, the $\overline{RMSE}$ of ${\widehat{L}}_{p\_cr}$ was not the minimum. For Jianghan Plain, Lower Gwydir Catchment and Yellow River Delta, the $\overline{RMSE}$ of ${\widehat{L}}_{p\_cr}$ reached the minimum with 4, 11 and 10 reference images, respectively (Figure 5a,b,). However, little difference was found between the minimum $\overline{RMSE}$ and the $\overline{RMSE}$ resulted from the selected set of reference images (1.43 × 10⁻⁴ for Jianghan Plain, 0.3 × 10⁻⁴ for Beijing, and 4.5 × 10⁻⁴ for Yellow River Delta), indicating that good accuracies are still achieved with these reference images.

Figure 5. $\overline{RMSE}$s of the virtual image in the buffer zone 2 ($L_{v\_buf2}$), the virtual image in the cloud region ($L_{v\_cr}$), and the reconstructed image in the cloud region ${\widehat{L}}_{p\_cr}$ resulted from one to twelve reference images. The minimum RMSEs of $L_{v\_buf2}$ and ${\widehat{L}}_{p\_cr}$ are marked with triangles.

Figure 5 shows that the first few reference images temporally close to the target date had more contribution than those images farther from the target date. The $\overline{RMSE}$s first decreased substantially and then had little change or even increased slightly. It seems that more reference images do not necessarily lead to higher accuracies. For instance, for Beijing the $\overline{RMSE}$s of ${\widehat{L}}_{p\_cr}$ showed slight increase when more than 3 reference images were used as input; for Lower Gwydir Catchment, the $\overline{RMSE}$ from 12 reference images was higher than that from 6 and 11 reference images. In addition, we observe that the $\overline{RMSE}$ of ${\widehat{L}}_{p\_cr}$ was much lower than that of $L_{v\_cr}$ regardless of the number of reference images, suggesting that the estimation of residual image based on similar pixels is critical for the estimation of ${\widehat{L}}_{p\_cr}$.

4.2 Experiment II: Assessment on TWSD

4.2.1 Experiment design

VICR uses TWSD to select similar pixels in buf1 and to estimate residuals for the target pixels in cloud region. The newly proposed TWSD considers the spectral similarity between pixels in temporal domain. To test the efficacy of TWSD, we compared the TWSD to spectral distance ($SD$) that has been widely used for similar pixel selection in cloud removal algorithms such as NSPI, mNSPI, STWR and the ARRC short-term component. Apart from similar pixel selection, these algorithms also used $SD$ to quantify the contribution of the similar pixels to the target pixel. In this experiment, we designed four scenarios.

Scenario 1: $SD$ was used for both similar pixel selection and residual allocation (substitute TWSD in Eq. 9(9) $S{T}_{ij}=\frac{\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_{ij}-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i\right)}{max\left(\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{T}\mathrm{W}\mathrm{S}{\mathrm{D}}_i\right)}+1$(9) with $SD$) based on virtual images in cloud region and buf1. The formulation of $SD$ is listed in Eq. 12(12) $S{D}_{ij}\text{=}\sum n=1N{L}_{v\_cr}\left(x_i,y_i,n\right)-L_{v\_buf1}\left(x_j,y_j,n\right)2N$(12) .

Scenario 2: $SD$ was used to select similar pixels, and TWSD was used to allocate residuals.

Scenario 3: TWSD was used to select similar pixels, and $SD$ was used to allocate residuals.

Scenario 4: TWSD was used for both similar pixel selection and residual allocation (i.e. VICR).

(12) $S{D}_{ij}\text{=}\sum n=1N{L}_{v\_cr}\left(x_i,y_i,n\right)-L_{v\_buf1}\left(x_j,y_j,n\right)2N$(12)

For each site and each band of Landsat image, the cloud-removed image in cloud region was compared with the actual image, and RMSE was calculated. By comparing the RMSEs of reconstruction results, the role of TWSD in similar pixel selection and residual redistribution was analyzed.

4.2.2 Experiment results

Figure 6 shows that TWSD for both similar pixel selection and residual allocation (Scenario 4) achieved the lowest RMSE for the four testing sites (except for the NIR and SWIR1 image in Beijing), while Scenario 1 produced the highest RMSE. When TWSD was used for similar pixel selection (Scenario 3 and 4), the reconstruction accuracies were significantly improved regardless which metrics were used for residual allocation (Scenario 3 vs. Scenario 1; Scenario 4 vs. Scenario 2). When $SD$ was used for similar pixel selection, TWSD used for residual estimation slightly improved the reconstruction accuracy (Scenario 2 vs. Scenario 1). Compared to the other three testing sites, the Lower Gwydir Catchment which experienced abrupt land change witnessed greater accuracy improvement when TWSD was used for similar pixel selection (Figure 6b). This suggests that the similarity between similar pixels and target pixels played a critical role in reconstruction, especially for the images with great land change.

Figure 6. Rmses of the reconstructed cloud image at each band for (a) Jianghan Plain, (b) Lower Gwydir Catchment, (c) Beijing and (d) Yellow River Delta at four scenarios where spectral distance (SD) and/or temporally weighted spectral distance (TWSD) for similar pixel selection and/or residual allocation.

4.3 Experiment III: Quantitative evaluation and comparison with other cloud-removal methods

4.3.1 Experiment design

To evaluate the performance of VICR, we compared the reconstruction accuracies of VICR with those of mNSPI, WLR, and ARRC. We conducted two sub-experiments. The first one was performed on the cloud-simulated images in Table 1. For mNSPI, WLR, and ARRC short-term component, we used the cloud free image that was closest to the cloud contaminated image as reference image (Table 1, Fig. S1). The second sub-experiment was performed for each site on another cloud simulated image (Fig. S2). On these images, three cloud regions with a total of over 120, 000 pixels were generated and randomly distributed on the image (Fig. S2). For both sub-experiments, 20 similar pixels were selected for VICR, mNSPI, WLR and ARRC short-term component for fair comparison. CC, RMSE and SSIM in the cloud regions were used as quantitative metrics. For the second sub-experiment, the values of the quantitative metrics for the three cloud regions were averaged.

4.3.2 Experiment results

Due to space limitation, here we display the results of the first sub-experiments. The results of the second sub-experiment are shown in supplementary materials (Fig. S2 and Fig. S3). Generally, the reconstructed images from VICR are visually more similar to the actual image than those from mNSPI, WLR and ARRC (Figure 7). For example, for the patchy croplands in Jianghan Plain, the cloud-removed images from WLR and ARRC showed poorer spatial heterogeneities and obviously smoother effects than VICR. Some fallow croplands were not well reconstructed. Both VICR and mNSPI better restored spatial details, while VICR obviously showed higher spectral similarities with the actual image (Figure 7a). For Lower Gwydir Catchment site where flooding caused significant land cover change, the cloud-removed images from mNSPI, WLR and ARRC were more similar to the reference image (Fig. S1). In comparison, the predicted image from VICR showed better spectral similarity with the actual image in part of the flood-inundated region (Figure 7b). For Beijing site where the land change in urban areas was relatively small, the results from the four methods were close (Figure 7c). For the Yellow River Delta, the results from mNSPI, WLR and ARRC demonstrated obvious spectral and texture discrepancies from the actual image and the tidal creeks were not well restored. In contrast, VICR generated more visually satisfactory results where tidal creeks and tidal flats kept intact (Figure 7d).

Figure 7. Visual comparisons between mNSPI, WLR, ARRC, and VICR the cloud-simulated images over (a) Jianghan Plain, (b) Lower Gwydir Catchment, (c) Beijing and (d) Yellow River Delta.

Quantitative assessment results are shown in Figure 8. In general, VICR presented much better performances (lower RMSE, higher CC and SSIM) than mNSPI and WLR for all sites at all bands. The decrease in RMSEs ranged from 0.0014 to 0.0223 compared to mNSPI, and ranged from 0.0017 to 0.0193 compared to WLR. For all sites, VICR performed better than ARRC in terms of CCs and RMSEs, especially for the sites with abrupt land change. For Lower Gwydir Catchment, VICR showed considerably higher CCs and lower RMSEs at red (0.0368 vs. 0.0400), green (0.0312 vs. 0.0334), SWIR1 (0.0422 vs. 0.0463) and SWIR2 (0.0302 vs. 0.0342) bands, indicating that VICR better captured the flood inundation than ARRC (Figure 8b). At the SWIR1 and SWIR2 bands, VICR showed higher SSIMs than ARRC. For the Yellow River Delta, the improvements in CCs and RMSEs were more significant compared to the other sites (Figure 8d). For example, CCs of VICR at NIR, SWIR1 and SWIR2 bands were 0.9507, 0.9221 and 0.9124, respectively, while those of ARRC were 0.9260, 0.8995 and 0.8962. The RMSEs of VICR at NIR, SWIR1 and SWIR2 bands were 0.0234, 0.0255 and 0.0169 respectively, while those of ARRC were 0.0286, 0.0290 and 0.0185, respectively. In most cases, VICR achieved comparable or slightly higher SSIM than ARRC. For Beijing and Jianghan Plain, VICR produced slightly lower SSIMs than ARRC at SWIR1 and SWIR2 bands (Figures 8a and 8c). In Beijing site, SSIM values of ARRC at SWIR1 and SWIR2 bands were 0.9485 and 0.9514 respectively, while those of VICR were 0.9454 and 0.9485 respectively. The second sub-experiments obtained similar comparison results as the first sub-experiment (Fig. S2 and Fig. S3), i.e. VICR achieved much better performance than mNSPI and WLR at all sites and all image bands, and in most cases VICR achieved higher CC, SSIM and lower RMSE than ARRC.

Figure 8. Comparisons of CC (left column), RMSE (center column) and SSIM (right column) between mNSPI, WLR, ARRC, and VICR for the cloud-simulated images over (a) Jianghan Plain, (b) Lower Gwydir Catchment, (c) Beijing and (d) Yellow River Delta.

4.4 Experiment IV: Sensitivity to reference images

4.4.1 Experiment design

We hope that a robust cloud removal method is little affected by the selection of reference images. To evaluate the sensitivity of mNSPI, WLR, ARRC, and VICR to different reference images, we conducted two sub-experiments by altering reference images and analyzing the change in reconstruction performances. The first sub-experiment assumes that the reference images in Experiment III were not available. For mNSPI, WLR, and ARRC short-term component, we replaced the reference image in Experiment III (images with * in Table 1 and Fig. S1) with another cloud-free image (images with ** in Table 1 and Fig. S1). For ARRC long-term component and VICR, we removed the reference image from the time-series reference image list. Note that the time intervals between the target image and the new reference image were longer than or same as the original reference image. In the second sub-experiment, we assumed the revisit period of Landsat satellite is 32, 48 and 64 days and thus lengthened the time interval between the reference image and the target image, and among the time-series reference images. For each site we predicted the image in the cloud region and compared the performances of the four algorithms.

4.4.2 Experiment results

Figure 9 displays the results from the first sub-experiment. In Jianghan Plain, Lower Gwydir Catchment and Beijing, the RMSEs from all algorithms increased at most bands when using the new reference image (Figure 9). Compared to mNSPI and WLR, ARRC and VICR obviously showed smaller increment of RMSEs, especially at the SWIR1 and SWIR2 bands for Lower Gwydir Catchment (Figure 9b), and at all bands for Beijing (Figure 9c). For the Beijing site, assuming that the reference image on 09/28/17 (MM/DD/YY) was not available, the next available reference image was on 10/30/17. From the target date (09/12/17) to the end of October, the land surface reflectance had changed a lot due to crop harvesting and vegetation phenological change (Figure S1A). In Lower Gwydir Catchment, although the time intervals between the two reference images (11/26/04 and 12/28/04) and the target image (12/12/04) were both 16 days, the accuracies of mNSPI and WLR decreased significantly because flood event occurred during mid-December. This indicates that large surface reflectance change has substantial negative effects on the performances of mNSPI and WLR.

Figure 9. Comparisons of RMSEs between mNSPI, WLR, ARRC and VICR for the cloud-simulated images for (a) Jianghan Plain, (b) Lower Gwydir Catchment, (c) Beijing and (d) Yellow River Delta in Experiment IV (colored bars) and in Experiment III (transparent bars with black border).

In Yellow River Delta, it is interesting that using the reference image temporally farther from the target date (more than a month from 06/29/18 to 07/31/18) achieved better performance than the original reference image (08/16/18) for all methods (Figure 9d). mNSPI, WLR and ARRC all showed considerable decrease in RMSEs at all bands, while the VICR presented more stable performance. Examination of the reference images showed that the wetlands at this site experienced tidal inundation due to high tide level on 08/16, resulting in great surface reflectance change. The reference image on 06/29 is spectrally more like the target image. Although ARRC long-term component utilized time-series Landsat images, its short-term component is still based on a single reference image. When surface reflectance experienced abrupt change, the short-term component is given higher weight when combining with long-term components (Cao et al. 2020). As a result, changing the reference image still had considerable impacts on the accuracies of ARRC in the cases where surface reflectance has abrupt change. Compared to ARRC, the change in RMSEs of VICR was negligible in Yellow River Delta (Figure 9d). In addition, even if the time interval between the reference image and the target image becomes longer, the VICR still ensured the lowest RMSE compared to other methods for all four sites.

Figure 10 shows the RMSEs averaged over all bands ($\overline{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$) from the second sub-experiment. All four algorithms showed increasing $\overline{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$ with increasing time intervals among images, while VICR generally showed lower magnitude of $\overline{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$ variation (except at Jianghan Plain). The reconstruction errors of VICR were significantly lower than mNSPI and WLR, and slightly lower than ARRC (except for revisit cycle = 32 days at Jianghan Plain). In Lower Gwydir Catchment and Yellow River Delta, VICR presents obvious improvement in accuracies (Figures 10b and 10d).

Figure 10. Comparisons of $\overline{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}$s between mNSPI, WLR, ARRC and VICR for the cloud-simulated images for (a) Jianghan Plain, (b) Lower Gwydir Catchment, (c) Beijing and (d) Yellow River Delta with revisit periods of 16, 32, 48 and 64 days.

4.5 Experiment V: Sensitivity to different cloud sizes

4.5.1 Experiment design

Intuitively, larger cloud region suggests less auxiliary information in the surrounding area can be provided for the reconstruction, especially for the central area of the cloud region. To evaluate the sensitivity of VICR to different cloud sizes, we generated clouds with sizes of 50 × 50, 100 × 100, 150 × 150, 200 × 200, 250 × 250, 300 × 300, 350 × 350, and 400 × 400 pixels (Figure 11) and reconstructed images in the cloud region based on the same reference images in Experiment III. For each method and cloud size, RMSEs at each band in the center 50 × 50 pixels was calculated and compared. Note that for Yellow River Delta, we used the image on 06/29/18 as reference image for mNSPI, WLR and ARRC short-term component to avoid the impacts of tidal inundation and to ensure good results.

Figure 11. Simulated clouds with 50×50, 100×100, 150×150, 200×200, 250×250, 300×300, 350×350 and 400×400 pixels on the actual image in the four sites.

4.5.2 Experiment results

Figure 12 illustrates RMSEs at NIR band of the four methods for different cloud sizes. RMSEs at other bands are shown in Fig. S4-S8 due to space limitation. VICR achieved considerably lower RMSE than mNSPI and WLR at all bands regardless of the cloud size for Jianghan Plain, Beijing and Yellow River Delta. In most cases, VICR presented slightly lower RMSEs than ARRC. In addition, RMSEs of VICR showed smaller fluctuations than ARRC with varying cloud sizes. This was particularly noticeable at NIR band in Lower Gwydir Catchment (Figure 12b), and at SWIR1 and SWIR2 bands in Yellow River Delta (Fig. S6). In Lower Gwydir Catchment, the RMSEs of all algorithms showed upward trend with increasing cloud size (Figure 12b). As the cloud covered the area with flood inundation, smaller cloud region indicated more flood-inundated neighboring pixels can be found. In an extreme case that the cloud covers the entire flood-inundated area, none of the neighboring pixels is spectrally similar to the target region, which certainly results in low reconstruction accuracies. Compared to mNSPI and WLR, VICR presented smaller increment of RMSE at SWIR1 and SWIR2 bands (Fig. S6), which are both representative bands showing water signal. Likewise, in Yellow River Delta, VICR achieved noticeably lower and more stable RMSEs than other methods at SWIR1 and SWIR2 bands for all cloud sizes.

Figure 12. The RMSEs of the reconstructed images from mNSPI, WLR, ARRC, and VICR at NIR band for cloud sizes ranging from 50×50 pixels to 400×400 pixels. Note that the RMSEs were calculated for the 50×50 pixels in the center of the clouds.

4.6 Experiment VI: Performance on real cloud-contaminated images

4.6.1 Experiment design

In this experiment, we performed cloud removal for real cloud-contaminated images and visually compared the reconstruction results of different methods. Considering practical applications, we also used VICR to remove clouds on the time-series Landsat 8 OLI SR images during 2015–2019 at Jianghan Plain and during 2016–2020 at Yellow River Delta. We calculated Normalized Difference Vegetation Index (NDVI) based on each reconstructed image and compared the time series NDVIs with the NDVIs from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices Version 6 product (MOD13Q1). We sorted the time-series Landsat images in order of the cloud cover percentage and processed them one by one from low to high cloud cover. Each cloud-removed image was then used as the reference image for the images with larger cloud cover. We call this process “time-series cloud removal procedure.”

4.6.2 Experiment results

As illustrated in Figure 13, the images reconstructed by VICR generally recovered better spatial details and maintained more realistic spectral characteristics compared to other algorithms. In Jianghan Plain, the reconstructed image by VICR presented greater spectral heterogeneity in cropland patches compared to other algorithms. The boundaries between the cropland patches were more visible in the image reconstructed by VICR. In Lower Gwydir Catchment, the images reconstructed by VICR, mNSPI and WLR presented similar spatial details, while that by ARRC presented smoother effect. In Yellow River Delta, it is obvious that the QA band layer detected larger clouds/cloud shadow areas than those in the actual image. The zoomed-in actual image covers a Spartina Alterniflora patch (Wang et al. 2021) which is partially contaminated by clouds in the southern coast of Yellow River Delta (yellow box in Figure 13d). The reconstructed image by VICR restored more accurate spectral information than the other methods as it is spectrally more similar to the actual image.

Figure 13. Visual effects of reconstructed images from VICR, mNSPI, WLR and ARRC on real cloud-contaminated images over (a) Jianghan Plain, (b) Lower Gywdir Catchment, (c) Beijing and (d) Yellow River Delta.

Figure 13. (Continued).

Figure 14 shows that the time series images recovered by VICR effectively supplemented the number of NDVI observations. The number of observations increased from 31 to 47 during 2015–2019 at the pixel in Jianghan Plain. The NDVI time series presented phenological dynamics well at both sites. In Jianghan Plain, rapeseed was normally harvested around May and then rice was planted. The crop rotation is observed in MOD13Q1 time series as it showed a local peak around March and a higher local peak around August. The local peak around March 2016 is not observed from the original Landsat NDVI time series but is well captured by the cloud-removed image (Figure 14a). Likewise, the growing peak of paddy rice in August 2018 was better captured by the cloud-removed images (Figure 14a). For the Phragmites australis marsh wetland pixel in Yellow River Delta, the reconstructed NDVI time series demonstrated greater variations than MOD13Q1 (Figure 14b). Negative NDVIs were obtained during December to February because the wetland pixel was covered by water in leaf-off season. In addition, the reconstructed NDVI time series better demonstrated growing and senescence patterns of Phragmites australis compared to the original time-series NDVI (e.g. in 2019).

Figure 14. NDVIs from cloud-free and cloud-removed Landsat images at (a) the pixel (30.2477°N, 112.8263°E) at a crop field in Jianghan Plain and (b) the pixel (37.8163°N, 119.0197°E) at Phragmites australis marsh wetland in Yellow River Delta. For comparison, MOD13Q1 NDVI time-series smoothed by spline function are plotted as blue line.

5. Discussion

VICR is different from previous algorithms in that it is based on virtual image construction. It is composed of patch-based virtual image prediction based on time-series reference images and pixel-based residual allocation incorporating spectral, temporal, and spatial information. Advantages of this algorithm are discussed in Section 5.1-5.4. The uncertainties are discussed in Section 5.5.

5.1 Virtual image construction for optimal reference images selection

The first improvement of VICR is that it allows automatic selection of optimal set of reference images. It is widely recognized that the selection of reference image has great impacts on the performance of cloud removal. Generally, the temporally closer reference image has more similar spectral characteristics to the target image so that higher reconstruction accuracy can be achieved (Zhu et al. 2012; Wang et al. 2021). However, this might not be true when sudden land change occurred. In the coastal wetlands in Yellow River Delta, the closest reference image resulted in lower accuracies of mNSPI, WLR and ARRC than using reference image with longer time interval to the target image. Several approaches have been proposed to address this issue. Some studies considered using CC or SSIM to select the most similar reference image (Lin et al. 2013; Kalkan and Derya Maktav 2018), but it was also proved that this method did not improve the reconstruction accuracy considerably (Cao et al. 2020). STWR algorithm sorted the reference image patches based on the spectral similarity and selected the most similar three patches as reference images (Chen et al. 2017). It is hard to determine whether the three selected images can achieve the highest accuracy. Another solution is to incorporate long-term time-series reference images such as the ARRC long-term component. However, the ARRC long-term component is advantageous in capturing gradual temporal patterns of a pixel rather than capturing abrupt change (Cao et al. 2020). Our algorithm balances between ARRC and those algorithms depending on a fixed number of reference images because it can flexibly select the optimal set of reference images by virtual image construction strategy. In virtual image construction, the coefficient a_t represents the contribution of each the reference image to the target image. As shown in Fig. S9, when the second reference image was added for Lower Gwydir Catchment and Yellow River Delta, the a_t value of the first reference image has a substantial decrease and becomes lower than that of the second image, indicating the second reference image has greater contribution in prediction. The varying a_t values when adding input reference images ensure that the land change pattern can be captured. In addition, based on our results in Section 4.1, we recommend that the optimal number of reference images for each cloud region need not exceed 12 because more reference images at this time can no longer significantly improve the reconstruction accuracy.

The idea of virtual image refers to the virtual image pair-based spatio-temporal fusion (VIPSTF) algorithm proposed by Wang et al. (2020), which constructs virtual Landsat-MODIS image pairs to predict image with the same spatial resolution as Landsat and the same temporal resolution as MODIS. The virtual image pair was generated by transplanting the linear relationship between time-series MODIS images to Landsat images. It was demonstrated theoretically that the virtual image pair is closer to the data at the prediction time than the known Landsat-MODIS image pairs (Wang et al. 2020). Inspired by this concept, in this study we constructed virtual image within cloud region to provide an initial prediction of the cloud image. Although the equation of building virtual image (Eq. (3))(3) $L_{buf1}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_buf1}(n)+b\left(n\right)+{\gamma }_{buf1}\left(n\right)$(3) is similar as the auto-regression used in ARRC [Eq. 2(2) $L_{v\_cr}(n)\text{=}{\displaystyle {\sum }_{t=1}^T{a}_t}\left(n\right)\times L_{t\_cr}(n)+b\left(n\right)$(2) in Cao et al. (2020)], the essential difference is that ARRC establishes regression model for each individual pixel, while VICR build one model for the whole image patch within the cloud region. Furthermore, we innovatively proposed to build two buffer zones around the cloud region, with each playing different roles. Buf1 helps to establish linear relationship between time-series reference images and the target image, and buf2 helps to find the optimal set of reference images. Results in Experiment I verified that this strategy is reasonable, and the results are reliable. More importantly, it eliminates the need for human intervention of reference image selection, which makes the algorithm more operational.

5.2 Advantages of TWSD and its potential applications

The second improvement of VICR is that the newly proposed TWSD can reduce the reconstruction errors by incorporating temporal information in similar pixel selection and residual allocation compared to the widely used SD (as illustrated in Section 4.2). The major contribution of TWSD to the accuracy improvement is in the process of similar pixel selection. Most existing similar pixel selection methods assumed that the similar pixels on the reference image remain similar on the cloud image (Chen et al. 2011), which is probably not true when land cover experienced abrupt change. For those methods based on multi-temporal or time-series reference images (such as ARRC), similar pixels were defined as those with similar temporal trend of spectral reflectance (Chen et al. 2017; Cao et al. 2020; Yan and Roy 2020; Chen et al. 2021). The similarity metrics treats the multi-temporal (or time-series) reference images as equally important. However, in real cases, the greater the differences between the reference images and the cloud image, the less reliability of the selected similar pixels based on these reference images. It should be noted that geostatistical methods can also be used to calculate the weights. Such methods include Geostatistical Neighborhood Similar Pixel Interpolator (GNSPI) (Zhu et al. 2012), which adopted ordinary Kriging to interpolate the residuals. Although GNSPI has shown to perform better than the empirical methods such as NSPI and mNSPI, the calculation of semi-variogram definitely increased the computational burden, especially when using multiple reference images. In this study, we developed TWSD by using the coefficient $a_t$ obtained from linear regression model as the weight in order to consider different contributions from the reference images. A larger $a_t$ value represents that the reference image at the corresponding time t contributes more to the prediction of the cloud image. This indicates that the similar pixels calculated based on this reference image are more reliable. A smaller $a_t$ indicates that the similar pixel selected based on this reference image is less relevant.

TWSD has potential to be applied to other tasks requiring similar pixel selection such as spatial-temporal fusion (STF) algorithms, where spectral change in coarse resolution pixels need to be interpolated according to similarity between neighboring pixels and target pixels to reconstruct fine-resolution image (Zhu et al. 2018; Liu et al. 2019; Zhou et al. 2021). Some STF algorithms utilized multi-temporal Landsat-MODIS image pairs as references. Typical examples include ESTARFM (Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model) (Zhu et al., 2010), STAIR (SaTellite dAta IntegRation) (Luo et al. 2020), and VIPSTF (Virtual image pairs-based spatio-temporal fusion) (Wang et al. 2020). In these algorithms, TWSD has potential be applied to calculate the spectral similarity in temporal domain to better select similar pixels and capture land change.

5.3 Robustness and efficiency of VICR method

Experiment III showed that VICR achieved better accuracies than mNSPI, WLR and ARRC, and the improvement in the accuracies was more obvious in the cases with abrupt land cover change. Compared to the other three algorithms, the performances of VICR were less affected when removing the closest reference image or when increasing the revisiting period (Experiment IV) and were less affected with varying cloud sizes (Experiment V). This is understandable because VICR is not dependent on one reference image, and it flexibly selects the optimal set of time-series reference images. In addition, it captures temporal change in the key processes including virtual image construction, similar pixel selection, and residual distribution by taking into account temporal contribution of each reference image.

The results of Experiment VI showed that VICR achieved good visual effects for real cloud removal and obtained reasonable time-series NDVI. It costs much less processing time than ARRC, and the cumulative processing time of cloud removal for five-year Landsat time-series images is acceptable (Fig. S10). The good computational efficiency can be explained by that the linear regression for virtual image construction is conducted based on images instead of individual pixels. Different from ARRC that performs autoregression for all neighboring pixels and calculates the contribution of all neighboring pixels to the target pixel, the virtual image construction based on image patches can substantially improve the computation efficiency. The combination of patch-based virtual image construction and pixel-based residual distribution enables both good accuracy and efficiency. The promising performance and high computing efficiency highlighted the operationality of VICR for cloud removal from a large volume of Landsat datasets.

5.4 Computation efficiency of VICR method

We expect a good cloud removal algorithm to have reasonable computation efficiency besides high accuracy, so that it can be applied to large number of time-series images in an operational way. In our experiments, we compared the computing time of mNSPI, WLR, ARRC and VICR in Experiment III (the first sub-experiment) and Experiment V. We also calculated the total computing time of VICR for cloud removal from time-series Landsat images in Experiment VI. Although mNSPI is provided in IDL code and ARRC and VICR are written on MATLAB, previous comparisons showed equivalent efficiency between IDL and MATLAB (Liu et al. 2019).

Table 3 lists the computation time of mNSPI, WLR, ARRC, and VICR when performing on the simulated cloud region in Experiment III. Figure 15 illustrates the computation time averaged on each pixel in cloud regions with different sizes in Experiment V. In Experiment III, the computation time is: ARRC > mNSPI > VICR > WLR (Table 3). In most cases in Experiment V, the computation time per pixel of VICR is slightly longer than mNSPI and WLR. In both experiments, ARRC required much longer time than the other three methods. It is interesting that mNSPI took much longer time in Experiment III than in Experiment V even if the cloud sizes are similar (e.g. 93,312 pixels in Jianghan Plain site in Table 3 vs. the 300 × 300 cloud region in Figure 15a). This is probably because mNSPI requires searching for the central pixel of the cloud region, and this process is much less efficient for irregular-shaped cloud than that for regular-shaped cloud. Therefore, it is expected that mNSPI has large variation in computation efficiency for reconstructing real cloud images. The total time of VICR for the largest cloud (160,000 pixels for 400 × 400 cloud) was only 480 seconds (maximum 0.003 seconds per pixel in Figure 15). When processing time-series Landsat images in 5-year periods in Experiment VI, the cumulative processing time was around 15.4 hours (55,000 seconds) for Jianghan Plain (47 images) and 32.8 hours for Yellow River Delta (72 images) (Fig. S9).

Figure 15. Elapsed time per cloud pixel (in seconds) for mNSPI, WLR, ARRC and VICR for cloud size of 50$\times$50, 100×100, 150×150, 200×200, 250×250, 300×300, 350×350 and 400×400 pixels in Experiment V.

Table 3. Elapsed time in seconds for different methods in Experiment III. The cloud pixels are 93312, 107529, 89813, and 103,174 in Jianghan Plain, Lower Gwydir Catchment, Beijing and Yellow River Delta, respectively.

Testing sites	mNSPI	WLR	ARRC	VICR
Jianghan Plain	1892.44	40.89	3464.405	208.97
Lower Gwydir Catchment	1157.27	79.19	5923.87	347.70
Beijing	2407.62	49.14	3305.695	206.21
Yellow River Delta	1615.19	99.49	3883.475	215.64

5.5 Uncertainties and implications for future research

In practical applications of VICR, there are some factors that lead to uncertainties in the cloud removal results. First, VICR still requires reference images to be cloud free within a cloud region. When the cloud region is too large, it may be difficult to select sufficient reference images. This problem can be addressed by time-series cloud removal strategy in Experiment VI. The images with low cloud cover are given higher priority for cloud removal, and these reconstructed images can be used as reference images to remove large cloud. One concern is that the errors produced in the images with small cloud cover may magnify the errors for cloud removal in the image with large cloud. To address this concern, we further conducted an experiment. We simulated a large cloud on the target images in Table 1. The shape of the cloud was taken from the real images in Figure 13 for each site. We removed the cloud following the same procedure in Experiment III using VICR (we call single image reconstruction). We also put the cloud-simulated image into the time-series Landsat images and reconstruct each image in order of the cloud cover percentage (same as time-series cloud removal procedure in Experiment VI). The time-series Landsat images are listed in Table S1 and the reconstructed images are shown in Figure S11. Figure 16 shows that the RMSEs of the reconstructed images from the time-series cloud removal procedure are comparable or even lower than those from the single image reconstruction. This is because the cloud-removed images provide better temporally closer information for reconstruction of the image with larger cloud cover. Nonetheless, both results show acceptable RMSEs and good visual effects (Figure S11), indicating the time-series cloud removal strategy is suitable in the extreme cases where satisfactory reference images are hard to be found for the large cloud region.

Figure 16. Rmses of the reconstructed images from VICR for single image cloud removal procedure and time-series image cloud removal procedure at (a) Jianghan Plain, (b) Lower Gwydir Catchment, (c) Beijing and (d) Yellow River Delta.

Second, in this study we only tested VICR on an identical Landsat sensor. Therefore, the shortest time interval from the reference image to the target image was 16 days. We didn’t incorporate Landsat ETM+ images mainly because the ETM+ sensor had missing strips owing to the failure of the scan-line corrector after May 2003. As Landsat 9 OLI-2 sensor has similar spectral characteristics as Landsat 8 OLI and provides images with 8-day offset with Landsat 8, the performances of VICR for cloud removal have great potential for enhancement when images from both sensors are included. Future work may consider incorporating Harmonized Landsat-8 Sentinel-2 surface reflectance (Claverie et al. 2018) to generate time-series cloud-free images with high temporal frequency.

6. Conclusion

This paper proposed a new thick cloud removal algorithm, i.e. VICR, for Landsat images based on the concept of virtual image. Six experiments were performed at four sites representing different landscape patterns and land change dynamics, and the proposed VICR was compared to existing mNSPI, WLR and ARRC algorithms. The main findings are summarized as follows:

1. VICR achieved better reconstruction accuracies than mNSPI, WLR and ARRC. Results on the cloud-simulated images showed that VICR yielded higher CC, lower RMSE and higher SSIM than mNSPI and WLR (e.g. RMSEs of VICR vs. MNSPI: 0.0288 vs. 0.0405, 0.0286 vs. 0.0382, 0.0193 vs. 0.0227, 0.0234 vs. 0.0317 in the NIR bands for four sites, respectively). VICR yielded higher CC and lower RMSE than ARRC and in most cases yielded higher SSIM than ARRC (e.g. RMSE of VICR vs. ARRC: 0.0193 vs. 0.0211, 0.0286 vs. 0.0302 and SSIM of VICR vs. ARRC: 0.9374 vs. 0.9375, 0.8860 vs. 0.8898 in the NIR bands for Beijing and Lower Gwydir Catchment, respectively in the first sub-experiment in 4.3). Experiments on both cloud-simulated image and real cloud-contaminated images showed that VICR generated visually better reconstruction results than the other algorithms. The accuracy improvement was especially obvious for the sites with abrupt land change, i.e. Lower Gwydir Catchment and Yellow River Delta.

2. VICR is more robust than mNSPI, WLR, and ARRC in terms of the sensitivities to reference image and different cloud sizes. When the time interval between images increased, VICR still achieved the lowest RMSEs (e.g. RMSE: 0.0321, 0.0309, 0.0266 and 0.02496 at NIR band for the four sites in the first sub-experiment in section 4.4). In addition, VICR shows smaller fluctuations in RMSEs with increasing cloud size from 50 × 50 to 400 × 400 pixels.

3. VICR requires much less computation time than ARRC, and less computation time than mNSPI for irregular-shaped clouds (e.g. 32.8 hours for 5-year images in Yellow River Delta).

4. The major advantages of VICR are that it allows automatics determination of optimal set of time-series reference images, and it better incorporates temporal change information. The robustness and efficiency of VICR can promote its operational use in large datasets processing.

Credit author statement

Zhanpeng Wang: Methodology, Investigation, Data processing, Software, Writing-Original Draft. Yinghai Ke: Conceptualization, Methodology, Investigation, Financial acquisition, Writing – Original Draft, Writing-Review & Editing. Demin Zhou: Supervision, Financial acquisition. Xiaojuan Li: Investigation, Project administration. Huili Gong: Investigation, Project administration.

Acknowledgments

The work is supported by the National Natural Science Foundation of China [grant number 42071396].

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement

The data that support the findings of this study are available upon request by contact with the corresponding author, or accessed through https://earthexplorer.usgs.gov.

Supplemental data

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2022.2160411

Additional information

Funding

The work was supported by the National Natural Science Foundation of China [42071396].