Cloud removal using scattering model and evaluation via semi-realistic simulation

ABSTRACT

Cloud removal is an essential task in remote sensing data analysis. As the image sensors are distant from the earth ground, it is likely that part of the area of interests is covered by cloud. Moreover, the atmosphere in between creates a constant haze layer upon the acquired images. To recover the ground image, we propose to use scattering model for temporal sequence of images of any scene in the framework of low rank and sparse models. We further develop its variant, which is much faster and yet more accurate. To measure the performance of different methods objectively, we develop a semi-realistic simulation method to produce cloud cover so that various methods can be quantitatively analysed, which enables detailed study of many aspects of cloud removal algorithms, including verifying the effectiveness of proposed models in comparison with the state-of-the-arts, including deep learning models, and addressing the long standing problem of the determination of regularization parameters. Theoretic analysis on the range of the sparsity regularization parameter is provided and verified numerically.

1. Introduction

As the imaging sensors are deployed kilometres above the earth ground, cloud usually appears in the acquired images, which complicates data analysis tasks. It is desirable to remove the cloud totally to recover clean ground scene, which gives rise to cloud removal. Due to the versatility of remote sensing imagery, cloud removal methods have to align with the characteristics of the sensors, for example, multiple channels or single band. Meanwhile, the platform is a decisive factor for the design of the algorithm, for example, the computation limitation and power consumption restriction. Furthermore, the analysis tasks after cloud removal have some influence as well. So one has to consider all possible contributing factors in the modelling process.

Our data is single band satellites images of the same scene sampled from different time points which are subjected to light to moderate random cloud coverage at various regions. The aim is to recover images without cloud, i.e. the clear images revealing the ground scene so that subsequent analysis can be performed reliably, for example, object detection and tracking. Therefore the fidelity is the most important factor to be considered, in other words, the recovered must be as close as possible to the truth, not just simply ‘visually fit’ (look plausible from afar). Unfortunately, there is no objective assessment except visual checking, and one of the goals of this paper is to fill this gap.

We focus on non-deep-learning-based methods for cloud removal, although latest deep learning methods were used as contenders in our empirical studies subject to code availability (e.g. Zheng, Liu, and Wang 2021; Sarukkai et al. 2020). The reason for this is that the fidelity of the recovered images is a concern for deep learning based methods. The workflow of these methods consists of two steps. The first is to identity cloud covered areas and remove them. The second is to apply generative models to fill the removed pixels. Generalised adversial networks (GAN) based models are popular choice for image completion. However, the working mechanism of GAN and its variants heavily rely on the training data on which the distribution is modelled by transforming a specified random distribution, e.g. uniform distribution or multivariate Gaussian distribution. Essentially, GAN is some sort of density estimator. Then the question is, what if the scene that the satellite sampled never appears in the training data? GAN will certainly generate something for the missing areas but will not be able to stretch outside its modelled distribution. Therefore we consider other alternatives, for example, temporal mosaicing (Guo et al. 2016, 2017). Although enforcing spatial smoothness is the most time-consuming component, the fidelity can be reassured that no ‘alien pixels’ will be inserted into the images like GAN based methods do. Another possibility is matrix completion methods for missing pixel filling, for example Wen et al. (2018), and its later development Zhang et al. (2019). The main model behind these methods is the low-rank robust principal component analysis (Emmanuel et al. 2011) coming from a long development of robust PCA (RPCA) (De la Torre and Black 2001; Junbin Gao, Guo, and Guo 2009) that improves the robustness of the classic linear PCA model by reducing the sensitivity to outliers. The elegance of RPCA comparing to its peers is the simplicity in its formation as well as its theoretical guarantee for the recovery of the low-rank signals and sparse noise. The application of RPCA implies that the observed images are the summation of low rank ground images and sparse cloud cover images (images with cloud only without the ground scene). It makes sense for such arrangement assuming that the ground scene changes little after excluding misalignment and geometric distortion, and clouds cover only small portion of the scene. The low-rank condition on ground component signifies the way of filling missing pixels and hence RPCA has better interpretability than GAN methods.

It seems that the aforementioned two-step workflow can be consolidated to a single one using RPCA. Nonetheless this two-step strategy was still adopted (e.g. in Wen et al. 2018; Zhang et al. 2019), for no obvious reason, in which RPCA is only used for cloud identification and a low-rank matrix completion follows after those cloud affected areas masked out. Two questions remain though. Firstly, where is the atmosphere modelled in the image data?The atmosphere is reflected as a thin haze layer in the acquired images which may not be negligible. Secondly, is the simple additive model in RPCA really the right description of the physics? Apparently not. The most realistic model so far is the so-called atmosphere scattering model (Narasimhan and Nayar 2002) for satellite images. Therefore one should build atmospheric affect into the model for cloud removal and ground scene recover.

2. Models considering atmosphere effects

We first describe RPCA-based methods here in the setting of imagery applications. Let $I_i\in {\mathbb{R}}^{d_1\times d_2}$ be the $i$-th sampled image of size $d_1\times d_2$ and $i=1,\dots ,n$; $D=[vec(I_1),\dots ,vec(I_n)]$ where $vec(X)$ is the vectorization of matrix $X$ to be a column vector, and hence $D\in {\mathbb{R}}^{d\times n}$ ($d=d_1{d}_2$). The RPCA model used in Wen et al. (2018) and Zhang et al. (2019) is the following,

(1) $\underset{L,\mathit{C}}{min}\parallel \mathit{L}{\parallel }_{\ast }+\lambda \parallel \mathit{C}{\parallel }_1$(1)

$\mathrm{s}.\mathrm{t}.\mathit{D}=\mathit{L}+\mathit{C}$

where $\parallel \mathit{X}{\parallel }_{\ast }$ is the nuclear norm of $\mathit{X}$, i.e. the summation of all singular values of $\mathit{X}$, which is the convex envelope for matrix rank, $\parallel \mathit{X}{\parallel }_1$ is the ${\mathrm{\ell }}_1$ norm of $\mathit{X}$, $\mathit{L}$ is the initial recovered ground images, $\mathit{C}$ is the cloud cover images, and both are the same size as $\mathit{D}$. $\lambda$ is the regularization parameter usually fixed to be $\frac{1}{\sqrt{d}}$ as recommended in Emmanuel et al. (2011). By incorporating group sparsity (defined by super-pixels) into (1) Zhang et al. (2019) claim slightly better performance. After solving (2.1), both methods proceed to matrix completion with the mask derived from $C$ as follows

(2) $\underset{\mathit{B},\mathit{S}}{min}\parallel \mathit{L}{\parallel }_{\ast }+\alpha \parallel {\mathit{S}}_{\mathrm{\Omega }}{\parallel }_1+\beta \parallel {\mathit{S}}_{\bar{\mathrm{\Omega }}}{\parallel }_1$(2)

$\mathrm{s}.\mathrm{t}.\mathit{D}=\mathit{B}+\mathit{S}$

where $\mathrm{\Omega }$ is the mask matrix of size $d\times n$ with 0’s for masked out elements and 1’s for others, $\overline{\Omega }$ is the negated version of $\mathrm{\Omega }$, i.e. flipping 0’s and 1’s, and $S_{\mathrm{\Omega }}$ is the projection of $S$ on $\mathrm{\Omega }$, i.e. masking out elements indicated by 0’s in $\mathrm{\Omega }$. The $ij$th element in the mask matrix, $[\mathrm{\Omega }{]}_{ij}=1$ if $[C{]}_{ij}>\gamma \sigma (vec(C))$ and $[\mathrm{\Omega }{]}_{ij}=0$ otherwise, where $\sigma (v)$ is the standard deviation of $v$ and $\gamma \in [0,1]$ is a pre-set ratio. $B$ is the final recovered ground images, which are supposed to be cloud free. $S$ is the noise. In implementation, $\gamma =0.8$, $\alpha =\frac{0.1}{\sqrt{d}}$ and $\beta =1$. Both problems are convex with two blocks of variables. There are many gradient projection based solvers/optimizers for them under the ADMM (Alternating Direction Method of Multipliers) framework (Boyd and Vandenberghe 2004). They all work reasonably well for moderate size of images, for example, $d_1=d_2=1024$ and $n=7$.

The critical step is in (1) where cloud cover $C$ is supposed to be separated. Note that the decomposition of the observed data $D=L+C$ reflects the basic model assumption. As mentioned earlier, this departures from the reality by ignoring atmosphere effect. So instead of simple additive model we propose to use atmosphere scattering (Narasimhan and Nayar 2002), $\mathit{D}=\mathit{L}\circ (1-\mathit{C})+\mathit{C}$, in the modelling, and hence optimize the following

(3) $\underset{\mathit{L},\mathit{C}}{min}\parallel \mathit{L}{\parallel }_{\ast }+\lambda \parallel \mathit{C}{\parallel }_1$(3)

$\mathrm{s}.\mathrm{t}.\mathit{D}=\mathit{L}\circ (1-\mathit{C})+\mathit{C}$

$[L{]}_{ij}\in [0,1],[C{]}_{ij}\in [0,1]$

$\mathit{X}$ where $\mathit{X}\circ \mathit{Y}$ is the element-wise product of matrix $\mathit{X}$ and $\mathit{Y}$ of the same size. In the above formulation, it is assumed that the pixels in observed images are rescaled to $[0,1]$, which is easily done by dividing the maximum digital number of the sensor, but not the maximum of the observed values. Note that (3) is no longer a convex problem as the equality condition is not affine. It is supposed to be much difficult to solve on itself, let alone the boxed conditions clamping the elements in both $\mathit{L}$ and $\mathit{C}$ within $[0,1]$. Nonetheless, there are still some strategies for the optimization. For example, introducing a dummy variable $\mathit{X}$ to untangle the interaction between $\mathit{L}$ and $\mathit{C}$

(4) $\underset{\mathit{L},\mathit{C}}{min}\parallel \mathit{L}{\parallel }_{\ast }+\lambda \parallel \mathit{C}{\parallel }_1$(4)

$\mathrm{s}.\mathrm{t}.\mathit{D}=\mathit{X}\circ (1-\mathit{C})+\mathit{C}$

$\mathit{L}=\mathit{X}$

$[L{]}_{ij}\in [0,1],[C{]}_{ij}\in [0,1],[X{]}_{ij}\in [0,1]$

and proceed with the normal ADMM. However, we observed that this does not converge well enough to be practically useful. Instead, we employ linearization using primal accelerated proximal gradient method (Kei Pong et al. 2010) for its ease in handling entangled nuclear norm optimization and stability. The Lagrange of (3) with proximity is

(5) $\mathcal{L}=\parallel \mathit{L}{\parallel }_{\ast }+\lambda \parallel \mathit{C}{\parallel }_1+⟨\mathit{Y},\mathit{D}-\mathit{L}\circ (1-\mathit{C})-\mathit{C}⟩$(5)

$+\frac{\mu }{2}\parallel \mathit{D}-\mathit{L}\circ (1-\mathit{C})-\mathit{C}{\parallel }_{\mathrm{F}}^2$

leading to

(6) $\mathcal{L}=\parallel L{\parallel }_{\ast }+\lambda \parallel C{\parallel }_1+\frac{\mu }{2}\parallel D-L^{\circ }(1-C)-C+\frac{Y}{\mu }{\parallel }_{\mathrm{F}}^2$(6)

by ignoring constants, where $\parallel \mathit{X}{\parallel }_{\mathrm{F}}$ is the Frobenius norm of $\mathit{X}$, $\mathit{Y}\in {\mathbb{R}}^{d\times n}$ is Lagrangian parameters for the equality condition and $\mu \ge 0$ is the proximity coefficient. Note that (5) and (6) are the proximal form, and the boxed conditions in (3) are ignored at this stage, which will be handled later by feasibility projection after updating all unknowns. Alternating the minimization w.r.t. $\mathit{L}$ and $\mathit{C}$ is adopted here. Apparently minimizing $\mathcal{L}$ with respect to $\mathit{L}$ is difficult due to the term $\mathit{L}\circ (1-\mathit{C})$ although no much trouble for $\mathit{C}$. The gradients are shown below.

(7) $\frac{\mathrm{\partial }\mathcal{L}}{\mathrm{\partial }\mathit{C}}=\lambda \mathrm{\partial }\parallel \mathit{C}{\parallel }_1-\mu (\mathit{D}-\mathit{L}\circ (1-\mathit{C})-\mathit{C}+\frac{\mathit{Y}}{\mu })\circ (1-\mathit{L})$(7)

(8) $\frac{\mathrm{\partial }\mathcal{L}}{\mathrm{\partial }\mathit{L}}=\mathrm{\partial }\parallel \mathit{L}{\parallel }_{\ast }-\mu (\mathit{D}-\mathit{L}\circ (1-\mathit{C})-\mathit{C}+\frac{\mathit{Y}}{\mu })\circ (1-\mathit{C})$(8)

where $\mathrm{\partial }\parallel \mathit{X}{\parallel }_1$ and $\mathrm{\partial }\parallel \mathit{X}{\parallel }_{\ast }$ are subgradients of ${\mathrm{\ell }}_1$ norm and nuclear norm respectively. The stationary point of (7) gives closed form solution ${\mathit{C}}^{\ast }$

(9) $[{\mathit{C}}^{\ast }{]}_{ij}=\left\{\begin{array}{cc}0,& b_{ij}\in [-\lambda ,\lambda ]\\ \frac{\lambda -|b_{ij}|}{a_{ij}}sign(b_{ij}),& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$(9)

where $a_{ij}=\mu [(\mathit{L}-1)\circ (\mathit{L}-1){]}_{ij}$, $b_{ij}=\mu [(\mathit{D}-\mathit{L}+\frac{\mathit{Y}}{\mu })\circ (\mathit{L}-1){]}_{ij}$ and $sign(x)$ is the sign function of $x$ which takes 1 when $x>0$ and $-1$ when $x<0$. It is a straightforward soft thresholding for ${\mathrm{\ell }}_1$ norm minimization. The only difference is the regularization is not global but local or adaptive as the regularization parameter $\lambda$ is rescaled by each $1/a_{ij}$ as shown in (9). Whereas there is no closed form solution for $\frac{\mathrm{\partial }\mathcal{L}}{\mathrm{\partial }\mathit{L}}=0$ because the singular value thresholding (SVT) (Cai, Candès, and Shen 2010) only works for the following general form

$\tau \mathrm{\partial }\parallel \mathit{X}{\parallel }_{\ast }+\mathit{X}-\mathit{A}=0$

where $\tau >0$ is an arbitrary scaler (normally regularization parameter) and $A$ is a matrix the size as $X$. The solution $X^{\ast }$ to above is $X^{\ast }={\mathcal{S}}_{\tau }(A)$ and ${\mathcal{S}}_{\tau }(A)$ is the so-called SVT operator defined as

(10) ${\mathcal{S}}_{\tau }(\mathit{A})=\mathit{U}\left(\begin{array}{ccc}max({\sigma }_1-\tau ,0)& \hspace{0.17em}& \hspace{0.17em}\\ \hspace{0.17em}& \ddots & \hspace{0.17em}\\ \hspace{0.17em}& \hspace{0.17em}& max({\sigma }_n-\tau ,0)\end{array}\right){\mathit{V}}^{\mathrm{T}}$(10)

where $\mathit{U}$ and $\mathit{V}$ are from SVD of $\mathit{A}$, i.e. $\mathit{A}=\mathit{U}\mathrm{\Sigma }{\mathit{V}}^{\mathrm{T}}$ and $\mathrm{\Sigma }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\sigma }_1,\dots ,{\sigma }_n)$.

To work around it, we linearize the smooth part in (6)

(11) $\frac{\mu }{2}\parallel \mathit{D}-\mathit{L}\circ (1-\mathit{C})-\mathit{C}+\frac{\mathit{Y}}{\mu }{\parallel }_{\mathrm{F}}^2\equiv f_s(\mathit{L})$(11)

by the first-order Taylor expansion with a proximal term w.r.t ${\mathit{L}}^k$, the $k$th value of $\mathit{L}$ in the iterative optimization for (6), and optimize $L$ while holding other variables constant as

$\underset{\mathit{L}}{min}\parallel \mathit{L}{\parallel }_{\ast }+⟨\frac{\mathrm{\partial }{f}_s}{\mathrm{\partial }\mathit{L}}{|}_{{\mathit{L}}^k},\mathit{L}-{\mathit{L}}^k⟩+\frac{{\mathrm{\ell }}_p}{2}\parallel \mathit{L}-{\mathit{L}}^k{\parallel }_{\mathrm{F}}^2.$

In above, $\frac{\mathrm{\partial }{f}_s}{\mathrm{\partial }\mathit{L}}{|}_{{\mathit{L}}^k}=\mu (\mathit{D}-{\mathit{L}}^k\circ (1-\mathit{C})-\mathit{C}+\frac{Y}{\mu })\circ (1-\mathit{C})$. ${\mathrm{\ell }}_p$ is the Lipschitz constant of $f_s(L)$, which is the operator norm of $⟨\frac{\mathrm{\partial }{f}_s}{\mathrm{\partial }L}{|}_{L^k},\cdot ⟩$ that maps a matrix of the same size of $L$ to $\mathbb{R}$

$⟨\frac{\mathrm{\partial }{f}_s}{\mathrm{\partial }\mathit{L}}{|}_{{\mathit{L}}^k},\cdot ⟩:{\mathbb{R}}^{n\times d}\to \mathbb{R}$

$\mathit{A}\mapsto ⟨\frac{\mathrm{\partial }{f}_s}{\mathrm{\partial }\mathit{L}}{|}_{{\mathit{L}}^k},\mathit{A}⟩$

It is straightforward to see that ${\mathrm{\ell }}_p=1$ due to the box conditions of $C$ and $L$. This leads to

(12) $\underset{\mathit{L}}{min}\parallel \mathit{L}{\parallel }_{\ast }+\frac{{\mathrm{\ell }}_p}{2}\parallel \mathit{L}-{\mathit{L}}^k+\frac{\mathrm{\partial }{f}_s}{\mathrm{\partial }\mathit{L}}{|}_{L^k}/{\mathrm{\ell }}_p{\parallel }_{\mathrm{F}}^2$(12)

The above linear approximation results is very convenient as the interaction between $L$ and $C$ has been removed and therefore (12) has closed form solution using SVT. We apply Nesterov acceleration to speed up the process, which is proven to be convergent for (12) with carefully chosen optimization parameters (Nesterov 2003). This iterative procedure for $L$ has to be embedded into the optimization for (6) and hence there are two loops in the entire algorithm. The detailed optimization algorithm for solving (3) is listed in Alg. 1. Note that the boxed conditions are satisfied by clamping in Alg. 1, which is the feasibility projection commonly used in many implementations (Jun Liu and Jieping 2009). We call the model in (3) and its realization in Alg. 1 atmosphere cloud removal model and ATM for short.

Algorithm 1 Solving atmosphere scattering model for cloud removal in (3)Require: $\mathit{D}$, $\lambda$, $\text{isin}$1: $\mathit{Y}=\mathit{D}/\lambda$, $\mathit{L}=0$, $\mathbf{C}=0$2: $\rho =1.5$, $\mu =1.25/\parallel \mathit{D}\parallel$, ${\mu }_{max}=\mu \ast {10}^7$3: _while $\frac{\parallel \mathit{D}-\mathit{C}-(1-\mathit{C})\circ \mathit{L}{\parallel }_{\mathrm{F}}}{\parallel \mathit{D}{\parallel }_{\mathrm{F}}}>\text{isin}$ do4: update $\mathit{C}$ using (2.9)5: $C(C<0)=0$, $\mathit{C}(\mathit{C}>1)=1$ ▷ Feasibility projection for $\mathit{C}$6: $\theta =\underset{\_}{\theta }=1$, $\underset{\_}{L}=L$, ${\mathrm{\ell }}_p=1$7: $\underset{\_}{f}=0$, $f=\mathrm{\infty }$, $k_{max}=100$8: while $k<k_{max}$ and $|f-\underset{\_}{f}|>{10}^{-3}$ do9: $W=L+(\theta /\underset{\_}{\theta }-\theta )(L-\underset{\_}{L})$10: $\underset{\_}{L}=L$, $\mathit{L}={\mathcal{S}}_{\frac{1}{\mu {\mathrm{\ell }}_p}}(\mathit{W})$11: $\underset{\_}{\theta }=\theta$, $\theta =\frac{\sqrt{{\theta }^4+4{\theta }^2}-{\theta }^2}{2}$12: $\underset{\_}{f}=f$, $f=\frac{\mu }{2}\parallel \mathit{D}-\mathit{C}-(1-\mathit{C})\circ \mathit{L}{\parallel }_{\mathrm{F}}^2+\parallel L{\parallel }_{\ast }$13: end while14: $\mathit{L}(\mathit{L}<0)=0$, $\mathit{L}(\mathit{L}>1)=1$ ▷ Feasibility projection for $\mathit{L}$15: $Y=vY+D-C-(1-C)°L$16: $\mu =min(\mu \rho ,{\mu }_{max})$17: end while

Due to the double iterative procedure for solving $\mathit{L}$ in ATM, it is expected to be slow. However, the recovered cloud component, i.e. $\mathit{C}$ is closer to reality than that obtained from RPCA as shown in Figure 1, where the source images are from GaoFen4 satellite captured at the same scene at seven time points. GF-4 is the first Chinese geosynchronous orbit remote sensing satellite and equipped with one array camera with resolution of 50 m for the panchromatic band and 400 m for the MWIR band. The test panchromatic data were captured on 20 August 2016, and the geographic coordinates are ${18}^{\circ }{09}^{\prime }{34}^{{ }^{\prime \prime }}-{18}^{\circ }{37}^{\prime }{27}^{{ }^{\prime \prime }}$ north latitude and ${108}^{\circ }{56}^{\prime }{30}^{{ }^{\prime \prime }}-{109}^{\circ }{48}^{\prime }{28}^{{ }^{\prime \prime }}$ east longitude (Hainan Province in China). Since the width of a single-field GF-4 satellite image is 10,240 × 10240 pixels, we cropped the west part of the Hainan Province with lots of clouds as the test data.

Figure 1. Cloud detection for GaoFen4 images of the same scene at seven time points. (a) original GF4 images, (b) ATM detected clouds, (c) RPCA detected clouds. Note that all the images are in vector form, including the ones in later sections. Please zoom in to see details in the PDF file of this paper.

It is clear that the ATM detected clouds are much brighter than those detected by RPCA thanks to its detailed atmosphere model, at the cost of much higher computational load as shown in Figure 10. This motivates us to reduce its computational cost while maintaining model capacity. The key is to disentangle the interaction between $L$ and $C$ that breaks the convexity. Let us take a closer look at the core in ATM model in (2.3), i.e. $D=L\circ (1-C)+C$. We decompose $L$ as $L=(1-P)\circ L+P\circ L$ for $P\in [0,1]$ and $PC$, where $\ge$ means element-wise $\ge$, i.e. $[P{]}_i\ge [C{]}_i$. We proceed using this decomposition

$D=L+C-L\circ C=(1-P)\circ L+C+(P-C)\circ L.$

Under the choice of $P$, $N\equiv (P-C)\circ L0$. The above can be written as

(13) $D=\tilde{\mathit{L}}+C+N$(13)

for $\tilde{\mathit{L}}=(1-P)\circ L$. We can easily write out an equivalent optimization problem to (3) using (13) with many coupling conditions, which complicate the optimization. However, if we drop some coupling conditions via relaxation and approximation, it will be much easier to solve, and yet the coupled problem is still a special case of the relaxed version. So we optimize the following

(14) $\underset{L,C,N}{min}\parallel L{\parallel }_{\ast }+\lambda \parallel C{\parallel }_1+\beta \parallel N{\parallel }_{\mathrm{F}}^2$(14)

$\mathrm{s}.\mathrm{t}.D=L+C+N$

$[L{]}_{ij}\in [0,1],[C{]}_{ij}\in [0,1],[N{]}_{ij}\in [0,1]$

Note that in above $L$ replace $\tilde{\mathit{L}}$ which is an approximation. We highlight this is a relaxed version of (3) with its own interpretation, that is $N$ acts as a thin haze layer accounting for the atmosphere.

In (14) the values in $N$ are controlled by the Frobenius norm. It is well known that the Frobenius norm will not encourage sparsity, but compress the values towards zeros uniformly. Depending on the value of $\beta$, $C+N$ can reach the so-called $\alpha$-sparsity (Yonina and Kutyniok 2012), i.e. sparsity beyond value $\alpha$. Note that we fix $\beta =1$ throughout this paper.

Equations (14)(14) $\underset{L,C,N}{min}\parallel L{\parallel }_{\ast }+\lambda \parallel C{\parallel }_1+\beta \parallel N{\parallel }_{\mathrm{F}}^2$(14) is significantly easier to solve than Equation (3)(3) $\underset{\mathit{L},\mathit{C}}{min}\parallel \mathit{L}{\parallel }_{\ast }+\lambda \parallel \mathit{C}{\parallel }_1$(3) for being convex with no interaction terms in the low-rank component. Although direct generalization of ADMM to more than two blocks of variables like those in (14) may not converge as discussed in (Chen et al. 2014) with hand-crafted counter examples, from many other applications, and vast amount of experiments we carried out, the optimization converged quite quickly. Detailed optimization algorithm is listed in Alg. 2. We call the model in (14) and its optimization algorithm in Alg. 2 alternative ATM, or aATM for short.

Algorithm 2 Solving (2.14)

Require: $\mathit{D}$, $\lambda$, $\beta$, $\text{isin}$

$\mathit{Y}=\mathit{D}/\lambda$, $\mathit{L}=0$, $\mathit{C}=0$, $\mathit{N}=0$

$\rho =1.5$, $\mu =1.25/\parallel \mathit{D}\parallel$, ${\mu }_{max}=\mu \ast {10}^7$

_while $\frac{\parallel \mathit{D}-\mathit{C}-\mathit{L}-\mathit{N}{\parallel }_{\mathrm{F}}}{\parallel \mathit{D}{\parallel }_{\mathrm{F}}}>\text{isin}$ do$[C{]}_{ij}=sign([C{]}_{ij})max(|[C{]}_{ij}|-\frac{\lambda }{\mu },0)$

$\mathit{C}(\mathit{C}<0)=0$, $\mathit{C}(\mathit{C}>1)=1$ ▷ Feasibility projection for $\mathit{C}$$\mathit{L}={\mathcal{S}}_{\frac{1}{\mu }}(\mathit{D}-\mathit{C}-\mathit{N}+\frac{\mathit{Y}}{\mu })$

$\mathit{L}(\mathit{L}<0)=0$, $\mathit{L}(\mathit{L}>1)=1$▷ Feasibility projection for $\mathit{L}$$\mathit{N}=\frac{\mu }{(\beta +\mu )(\mathit{D}-\mathit{L}-\mathit{C}+\frac{\mathit{Y}}{\mu })}$

$\mathit{N}(\mathit{N}<0)=0$, $\mathit{N}(\mathit{N}>1)=1$▷ Feasibility projection for $\mathit{N}$ but not necessary$Y=Y+D-C-L-N$$\mu =min(\mu \rho ,{\mu }_{max})$end while

For the regularization parameters, we provide theoretic analysis on the range of the main regularization parameter $\lambda$ in Appendix A. The results align with the empirical study outcomes presented in the next section. Furthermore, we will present an empirical equation based on numerical method to determine the best value for $\lambda$ as a guidance for practical use. Figure 2 shows the flow chart of our method. It is indeed very straightforward. The input stack of images is first normalized by dividing the maximum digital number of the instrument. Then determine regularization parameter $\lambda$ by the procedure provided in 3.4. Finally run Algorithm in 2 to obtain the clean images and cloud cover. We need to point out that all our models can be used directly to recover ground images unlike the main contenders (Wen et al. 2018; Zhang et al. 2019) where a matrix completion (MC) step has to follow although it is debatable whether MC is necessary. However, without some sort of ground truth, it would be a myth and the arguments would be meaningless. To address this long standing issue, we design a semi-realistic simulation of cloud covered images so that all methods can be analysed quantitatively.

Figure 2. Flowchart of the proposed method. Note that $L{\oplus }_{}C\equiv D=L\circ (1-C)+C$.

3. Quantification of performance

3.1. Simulation and performance indicator

Cloud removal experiments are normally conducted on real images from satellites and the evaluation of the effectiveness of the recovery is based on visual checking and cloud cover by IoU (Intersection over Union) originated from computer vision (Liggins et al. 1997) which is basically Jaccard index (Jaccard 1912). The ground truth of cloud cover is obtained by time-consuming manual labelling of cloud clusters. Due to the complexity of the nature of cloud, it is extremely difficult to delineate the boundary of cloud clusters accurately, especially for thin clouds, and hence there exists large amount of errors when segmenting clouds manually. An ideal solution is to build cloud model to capture the shape and formation of all sorts of cloud, thick or thin. Unfortunately it is quite involved in physics and mathematics and it is a multi-facet problem (Dobashi et al. 1999, 2017; Yuan et al. 2014; Yumeihui Xing et al. 2017). Even if the cloud cover is known, the other side of the problem, way more important than cloud, is the ground truth of the ground scene. The ultimate goal of cloud removal is to recover ground scene accurately. Whereas current practice largely relies on subjective evaluation, i.e. ‘eye-balling’, which is apparently very vulnerable to bias. Therefore, an objective and robust evaluation is highly desirable. The work in Zheng, Liu, and Wang (2021) used an overly simplified method to train the Unet for cloud separation by simulating random strips of white rectangles or from brightest to darkest colour gradient boxes on top of clear ground images. This is a bit primitive. Not only are they far from real clouds, but most importantly the regular shape reduces the complexity of the problem. Inspired by the success of applying Perlin noise (Perlin 1985, 2002) in the simulation of virtual landscapes, we adopt Perlin noise to generate synthetic cloud. We take a cloud free image, say from the Inria aerial image labelling dataset (Maggiori et al. 0000), convert it to greyscale as true ground image $I$ (pixels rescaled to $[0,1]$), and generate multiple 2D Perlin noise the same size as the image, as $C_i$, $i=1,\dots ,n$. Then the observed image $I_i$ is

$I_i=C_i+(1-C_i)\circ \mathit{I}$

where pixels in $C_i$ are rescaled within $[0,1]$. Optionally one can apply any transformation $f$ to $\mathit{I}$ before combining to clouds, e.g. geometric distortion to study some aspects of the methods; or generate a base $C$ and apply dynamics to $C$ for cloud time series mimicking clouds movement. We leave these for future work. By varying the parameters in Perlin noise generator, we can control the density of the generated clouds, from sparse to heavy cover. We also apply some image correction, e.g. Gamma correction and histogram equalization, totally optional, to enhance the similarity to real clouds and haze.

As the ground truth is readily accessible, we can apply any suitable quantitative evaluation to the cloud removal methods for detailed study. Given the main focus is the fidelity of the recovered image, we use the following to quantify the goodness of recovery

(15) $r=\frac{\parallel {\hat{I}}_i-\mathit{I}{\parallel }_{\mathrm{F}}}{\parallel \mathit{I}{\parallel }_{\mathrm{F}}}$(15)

Where ${\hat{I}}_i$ is the recovered image from any method. The quantity $r$ defined in (15) is the normalized distance metric, which is not meant to be the best. Other sophisticated measures could be applied certainly. However, (15) is sufficient by virtual of equivalency of norms (Åžuhubi 2010, Ch6.6), although in modelling process, different norms affect model behaviours vastly.

3.2. Performance evaluation on simulations on single image

Thanks to the above semi-realistic simulation, we can now investigate another important aspect, that is the regularization parameters used in the models. Using the goodness of recovery $r$, we can determine the best values from large-scale randomized trials. Meanwhile we can also verify the necessity of the MC step.

Let us first visually check the outcomes of different methods on one set of simulated images. The true image is from Inria dataset named tyrol-w1 from Lienz in Austrian Tyrol resized to $1024\times 1024$ ($d=2^{20}$). It is a mixture of urbane and nature scene with some high intensity areas such as roads and roof tops shown in Figure 3. We simulate seven thin cloud covers. One simulated image and the cloud layer are also shown in Figure 3.

Figure 3. Simulated image. Left to right: true clear image, one of the simulated image, its cloud cover.

The cloud looks very nature. Note that the cloud cover image appears to be sparse as large dark areas exist as shown in the histograms in Figure 4 top panel where the right one is showing details in the range of $[0,0.2]$. However, they are not exactly zero and correspond to thin haze. If one thresholds them to zero, the cloud cover then becomes very artificial visually. The bottom panels in Figure 4 show thresholded results, by 0.1 and 0.2 respectively from left to right. The visible boundaries of clouds are unpleasant and against the intuition due to the lack of the critical smoothness commonly present in natural images with cloud. This also shows the tremendous difficulty to manually separate clouds in real images.

Figure 4. Details of a simulated cloud cover.

Figure 5 shows the recovered images by different methods obtained with the setting of the regularization parameters as $\lambda =\frac{1}{\sqrt{d}}=\frac{1}{1024}$ and $\beta =1$ in aATM. Simple visual checking tells us that aATM and RPCA are better than ATM as ATM results (with and without MC) contain fair amount of cloud pixels. This may be straightforward. However, it is not clear which one is the best. It appears that aATM is slightly better for less ‘washed away’ areas. It is also impossible to identify the effect of MC. These indicate the limit of visual examination. Nonetheless, the $r$ values of these methods are 0.1758,0.3195, 0.1754 in the order of aATM, ATM, RPCA with MC, and 0.1678, 0.3373, 0.1681 without MC. Now it is clear that aATM without MC is the best and MC does not do anything useful to enhance the results.

Figure 5. Recovered images. Top row: simulated images with cloud cover; 2nd-3rd row (in green box): aATM results with/without MC; 4th-5th row (in blue box): ATM results with/without MC; 6th-7th row (in yellow box): RPCA results with/without MC.

Figure 6 shows the cloud covers detected by these methods. The clouds separated by ATM are in better contrast, i.e. very bright and very dark although it appears very conservative, that is visually sparser than others. In contrast, aATM and RPCA seem to have more cloud pixels identified. Again, it is impossible to tell the difference between aATM cloud and RPCA cloud by visual examination. Note that for aATM the cloud is the summation of $C$ and $N$.

Figure 6. Cloud covers estimated by various methods. From top: known cloud covers, aATM results, ATM results, RPCA results.

One major benefit of simulation is to validate the sensitivity of the regularization parameter, mainly $\lambda$ in the models. We ran large-scale simulation with $n=7$ and 15, using the same true image. We tested 51 values of $\lambda$ equally spaced in log scale with the recommend value, $\frac{1}{\sqrt{d}}$, in the middle, i.e. from 9.7656e-05 to 0.0098, and for each value of $\lambda$, we ran 50 randomized trials. The results are collected in Figures 7 and 8, where each data point is the mean and $\pm 1$ standard deviation of the $r$ values across all trials for a given $\lambda$ value.

Figure 7. Randomised trials $\left(n=7\right)$. $\lambda$ in log scale.

Figure 8. Randomised trials ($n=15$). $\lambda$ in log scale.

Many things can be read out from the plots. The first is that ATM is not as good as competitors for small $n$, e.g. $n=7$, regardless of the choice of the $\lambda$ values. However, it begins to gain advantage when $n$ is larger. This will be investigated later. The second is that $\lambda$ has roughly three zones: 1) failure zone, where the sparsity is too weak and all methods fail with no recovered images; 2) clamping zone, where the sparsity is overwhelming such that sparse component is wiped out and all methods lose the capacity to identify cloud; 3) Goldilock zone, where the algorithms work reasonably well ($r\le 0.2$ for $\lambda \in [0.0007,0.0012]$), including their bests. Of course, these zones have different boundaries for different methods, and their $r$ values inside these zones have different shapes. For example, RPCA seems to have rather flat $r$ values in its Goldilock zone meaning that its performance varies just a little bit if $\lambda$ is from that zone. There exists a value for $\lambda$ which is better than the default recommended value. This holds true for all methods, interestingly with different margin of being true. For example, for RPCA, the margin is smaller, that is the optimal value of $\lambda$ brings 17.23% reduction of $r$ value on average in $n=7$ case, while that is 42.11% for aATM. Similar observation for $n=15$. When all methods take the default value of $\lambda$, aATM without MC works the best on average, which is 22.84% better than RPCA in expectation sense. The overall best performance of aATM against that of RPCA is 43.06% reduction in $r$ value, down from 0.1625 to 0.0941, which is very significant. This is verified by a one-side t-test with null hypothesis of no $r$ values reduction performed on the trials with the optimal and default $\lambda$ values where significance level $\alpha ={10}^{-5}$. The resulting p-value for null hypothesis is extremely low $1.7557\times {10}^{-29}$ strongly supporting the alternative hypothesis that the reduction is quite significant. A very interesting observation is that ATM without MC comes to the second when $n=15$ in terms of the overall best performance, better than RPCA. Figure 9 reveals the details of the $r$ values of both methods in the trials when holding $\lambda$ value constant, $\lambda =4.6741$e$-04$, the optimal value for both methods. The $r$ values of each method vary during the trials due to the randomness of the simulation. RPCA has higher values of $r$ almost constantly with greater variation than ATM. There is no doubt that ATM outperforms RPCA when $\lambda$ is optimal. The third is that MC does not bring much improvement even acts adversely when $\lambda$ is in the Goldilock zone. This claim is strongly supported by statistical evidence. Table 1 shows the one-side t-tests results performed on the trials of various methods with optimal $\lambda$ values for both $n=7$ and $n=15$ cases. The null hypothesis is that MC brings $r$ value reduction on average, i.e. the mean of $r_{MC}-r_{\stackrel{\sim }{MC}}$ is no greater than 0. $r_{MC}$ and $r_{\stackrel{\sim }{MC}}$ are the $r$ values of a method with and without MC respectively. The significance level $\alpha$ is set as low as ${10}^{-5}$. The p-values are extremely low suggesting that the null hypothesis should be rejected almost surely. The only exception is ATM when $n=7$, which favours the MC to further improve its performance. So clearly the recommendation is to omit MC step in cloud removal in these methods, which is extra computation with little benefit. However, we need to point out here though that there are regularization parameters as well in MC, for which we took the default/recommended values. See previous sections for detail.

Figure 9. $r$values of ATM and RPCA in randomised trials ($n=15$) with their best $\lambda$ values ($\lambda =4.6741$e$-04$ for both). X-axis is the repeat indices from 1 to 50. The coloured horizontal lines are the means of corresponding $r$ values.

Table 1. Null hypothesis $H_0$ is $\overline{r_{MC}-r_{\stackrel{\sim }{MC}}}\le 0$. Significance level $\alpha ={10}^{-5}$ in t-tests.

Methods	$n=7$		$n=15$
	p-value	Confidence interval	p-value	Confidence interval
aATM	1.4973e-58	$[0.0430,\mathrm{\infty })$	1.0251e-74	$[0.0629,\mathrm{\infty })$
ATM	1.0000	$[-0.0234,\mathrm{\infty })$	5.8915e-31	$[0.0155,\mathrm{\infty })$
RPCA	1.7878e-17	$[0.0039,\mathrm{\infty })$	7.7511e-44	$[0.0126,\mathrm{\infty })$

To see the comparison more clearly, we present the mean and standard deviation values of the results in some $\lambda$ range (in the Goldilock zone) in Tables 2 and 3 for clarity. The column in the middle of the tables with double column indicates the values when $\lambda$ is equal to the default value. The column-wise best (minimum among all methods) is highlighted by italic font and overall best is highlighted by bold font. They show clearly that aATM is the best in terms of both expected $r$ value and stability reflected by smaller standard deviations.

Table 2. Mean $r$ values of various methods from all trials for some $\lambda$ values.

$\lambda (\times {10}^{-3})$	0.6162	0.6756	0.7408	0.8123	0.8906	0.9766	1.0708	1.1741	1.2874
aATM+MC	0.2436	0.1392	0.1440	0.1506	0.1596	0.1712	0.1860	0.1990	0.2145
ATM+MC	0.4262	0.2358	0.1867	0.2152	0.2618	0.3118	0.3567	0.3847	0.4054
RPCA+MC	0.2543	0.1714	0.1771	0.1826	0.1896	0.1958	0.2058	0.2111	0.2217
aATM	0.1106	0.0941	0.1033	0.1172	0.1374	0.1625	0.1928	0.2203	0.2504
ATM	0.4063	0.2985	0.2075	0.2304	0.2820	0.3297	0.3687	0.3908	0.4078
RPCA	0.1700	0.1652	0.1729	0.1806	0.1905	0.1996	0.2142	0.2288	0.2531

Table 3. Standard deviations of $r$ values of various methods from all trials for some $\lambda$ values.

$\lambda (\times {10}^{-3})$	0.6162	0.6756	0.7408	0.8123	0.8906	0.9766	1.0708	1.1741	1.2874
aATM+MC	0.0176	0.0082	0.0104	0.0120	0.0144	0.0132	0.0149	0.0177	0.0188
ATM+MC	0.0360	0.0112	0.0191	0.0246	0.0381	0.0369	0.0348	0.0360	0.0329
RPCA+MC	0.0138	0.0127	0.0154	0.0166	0.0193	0.0173	0.0182	0.0204	0.0203
aATM	0.0059	0.0091	0.0126	0.0138	0.0181	0.0171	0.0195	0.0237	0.0255
ATM	0.0175	0.0118	0.0182	0.0278	0.0386	0.0348	0.0318	0.0330	0.0307
RPCA	0.0134	0.0153	0.0183	0.0192	0.0221	0.0198	0.0203	0.0239	0.0253

3.3. Computation costs comparison on simulations using single fixed image

We report the time for computation. Figure 10 shows the time consumed by various methods, with and without MC. Similar to previous plots, the data points in the plot are the means and $\pm 1$ standard deviation of the times (in seconds) across all trials for a given $\lambda$ value. Apparently they vary across simulations.

Figure 10. Time consumed in randomised trials ($n=7$). $\lambda$ in log scale.

It is quite obviously that MC is extra work. Given no extra benefit, the computation for MC should be saved. ATM is pretty difficult to solve indeed, reflected by the skyrocketed computational time compared with others. There are double optimization loops inside its solver. Interestingly, when $\lambda$ is correct, ATM takes the most of time to compute on average. When $\lambda$ is growing from the failure zone to the Goldilock zone, a huge jump of needed computation can be observed, which is statistically significant. As ATM’s performance turns very sharply along $\lambda$ values, its computation cost varies accordingly, peaking at where ATM works the best and jumping down quickly. This is a very interesting observation that may lead to a way of selection of its regularization parameter as well as a hypothesis of required computational cost vs $\lambda$ value. Along with the well known regularization path in sparse models (Hastie et al. 2004; Friedman, Hastie, and Tibshirani 2010; Tibshirani and Taylor 2011), this may be a useful route leading to optimal regularization selection in future. This was never possible previously without simulation. In general, aATM is more expensive to compute than RPCA because of the extra block of variables $N$, doubling the cost almost for all $\lambda$ values. However, the base is quite small. when $\lambda$ is in the Goldilock zone, aATM is doubling RPCA from about 6 seconds to 10 seconds. Therefore it is not dramatic.

Again, we present the mean and standard deviation values of the results in some $\lambda$ range (in the Goldilock zone) in Tables 4 and 5 for clarity. The column in the middle of the tables with double column indicates the values when $\lambda$ equal to the default value. The column-wise best (minimum time among all methods) is highlighted by italic font and overall best is highlighted by bold font. They show clearly that RPCA is the fastest and aATM is about 50% more expensive to run at this range of $\lambda$ values. Considering its superiority in recovery performance, this cost is absolutely worthwhile.

Table 4. Means of time consumed by various methods from all trials for some $\lambda$ values.

$\lambda (\times {10}^{-3})$	0.6162	0.6756	0.7408	0.8123	0.8906	0.9766	1.0708	1.1741	1.2874
aATM	9.84	9.56	9.06	9.18	8.82	8.35	5.77	5.71	6.25
ATM	127.17	130.38	104.43	97.66	90.77	80.95	51.69	46.87	45.81
RPCA	6.45	6.58	6.64	6.70	6.61	6.42	3.75	3.30	3.47
aATM+MC	14.59	14.32	13.85	13.97	13.54	12.73	8.69	8.63	9.49
ATM+MC	131.92	135.02	108.91	102.17	95.06	84.95	54.37	49.49	48.61
RPCA+MC	11.26	11.43	11.42	11.51	11.34	10.94	6.72	6.24	6.69

Table 5. Standard deviations of time consumed by various methods from all trials for some $\lambda$ values.

$\lambda (\times {10}^{-3})$	0.6162	0.6756	0.7408	0.8123	0.8906	0.9766	1.0708	1.1741	1.2874
aATM	0.36	0.40	0.39	0.38	0.42	1.07	0.13	0.07	0.84
ATM	6.75	8.16	5.50	5.19	3.93	10.38	1.28	2.66	7.29
RPCA	0.30	0.32	0.31	0.33	0.31	0.87	0.28	0.04	0.48
aATM+MC	0.48	0.51	0.47	0.44	0.50	1.62	0.13	0.11	1.31
ATM+MC	6.80	8.20	5.52	5.20	3.97	10.86	1.29	2.70	7.62
RPCA+MC	0.48	0.48	0.41	0.46	0.46	1.45	0.30	0.08	0.94

Both aATM and RPCA exhibit the same pattern observed from ATM but less pronounced. When $\lambda$ goes form failure zone to Goldilock zone, there is time cost leap and stabilizes for a while and then some up and downs. Again, the zone changing pattern of time cost is a good indicator of entering the Goldilock zone from failure zone. It is possible to exploit it for finding a better $\lambda$ value than the default one, although it is tricker than ATM where the pattern is very clear.

3.4. Determining the best $\lambda$ value

What is the best value for the regularization parameter $\lambda$? This is an inevitable and yet critical question in practice. It is almost impossible to address it without many assumptions and lengthy theoretic analysis. Please refer to Appendix A for Goldilock zone bounds to see the complexity. However, thanks to simulation, we can fit the data to derive some equation for the best $\lambda$ value. Different from drilling into the computational cost pattern suggested by previous section, we look at the best $\lambda$ values of different methods by stretching $n$ from 2 to 250. The ‘best’ is defined as the $\lambda$ value corresponding to the minimum average $r$ value across trials, which we denote as ${\lambda }^{\ast }$. Figure 11 shows the ratio of ${\lambda }^{\ast }$ of all methods to the suggested default value $\frac{1}{\sqrt{d}}$, i.e. $\sqrt{d}{\lambda }^{\ast }$. As $n$ becomes larger, $\sqrt{d}{\lambda }^{\ast }$ decreases exponentially. We turn this into almost linear by applying $log$ twice to $n$, as shown in Figure 11 right panel. From this data, we fit a linear model and derive the following ${\lambda }^{\ast }$ estimator

Figure 11. $\sqrt{d}{\lambda }^{\ast }$v.S $n=2\sim 250$.

(16) ${\hat{\lambda }}^{\ast }=\frac{-0.5682log(log(n))+1.0747}{\sqrt{d}},\mathrm{\forall }n\ge 2$(16)

The red curves in Figure 11 are the values of $\sqrt{d}{\hat{\lambda }}^{\ast }$ at different scales of $n$.

Remark 3.1. ${\hat{\lambda }}^{\ast }$ is estimated for all methods. However, it is possible to derive the estimator for individual method. ATM may be disadvantaged as the fit is not as good as others. ${\hat{\lambda }}^{\ast }$ should be lower bounded by the minimum value of $\lambda$ in Lemma A.1 so

${\hat{\lambda }}^{\ast }=max\{\frac{-0.5682log(log(n))+1.0747}{\sqrt{d}},\frac{1}{\sqrt{dn}}\}$

prevents ${\lambda }^{\ast }$ from being too small. Finally this ${\hat{\lambda }}^{\ast }$ is empirical and approximate with no model assumption. More sophisticated regression methods are possible.

3.5. Performance evaluation on simulations on multiple images

Now we are ready for more comprehensive tests. One last question is how these evaluations hold across different $n$ (image sequence length) and different scenes? To this end, we picked three other images from Inria data set, chicago1, kitsap1 and vienna1, and ran the same randomized trials with $n=5,7,10,12,15,20$, each with 50 repeats. The ground truth images are displayed in Figure 12. In this experiment, we bring in the state-of-the-art deep learning methods (Sarukkai et al. 2020) (called STGAN+Resnet and STGAN+Unet) and (Zheng, Liu, and Wang 2021) (called UNET and UNET+GAN) for a thorough comparison. STGAN provides two variants using Resnet and Unet backbone networks. UNET separates cloud and ground only and UNET+GAN uses GAN to fill thick cloud covered areas. The training of these deep learning models strictly followed the procedures in their code base repository, and were optimized for best performance as per instructions. For our models and RPCA, $\lambda$ was automatically determined by (16) for different $n$.

Figure 12. Four base images for randomised trials. (a) tyrol-w1,(b) chicago1,(c) kitsap1 and (d) vienna1.

Figure 13 visaulises all $r$ values in one place, where the height of the bars are the means and error bars on top show the standard deviation calculated from multiple trials. The results of the same methods are grouped together with different coloured bars showing the results for different $n$. The overall impression is that deep learning methods are not as good although they have quite stable performance across different $n$ values. They may have some advantages when $n$ is small, say $n=5$. STGAN is better than RPCA in kitsap1 although no match for aATM and ATM when $n>7$. Deep learning methods have large performance variations across different scenes, while others are rather consistent. All sparse models have better $r$ values when $n$ grows larger. This suggests that a strategy to boost performance is to increase the sampling frequency moderately. It makes perfect sense as more images provide more information for the missing pixels covered by cloud, and it is more likely that some areas covered in one image are not covered in another. The rank minimization in aATM/ATM/RPCA is designed to fully utilize this. The unnecessity of MC is once again verified in this test. The add-on value of MC is only observable for ATM when $n$ is small, i.e. $n<15$.

Figure 13. Randomised trials for various base images (in the order shown in Figure 12) and sequence length $n$. The legend shows $n$.

The above observations provide us a clear clue to the questions raised at the beginning of this paper and justify our motivation. Deep learning methods in general have lower fidelity (higher $r$ values). This uncertainty poses many questions for subsequent applications. They may have good performance on some specific scenes, for example, pure nature scene like kitsap1. However, it is not clear how GAN’s distribution transformation works. The large performance variation reveals their problems in dealing with different situations. While our models do not have these issues and are interpretable in terms of their working mechanism.

4. Conclusions

In this paper, we introduced atmosphere scattering model into cloud removal modelling process and proposed two ATM models as superior alternatives to RPCA-based model. Furthermore, we proposed a method to simulate controllable cloud cover scenes. This semi-realistic simulation enables detailed study of various cloud removal methods and provides valuable insights to several aspects of the algorithms as large-scale randomized trial and quantitative analysis become available. Examining the methods by using this powerful experimental tool, we saw clearly that the proposed aATM outperforms not only RPCA model, the state-of-the-art in this category of non-deep-learning-based cloud removal methods, but also latest deep learning models constructed on large scale backbone networks, by quite a large margin. There were many interesting findings in this process, for example the zoning of the regularization parameter $\lambda$, computational cost pattern across zones, automated regularization parameter determination and so on. These may be out of the question without the assistance of the simulation. We envisage a robust development of the cloud removal algorithms under this framework in near future. For this purpose, we released the Matlab code at https://github.com/yguo-wsu/clouldremoval.

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix

A Analysis of regularisation parameter $\lambda$

In this section, we focus on the theoretic analysis on the regularisation parameter $\lambda$ in the models, in particular its valid range. We first have the following minimum $\lambda$ value lemma.

Lemma A.1.

For any given data in $\mathit{D}\in {\mathbb{R}}^{d\times n}$ assuming $d>n$, the minimum and maximum value for $\lambda$ in ATM, aATM and RPCA model is $\frac{1}{\sqrt{dn}}$ and $\sqrt{n}$ respectively. The extremum is in the sense of bound for the models to generate non-trivial solutions.

Proof.

The optimality condition of both RPCA and ATM models requires

$\mathrm{\partial }\parallel {\mathit{L}}^{\ast }{\parallel }_{\ast }+\lambda \mathrm{\partial }\parallel {\mathit{C}}^{\ast }{\parallel }_1=0$

leading to

(A.1) $\mathit{U}{\mathit{V}}^{\mathrm{T}}+\mathit{W}+\lambda sign({\mathit{C}}^{\ast })=0$(A.1)

where superscribed star ${ }^{\ast }$ stands for the optimal value, ${\mathit{L}}^{\ast }=\mathit{U}\mathrm{\Sigma }{\mathit{V}}^{\mathrm{T}}$ is the skinny SVD of ${\mathit{L}}^{\ast }$ (i.e. $U\in {\mathbb{R}}^{d\times l}=[u_1,\dots ,u_l]$, $\mathit{V}\in {\mathbb{R}}^{n\times l}=[v_1,\dots ,v_l]$, and $\mathrm{\Sigma }=I_l$, identity matrix of size $l\times l$) and $\mathit{W}\in {\mathbb{R}}^{d\times n}$ is any matrix satisfying

$\parallel \mathit{W}{\parallel }_2\le 1,{\mathit{U}}^{\mathrm{T}}\mathit{W}=0,\mathit{W}\mathit{V}=0.$

It is easy to see that $\mathit{W}={\mathit{U}}_{\mathrm{\perp }}{\mathit{E}\mathit{V}}_{\mathrm{\perp }}^{\mathrm{T}}$, where ${\mathit{A}}_{\mathrm{\perp }}$ is the complementary components in ambient space that is orthogonal to matrix $\mathit{A}$, and $E$ is a diagonal matrix with all element $e_j\in [0,1]$ to satisfy $\parallel \mathit{W}{\parallel }_2\le 1$. With the re-wrting Eq. (A.1), we are seeking

(A.2) $\underset{\lambda }{inf}\{\lambda \ge 0:\sum _{i=1}^l{\mathit{H}}_i+\sum _{j>l}^n{e}_j{\mathit{H}}_j=\lambda sign(-{\mathit{C}}^{\ast })\}$(A.2)

where as $H_i=u_i{v}_i^{\mathrm{T}}$. $H_i$’s are orthogonal to each other and unitary in terms of Frobenius norm. Equation (A.2)(A.2) $\underset{\lambda }{inf}\{\lambda \ge 0:\sum _{i=1}^l{\mathit{H}}_i+\sum _{j>l}^n{e}_j{\mathit{H}}_j=\lambda sign(-{\mathit{C}}^{\ast })\}$(A.2) shows that the subgradient of $\parallel L\parallel$ at $L^{\ast }$ is clamped by $\lambda$ regardless $C^{\ast }$, meaning the elements in the left hand side of (A.2) have to be in the range of $[-\lambda ,\lambda ]$, a boxed condition. The largest norm within the $\lambda$ box is at one of its corners. Without loss of generality, we can choose the first orthant corner. According to Pethagorean, adding an orthogonal component to a vector will only increase the norm. Hence to allow as large as possible for $\lambda$, one can seek the vector with smallest norm, which reflected to the situation in (A.2) is to let $l=1$ and set $\mathrm{\forall }j{e}_j=0$. It is equivalent to choose $U$ and $V$ to be vectors of all $1/\sqrt{d}$ and all $1/\sqrt{n}$ respectively and let $W=0$. In this case, $U{V}^{\mathrm{T}}+W=\frac{1}{\sqrt{dn}}\mathbf{1}$ where $\mathbf{1}$ is matrix of all one’s with compatible dimensions. This is the smallest norm the subgradient of $L$ can fit in the $\lambda$ box. Therefore, the infimum in (A.2), i.e. what is required in this lemma is $\frac{1}{\sqrt{dn}}$.

Similarly the maximum is

(A.3) $\underset{\lambda }{sup}\{\lambda \ge 0:\sum _{i=1}^l{\mathit{H}}_i+\sum _{j>l}^n{e}_j{\mathit{H}}_j=\lambda sign(-{\mathit{C}}^{\ast })\}.$(A.3)

The only difference is that one has to consider all possibilities, i.e. the maximum of the norm. Therefore, (A.3) is equivalent to

$max\{\parallel \sum _{i=1}^l{\mathit{H}}_i+\sum _{j>l}^n{e}_j{\mathit{H}}_j{\parallel }_F\}=\sqrt{\mathrm{t}\mathrm{r}[{\mathit{V}}_n{\mathit{U}}_n^{\mathrm{T}}{\mathit{U}}_n{\mathit{V}}_n^{\mathrm{T}}]}=\sqrt{n}$

where ${\mathit{U}}_n=[u_1,\dots ,u_n]$ and ${\mathit{V}}_n$ likewise.

Remark A.2. From above, we can see that, when $\lambda <\frac{1}{\sqrt{dn}}$, the only allowed solution is to nullify the elements in $\mathit{L}$, in which case, ${\mathit{L}}^{\ast }=0$ and ${\mathit{C}}^{\ast }=\mathit{D}$. This is what we have seen in Fig. 7, where when $\lambda$ is very small, $r$ value is 1 as $L=0$. aATM has the same result although it has another regularisation because the infimum of $\lambda$ happens only when $\beta =\mathrm{\infty }$, otherwise it would further reduce the value of $\lambda$.

This matches the purpose of regularisation. When $\lambda$ is too small, the penalty to sparsity is next to null. Hence the sparse component is free. The sensible choice is of course to set the low rank component zeros, such that the objective is quite small, although this is a trivial solution. Similar logic for maximum value of $\lambda$.

Note that the extrema values of $\lambda$ deduced in Lemma A.1 is for general cases, in other words, no specific conditions. The minimum value of $\lambda$ is very close to the recommended value of $\lambda$ in RPCA, while the maximum is rather loose, due to the generality. Actually we can have the following tighter upper bound of $\lambda$.

Lemma A.3.

For a given $\mathit{D}\in {\mathbb{R}}^{d\times n}$ assuming $d>n$, the maximum value for $\lambda$ in ATM, aATM with $\beta =\mathrm{\infty }$ and RPCA model, written as ${\lambda }_m$, is $\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$ where $\mathit{U}$ and $\mathit{V}$ is from the skinny SVD of $\mathit{D}=\mathit{U}\mathrm{\Sigma }{\mathit{V}}^{\mathrm{T}}$ and $\parallel \mathit{A}{\parallel }_{\mathrm{\infty }}$ is the matrix infinity norm, i.e. the maximum absolution value of all its elements.

Proof.

The proof is similar to that of Lemma A.1 except now we consider only one data set $D$, i.e.

(A.4) $\underset{\lambda }{sup}\{\lambda \ge 0:\sum _{i=1}^l{\mathit{H}}_i\in [-\lambda ,\lambda ]\},$(A.4)

as the sparse component ${\mathit{C}}^{\ast }$ is totally zero. ${\sum }_{i=1}^l{\mathit{H}}_i=\mathit{U}{\mathit{V}}^{\mathrm{T}}$ and all its elements are surely less than or equal to $\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$. Therefore when $\lambda =\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$, the condition in (A.4) holds. This is equivalent to setting $\mathrm{\forall }j,e_j=0$ as in (A.3).

The upper bound of $\lambda$ in Lemma A.3 is much better than that in Lemma A.1, especially when $d\gg n$ and $\mathit{D}\in [0,1]$. However it is possible to further quantify $\lambda =\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$ without actual SVD. We give asymptotic results here of the upper bound of $\lambda =\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$ and hence ${\lambda }_m$. To proceed, we need the following proposition to bound $\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$.

Proposition A.4.

For any matrix $\mathit{A}$ of size $d\times n$ ($d>n$)and its skinny SVD as $\mathit{A}=\mathit{U}\mathrm{\Sigma }{\mathit{V}}^{\mathrm{T}}$, the following holds

$\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}\le \frac{\parallel \mathit{A}{\parallel }_{\mathrm{\infty }}}{{\sigma }_n(\mathit{A})}$

where ${\sigma }_i(\mathit{A})$ is the $i$th largest singular value of $A$ and then ${\sigma }_n(A)$ is the smallest singular value of $A$.

Proof.

Following the same way of thinking from previous lemmas, we see that $\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$ can only bound the Frobenius norm of the matrix spanned by the same bases up to $\sqrt{n}$. In other words, for any matrix of size $d\times n$, if $\parallel \mathit{X}{\parallel }_F\ge \sqrt{n}$, then $\parallel \mathit{X}{\parallel }_{\mathrm{\infty }}\ge \parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$. Also $\parallel {\sum }_i{\sigma }_i{u}_i{v}_i^{\mathrm{T}}{\parallel }_F^2={\sum }_i{\sigma }_i^2$. Therefore, if ${\sigma }_i\ge 1$ for all $i$, then $\parallel {\sum }_i{\sigma }_i{u}_i{v}_i^{\mathrm{T}}{\parallel }_F\ge \sqrt{n}$. Combining this observation with the fact that $\tilde{A}=\frac{1}{{\sigma }_n(A)}A=U\Sigma /{\sigma }_n{V}^{\text{T}}$, we obtain the claim in this proposition, as we have $\parallel \tilde{A}{\parallel }_{\mathrm{\infty }}>\parallel \mathit{U}{\mathit{V}}^{\mathrm{T}}{\parallel }_{\mathrm{\infty }}$ since ${\sigma }_1(\tilde{A})\ge {\sigma }_1(\tilde{A})\ge \dots \ge {\sigma }_n(\tilde{A})=1$.

Now we treat the observed image matrix $D$ as a random matrix whose elements are i.i.d from uniform distribution from $[0,1]$ and hence we do not assume any further structure. We are then concerned with the smallest singular value ${\sigma }_n(D)$ of a non-central random matrix, i.e. the mean is non-zero. There is limited results on smallest singular values. The closest one is (Tao and Van 2010), which deals with centralised random matrix. Fortunately, we only need a lower bound on ${\sigma }_n(D)$. We use the following theorem from (Bai and Yin 1993), which also appeared in (Vershynin 2012.

Theorem A.5

(Bai-Yin’s law). Let $A$ be a $d\times n$ random matrix whose entries are independent copies of a random variable with zero mean, unit variance, and finite fourth moment. Suppose that the dimensions $d$ and $n$ grow to infinity while the aspect ratio $n/d$ converges to a constant in $[0,1]$. Then

${\sigma }_n(\mathit{A})=\sqrt{d}-\sqrt{n}+o(\sqrt{n}),{\sigma }_1(\mathit{A})=\sqrt{d}+\sqrt{n}+o(\sqrt{n})$

almost surely.

Note in Bai-Yin’s law there is no assumption on the distribution but centrality. We then write $\mathit{D}=\mathit{A}+0.5\mathit{E}$ where $\mathit{A}$’s elements are i.i.d from uniform distribution from $[-0.5,0.5]$. We use the Courant-Fischer minimax characterisation of singular values to obtain the bound

${\sigma }_k(\mathit{A})=\underset{dim(\mathcal{S})=k}{max}\left\{\underset{\parallel v{\parallel }_2=1,v\in \mathcal{S}}{min}\parallel \mathit{A}v{\parallel }_2\right\}$

where $\mathcal{S}$ is any subspace of ${\mathbb{R}}^n$. This leads to

${\sigma }_n(\mathit{A})=\underset{\parallel v{\parallel }_2=1,v\in {\mathbb{R}}^n}{min}\parallel \mathit{A}v{\parallel }_2$

as $\mathcal{S}$ is just ${\mathbb{R}}^n$. We have the following lemma tailored to non-central uniform distribution.

Lemma A.6.

Let $\mathit{A}$ be a $d\times n$ random matrix whose entries are independent copies of a random variable from uniform distribution with mean $\epsilon$, unit variance. Suppose that the dimensions $d$ and $n$ grow to infinity while the aspect ratio $n/d$ converges to a constant in $[0,1]$. Then

${\sigma }_n(\mathit{A})\ge \sqrt{d}-\sqrt{n}+o(\sqrt{n})$

almost surely.

Proof.

Let $E$ be a all 1 matrix with compatible dimensions, then $\mathit{A}=\tilde{\mathit{A}}+\epsilon \mathit{E}$ where $\tilde{\mathit{A}}$ satisfy conditions in A.5. We have

${\sigma }_n(\mathit{A}{)}^2=\underset{\parallel v{\parallel }_2=1,v\in {\mathbb{R}}^n}{min}\parallel (\tilde{\mathit{A}}+\epsilon \mathit{E})v{\parallel }_2^2$

$\ge \underset{\parallel v_1{\parallel }_2=1,v_1\in {\mathbb{R}}^n}{min}\parallel \tilde{\mathit{A}}{v}_1{\parallel }_2^2+\underset{\parallel v_2{\parallel }_2=1,v_2\in {\mathbb{R}}^n}{min}\parallel \epsilon \mathit{E}{v}_2{\parallel }_2^2$

$+\underset{\parallel v_3{\parallel }_2=1,v_3\in {ℝ}^n}{\min }{v}_3^{\text{T}}{\tilde{A}}^{\text{T}}E{v}_3$

Note we use $v_j$ $j=1,2,3$ to highlight that these minimisations are separated and hence the above holds. ${\tilde{\mathit{A}}}^{\mathrm{T}}\mathit{E}$ in the third terms gives

${\tilde{\mathit{A}}}^{\mathrm{T}}\mathit{E}=(s_1,s_2,\dots ,s_n{)}^{\mathrm{T}}{1}_n^{\mathrm{T}}$

where $s_i={\sum }_j^d{a}_{ij}$, $a_{ij}$ is the $ij$th element in $\mathit{A}$ and $1_n$ is the vector with all 1 with length $n$. As $a_{ij}$’s are from centralised population, under asymptotic condition, $s_i=0$ almost surely and hence the third term vanishes. The first two terms are the square of the smallest singular values of corresponding matrices, i.e. $\tilde{\mathit{A}}$ and $\mathit{E}$. Since ${\sigma }_i(\mathit{E})=0$ for $i\ne 1$, we have the required inequality by using Theorem A.5.

Combining Lemma A and Proposition A.4, we have the following corollary.

Corollary A.7.

Assume the elements in data matrix $D$ of size $d\times n$ is uniformly distributed in $[0,1]$ and its size grows asymptoticly as in Lemma A. Given the conditions in Lemma A.3, the maximum value for $\lambda$ is almost surely

(A.5) $\frac{2\sqrt{3}}{\sqrt{d}-\sqrt{n}+o(\sqrt{n})}$(A.5)

Proof.

It is simply the rescaling result of Lemma A by recognizing the standard uniform distribution has variance $1/\sqrt{12}$ and also $\parallel \mathit{D}{\parallel }_{\mathrm{\infty }}=1$.

Remark A.8. Although ${\lambda }_m$ in Corollary A.7 is asymptotic result, as the images are quite large, say $2^{10}\times 2^{10}$ in our experiments, i.e. $d=2^{20}$, the bound of ${\lambda }_m$ is quite good. In practice, $n$ is relatively small, typically at the order of 10, $n\ll d$. Therefore we can further simplify (A.5) to

${\lambda }_m\approx \frac{2\sqrt{3}}{\sqrt{d}}$

That is what we see from Fig. 7 and 8 that when $\lambda$ is too large, precisely larger than 0.0032 as shown in the figure, the sparse component is erased, i.e. $\mathit{C}=0$. In this case, our theory predicted ${\lambda }_m=\frac{2\sqrt{3}}{\sqrt{2^{20}}}=0.0034$, very close to our observation. When $\mathit{C}=0$, $\mathit{L}=D$ and hence the observed images with clouds. The variations we see from the figures are due to the simulated clouds.

B

Deep learning models training details

We describe how the deep learning models were trained. Basically, we strictly followed the training instructions from the original papers [32] and [24].

B.1 U-Net+GAN

We used the training data from the original paper. We use the cloudless aerial imagery data set from the National Institute for Research in Computer Science and Automation in France (INRIA) as the ground images dataset, and use the greyscale infrared image data set from National Oceanic and Atmospheric Administration (NOAA) as cloud images dataset, respectively. A cloud covered satellite image was obtained by combining clear (cloud free) ground images with cloud only image by atmospheric scattering models as follows:

1. Clear ground images were taken from the INRIA aerial image data set (INRIA), which contains 180 training images. Each image has 5000 $times$ 5000 pixels with a spatial resolution of 0.3 m. As directed in [32], each ground image was cropped into non-overlapping patches with 256 $\times$ 256 pixels. Therefore, each ground image results in 361 non-overlapping patches of size 256$times$256. Finally, a dataset of 64,980 ground images was generated.

2. Following the instructions from [32], we collected cloud images from the NOAA, totally 150 cloud images from 2018 to form a training set. These images are multispectral (RGB + infrared) images with a resolution of 1920 $\times$ 1080 pixels. They were preprocessed into greyscale images. After preprocessing, each cloud image was downsampled and cropped into 16 non-overlapping patches of size 256$\times$256. Finally, a dataset of 900 cloud images is obtained.

3. Cloud and clear images were randomly combined by using atmospheric scattering model to form cloud covered image data with total of 64,980 images for model training.

Finally, model structure and hyperparameters including optimization parameters were set as required in [32].

B.2 STGAN

We used the training data from the original paper [24]. The dataset in the original paper was constructed using publicly available real Sentinel-2 images. When training STGAN, a total of 3130 sets (256$\times$ 256) of image pairs were constructed, and each set of datasets was obtained by pairing one cloudless ground image with three cloud covered images of the same area at different times. The training, verification and test datasets were split according to the ratio of 8:1:1, and there were 2504, 313 and 313 training, verification and test image groups respectively.

Once again, model structure and hyperparameters including optimization parameters were set as required in [24].