Non-local tensor completion for multitemporal remotely sensed images inpainting

IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 1 Non-local tensor completion for multitemporal remotely sensed images inpainting T eng-Y u Ji, Naoto Y okoya, Member , IEEE, Xiao Xiang Zhu, Senior Member , IEEE, and T ing-Zhu Huang Abstract — This is the pre-acceptance version, to read the final version please go to IEEE T ransactions on Geoscience and Remote Sensing on IEEE Xplore. Remotely sensed images may contain some missing areas because of poor weather conditions and sensor failure. Information of those areas may play an important role in the interpretation of multitemporal remotely sensed data. The paper aims at reconstructing the missing inf ormation by a non-local low-rank tensor completion method (NL-LR TC). First, non-local correlations in the spatial domain are taken into account by searching and grouping similar image patches in a large search window . Then low-rankness of the identified 4- order tensor groups is promoted to consider their correlations in spatial, spectral, and temporal domains, while reconstructing the underlying patter ns. Experimental r esults on simulated and real data demonstrate that the proposed method is effective both qualitatively and quantitatively . In addition, the proposed method is computationally efficient compared to other patch based methods such as the recent proposed PM-MTGSR method. Index T erms —Multitemporal remotely sensed images, missing information reconstruction, tensor completion. I . I N T R O D U C T I O N R EMO TEL Y sensed images are important tools for ex- ploring the properties of our living en vironment and hav e been used in many applications, such as hyperspectral unmixing [1]–[7], classification [8]–[14], and target detection [15]–[20]. Howe ver , these applications are largely limited by the missing information that is introduced when acquiring these data by both/either the defecti ve sensor and/or the poor atmospheric conditions. For e xample, three-quarters of the The work of T .-Y . Ji and T .-Z. Huang was supported by NSFC (61772003, 61402082). The work of N. Y oko ya was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI 15K20955 and Alexander von Humboldt Fellowship for postdoctoral researchers. The w ork of X. X. Zhu has recei ved funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innov ation programme (grant agreement No [ERC-2016-StG-714087]), as well as from Helmholtz Association under the framework of the Y oung In vestigators Group SiPEO (VH-NG-1018, www .sipeo.bgu.tum.de). ( Corr esponding authors: Xiao Xiang Zhu and T ing-Zhu Huang. ) T .-Y . Ji and T .-Z. Huang are with the School of Mathematical Sciences, Univ ersity of Electronic Science and T echnology of China, Chengdu 610054, China (e-mail: tengyu ji@126.com; tingzhuhuang@126.com). N. Y oko ya is with the Department of Advanced Interdisciplinary Studies, Univ ersity of T okyo 153-8904, Japan, the Remote Sensing T echnology In- stitute (IMF), German Aerospace Center (DLR) W essling 82234, German y , and Signal Processing in Earth Observation (SiPEO), T echnical Uni versity of Munich (TUM), Munich 80333, Germany (e-mail: yokoya@sal.rcast.u- tokyo.ac.jp). X. X. Zhu is with the Remote Sensing T echnology Institute (IMF), German Aerospace Center (DLR), 82234 W essling, Germany , and also with Signal Processing in Earth Observation (SiPEO), T echnical University of Munich (TUM), 80333 Munich, Germany (e-mail: xiao.zhu@dlr .de). detectors (in band 6) of the Aqua Moderate Resolution Imag- ing Spectroradiometer (MODIS) are inef fecti ve [21], [22], the scan line corrector (SLC) of the Landsat enhanced thematic mapper plus (ETM+) sensor has permanently failed [23], [24], and the ozone monitoring instrument (OMI) onboard the Aura satellite suffers a row anomaly problem. On the other hand, the clouds cover approximately 35% of the Earth’ s surface at any one time [25]. Owing to the abo ve two reasons, missing information is ine vitable in optical remotely sensed images, particularly in the multitemporal image analysis. Thus, reconstructing the missing information is highly desirable. Recently , many reconstruction methods for remotely sensed images have been proposed, which can be classified into four categories: spatial-based, spectral-based, temporal-based, and hybrid methods. The spatial-based methods take adv antage of the relationships between different pixels in the spatial dimension without any other spectral and temporal auxiliary images and include interpolation methods [26], propagated diffusion methods [27], [28], variation-based methods [21], [29]–[33], and ex emplar-based methods [34], [35]. This kind of method cannot reconstruct a large missing area because there is not enough reference information. The spectral-based methods borrow the correlati ve infor- mation from other spectral data to reconstruct the missing area. The basic idea of this kind of method is to estimate the relationship of the kno wn areas between the complete and incomplete bands and then use the relationship to reconstruct the missing area. The typical example of this kind of method is Aqua MODIS band 6 inpainting. F or e xample, W ang et al. [22], Rakwatin et al. [36], and Shen et al. [37] reconstructed the missing area of band 6 by considering the spectral re- lationship with band 7 because these two bands are highly correlated. Gladko va et al. [38] and Li et al. [39] took the relationships between band 6 and the other six bands into consideration. These methods can reconstruct a lar ge area and get a better result than the spatial-based methods. Howe v er , for the most remotely sensed images, all bands contain the same missing areas. For this case, the spectral-based methods fail in getting a promising result. The third class of methods is to reconstruct the missing area by making use of other data taken at the same location and different periods. T emporal-based methods hav e been widely studied for remotely sensed inpainting, especially cloud remov al. The clouds vary with time. Thus, the missing areas in different images are diverse. The basic temporal-based method is to replace the missing area with the same area of different periods [23], [25], [40]. Inspired by the temporal filter methods for the one-dimensional signal denoising, many researchers IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 2 dev eloped temporal filter methods by regarding the temporal fibers as signals [41]–[43]. Recently , temporal learning model- based methods e xploit the compressed sensing and re gression technologies to reconstruct the missing information [44], [45]. More recently , W ang et al. [46] proposed a temporally con- tiguous robust matrix completion model for cloud remov al by making the best use of the temporal correlations: lo w- rankness in time-space and temporal smoothness. As the method (ALM-IPG) considers the local temporal correlation, temporally contiguous property , it is good at processing the data whose adjacent temporal images are similar . The abov e three classes of methods make use of only one kind of relationship (spatial domain, spectral domain, or temporal domain). In some cases, the y are po werful, but sometimes they are not. T o get a better result, the hybrid meth- ods were introduced to extract the complementary information from two or three domains. This kind of method includes the joint spatio-temporal methods [23], [47], joint spatio-spectral methods [48], and joint spectral-temporal methods [49]. Re- cently , Li et al. [50] proposed the patch matching-based multitemporal group sparse representation (PM-MTGSR) that makes use of the local sparsity in the temporal domain and the non-local similarity in the spatial domain to reconstruct the missing information. Obviously , PM-MTGSR is a joint spatio-temporal method, namely , it also makes use of only two domains relationships. The hybrid methods perform better than each of the three classes of methods. This indicates that the results would be better if a method takes advantage of more latent structures information in the observed data. In this paper , we present a ne w methodology that makes full use of spatial, spectral, and temporal relationships for the reconstruction of missing data in multitemporal remotely sensed images. The proposed method is designed to be good at processing not only the temporally contiguous data but also the data that have a large difference between the adjacent temporal images. Lo w-rank tensor-based methods characterize the global correlations for each dimension. Inspired by this, the paper introduces a non- con v ex low-rank approximation for tensor rank to make the best use of the global correlations in the spatial, spectral, and temporal domains. Similar concept has been applied for time series analysis of radar data [51]. T o take advantage of the three domains similarities, we group similar patches and con- sider their low-rankness. Experimental results on both cloud remov al and destriping experiments show that our low-rank approach is capable of achieving more accurate reconstruction than other state-of-the-art approaches. This paper is organized as follows. Some notations are introduced in Section II. Section III describes the proposed algorithm for the multitemporal remotely sensed image re- construction. Section IV presents the experimental results and discussion, and the conclusion is giv en in Section V. I I . P R E L I M I N A RY In this paper , we use non-bold low-case letters for scalars, e.g., x , bold low-case letter for vectors, e.g., x , bold upper- case letters for matrices, e.g., X , and bold calligraphic let- ters for tensors, e.g., X . Moreo ver , we also use bold norm up-case letters for clusters of some variables, e.g., M = ( M 1 , · · · , M N ) . An N -order tensor is defined as X ∈ R J 1 ×···× J N , and x j 1 , ··· ,j N is its ( j 1 , · · · , j N ) -th component. Fibers are the higher -order analogue of matrix rows and columns. A fiber is defined by fixing e very index but one. For example, x : j 2 ··· j N = ( x 1 ,j 2 , ··· ,j N , · · · , x J 1 ,j 2 , ··· ,j N ) is one of mode-1 fibers of N -order tensor X ∈ R J 1 ×···× J N . The mode- n unfolding of a tensor X is denoted as X ( n ) ∈ R J n × Q j 6 = n J j , whose columns are the mode- n fibers of X in the lexico- graphical order . Fig. 1 shows the mode- n ( n = 1 , 2 , 3) fibers and unfoldings for a 3-order tensor . The inv erse operator of unfolding is denoted as “fold”, i.e., X = fold n ( X ( n ) ) . Then n -rank of an N -order tensor X , denoted as rank n ( X ) , is the rank of X n , and the rank of X based on n -rank is defined as an array: rank ( X ) = ( rank 1 ( X ) , · · · , rank N ( X )) . The tensor X is low-rank, if X ( n ) is low-rank for all n . Please refer to [52] for its extensiv e o vervie w . Fig. 1. Mode- n fibers and corresponding unfoldings of a 3-order tensor . I I I . M E T H O D O L O G Y Missing information is inevitable in the observation process for remotely sensed images. The existing methods that charac- terize the correlation are mostly interpolation [26], sparse [50], smooth [46], and low-rank technologies [46], no matter which of the four methods (spatial, spectral, temporal, or hybrid) is adopted. For e xample, PM-MTGSR [50] characterizes the lo- cal relationships in the temporal domain using the group sparse technology , and ALM-IPG [46] characterizes the local and global correlations in the temporal domain using the smooth and low-rank technology , respecti vely . Although these two methods take adv antage of spatial and temporal relationships, they prefer the relationships in the temporal domain. Recently , low-rank tensor based methods have attracted much attention regarding the completion of high-dimensional images because the tensors rank can characterize the correlations in different domains [53], [54]. Combined with the definition of n -rank that is a vector composed of ranks of mode- n unfoldings, it can be seen from the Fig. 1 that the rank of mode- n unfolding describes the correlations of mode- n fibers. T o present the motiv ation in detail, we analyze the low-rankness of some 4- order tensor groups stacked by the 3-order similar patches that are extracted from the 4-temporal cloud-free Landsat- 8 data (“Image 3” in Fig. 7, Section IV); see T ab . I. For example, the first group is of size 4 × 4 × 3 × 708 , where the four dimensions correspond to the numbers of pixels, observations, spectral channels, and patches, and it has 8496 mode-1 fibers, 8496 mode-2 fibers, 11328 mode-3 fibers, and 48 mode-4 fibers. The dimensions of the spaces (DimSpac) generated by mode-1, -2, -3, and -4 fibers are 2, 3, 2, and 13, respectiv ely . That means the mode- n ( n = 1 , · · · , 4) fibers are highly correlated. i.e., it is possible to reconstruct the IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 3 T ABLE I A NA L Y S IS O F L O W - R A N KN E S S O F 4 - OR D E R G RO U P S . Group mode-1 mode-2 mode-3 mode-4 1 # of fibers 8496 8496 11328 48 DimSpac 2 3 2 13 2 # of fibers 360 360 480 48 DimSpac 2 3 2 9 3 # of fibers 5424 5424 7232 48 DimSpac 2 3 2 12 4 # of fibers 1956 1956 2608 48 DimSpac 2 2 2 12 5 # of fibers 3444 3444 4592 48 DimSpac 2 3 2 14 6 # of fibers 11484 11484 15312 48 DimSpac 2 2 2 15 missing area using tensor low-rank technology . Inspired by this, we introduce the tensors rank to characterize the global correlations in the spatial, spectral, and temporal domains to reconstruct the missing information of remotely sensed images after grouping the similar patches. It should be noted that missing areas are detected before their reconstruction. The flowchart of the proposed NL-LR TC is shown in Fig. 2. The proposed NL-LR TC method consists of three parts. The method first reshapes the observed 4-order tensor into a 3-order tensor so that the pixels at the different periods but the same location become adjoining. Next, we search and group similar patches in a searching window . Last, the missing information of ev ery group is reconstructed using the low-rank tensor completion method. A. Data Rearrang ement The observed multitemporal remotely sensed data set Y ∈ R m × n × b × t is a 4-order tensor , where m × n denotes the num- ber of pixels of remotely sensed images, b denotes the number of spectral channels of remote sensors, and t is the number of time series. The indices set Ω ∈ R m × n × b × t is also a 4-order tensor , where the position ( i, j, k, l ) ∈ Z m × Z n × Z b × Z t is covered by cloud if Ω i,j,k,l = 0 and is cloud free if Ω i,j,k,l = 1 . T o make the best use of correlations between different peri- ods and find more accurate similar patches, it is necessary to reshape the observed data. The reshaping process is illustrated in the first step of Fig. 2. PM-MTGSR also reshaped the data before searching the similar patches. The difference between PM-MTGSR and our method is that PM-MTGSR reshapes the 4-order tensor into a matrix, while the result of our reshaping is a 3-order tensor . The dif ference leads to another difference between these two methods: the similar patches in our method are 3-order tensors but matrices in PM-MTGSR. As the description abov e indicates, the aim is to reshape the observed data into a 3-order tensor to take advantage of the temporal correlations. Detailedly , we stack the mode-1 slices at the same locations and different periods one by one, i.e., the observed 4-order tensor Y is reshaped into a 3-order tensor ˆ Y ∈ R m × tn × b which is defined by ˆ Y u,v ,w = Y i,j,k,l when u = i , v = ( j − 1) n + l , and w = k . Similarly , we reshape the indices tensor Ω ∈ R m × n × b × t into a 3-order tensor ˆ Ω ∈ R m × tn × b which is defined by ˆ Ω u,v ,w = Ω i,j,k,l when u = i , v = ( j − 1) n + l , and w = k . Let X be the recov ery data we are seeking and ˆ X ∈ R m × tn × b be the corresponding reshaped 3-order tensor . B. Gr ouping of Similar P atches This section is to search and group the similar patches for missing area pixels of the reshaped data ˆ Y . The second step of Fig. 2 shows the process of similar patch searching after reshaping. The red square denotes the target patch ˆ Y i,j = ˆ Y ( i : ( i + w − 1) , j : ( j + w − 1) , :) ∈ R w × w × b , where ( i, j ) denotes the coordinate of the target patch, and w × w × b denotes the patch size. When the target patch is gi ven, similar patches are searched for in the surrounding window with a radius of r in the reshaped data ˆ Y . According to the reshaping procedure only the similar information in the square window of size (2( r /t ) + 1) × (2( r /t ) + 1) in the original data Y is used. Given a metric of the similarity indicator between the target patch and the other patches and an indicator threshold γ 2 , a similar patch can be found when the indicator reaches the given condition. There are many indicators for measuring similarity between two vectors [50], [55], such as the Euclidean distance, the mean relativ e error , normalized cross-correlation, cosine coefficients, and so on. This work adopts the normalized cross-correlation defined as: Q = P j 1 ··· j N ( x j 1 , ··· ,j N − µ X )( y j 1 , ··· ,j N − µ Y ) r P j 1 ··· j N ( x j 1 , ··· ,j N − µ X ) 2 r P j 1 ··· j N ( y j 1 , ··· ,j N − µ Y ) 2 , where X , Y ∈ R J 1 ×···× J N , µ X , µ Y are the mean v alues of X and Y , respecti vely . The mean v alue of an N -order tensor X is defined as µ X := 1 N P j 1 , ··· ,j N x j 1 , ··· ,j N . After completing a search for similar patches, these 3- order-tensor patches are stacked into a 4-order tensor ˆ Y G ∈ R w × w × b × n , where n is the number of similar patches. The corresponding indices set ˆ Ω G can be obtained according to the coordinates of patches in ˆ Y G . C. Low-rank Reconstruction In this section, we propose a low-rank method to reconstruct the missing information in the 4-order tensor ˆ Y G obtained in the last subsection. Dif ferent from the low-rank matrix methods which consider only one mode correlation, e.g., [46], NL-LR TC studies the low-rankness of ˆ Y G from four aspects, i.e., NL-LR TC considers the correlations of mode- i ( i = 1 , 2 , 3 , 4) using rank ( ˆ Y G ) . In fact, the four dimensions of ˆ Y G denote spatial, temporal, spectral, and patch similarity , respectiv ely . As mentioned in the description about tensor rank previously with Fig. 1, rank ( ˆ Y G ) takes advantage of all of the spatial, spectral, and temporal relationships. This can be found in Fig. 3 where NL-LR TC considers the low-rankness of the group of similar patches by analyzing the lo w-rankness of four unfolding matrices that can describe the correlations of the spatial, temporal, and spectral domains. In contrast, the group of PM-MTGSR only describes the spatial and temporal domains relationships (seen Fig. 3 of [50]). This is another IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 4 Fig. 2. Illustration for the proposed NL-LRTC method. “Height” denotes one of the spatial mode, “W idth” denotes another mode of spatial domain, “W idth*T ime” means this mode contains the information of both spatial (Width) and temporal mode. The proposed method comprises three steps: (1) reshape the 4-order observed data into 3-order data; (2) search and group similar patches; (3) reconstruct each group using the low-rank tensor completion method. difference between NL-LR TC and PM-MTGSR. T ensor nu- clear norm and the corresponding algorithm (HaLR TC) were dev eloped in order to make it possible to minimize the tensor rank [53]. Howe ver , the tensor nuclear norm cannot treat the different singular values accurately according to their different importance. For the proposed non-con ve x surrogate of tensor rank, the lar ger singular values that are associated with the major projection orientations and are more important can be shrunk less to preserve the major data components [56], [57]. This is one of the dif ferences between HaLR TC and NL- LR TC. Another difference is that HaLR TC is without the patch strategy . In the last section, ˆ X G and ˆ Ω G hav e been obtained. Next, we reconstruct the missing areas in ˆ X G group by group using the follo wing model: min ˆ X G logDet ( ˆ X G , ε ) = 4 X i =1 α i L ( ˆ X G ( i ) ) s.t. ˆ X G ˆ Ω G = ˆ Y G ˆ Ω G , (1) where L ( ˆ X G ( i ) ) = log det (( ˆ X G ( i ) ˆ X T G ( i ) ) 1 / 2 + ε i I i ) , I i is the Identity matrix, α i s are constants satisfying α i ≥ 0 and P N i =1 α i = 1 , ε = ( ε 1 , · · · , ε N ) T > 0 , and ˆ X G ( i ) is the mode- i unfolding of ˆ X G . The L ( ˆ X G ( i ) ) can be approximated by using the first-order T aylor e xpansion: L ( ˆ X G ( i ) ) ≈ J i X j =1 σ j ( ˆ X G ( i ) ) σ k j ( ˆ X G ( i ) ) + ε i + constant =( ω k ) T σ + constant M = L ω k ( ˆ X G ( i ) ) , (2) where σ k j ( ˆ X G ( i ) ) s are the solutions obtained in the k -th itera- tion, ω k = (1 / ( σ k 1 ( ˆ X G ( i ) ) + ε i ) , · · · , 1 / ( σ k J i ( ˆ X G ( i ) ) + ε i )) T , σ = ( σ 1 ( ˆ X G ( i ) ) , · · · , σ J i ( ˆ X G ( i ) )) T , and ( J 1 , J 2 , J 3 , J 4 ) = ( w , w , b, n ) . From the approximate function (2), we can see that the proposed function logDet indeed shrinks the lar ger singular v alues less. Next, we present a computationally ef ficient algorithm that is based on the alternating direction method of multipliers (ADMM) [58]–[60] to solve the problem (1) by replacing L ( ˆ X G ( i ) ) with L ω k ( ˆ X G ( i ) ) and introducing some auxiliary values. Thus, the problem (1) can be rewritten as: min ˆ X G , M 1 ,..., M 4 1 ˆ Ω G ˆ Y G ( ˆ X G ) + 4 X i =1 α i L ω k ( M i, ( i ) ) s.t. M 1 = ˆ X G , · · · , M 4 = ˆ X G , (3) where 1 ˆ Ω G ˆ Y G ( ˆ X G ) = 0 if ˆ X G ˆ Ω G = ˆ Y G ˆ Ω G , otherwise 1 ˆ Ω G ˆ Y G ( ˆ X G ) = ∞ . By attaching the Lagrangian multiplier { Λ i } 4 i =1 that hav e the same size with ˆ X G to the linear constraint, the augmented Lagrangian function of (3) is giv en by: L ( ˆ X G , M 1 , . . . , M 4 , Λ 1 , . . . Λ 4 ) = 1 ˆ Ω G ˆ Y G ( ˆ X G )+ 4 X i =1  α i L ω k ( M i, ( i ) ) + β 2 k ˆ X G − M i + Λ i β k 2 F  , (4) where β is the penalty parameter for the violation of the linear constraints, and h· , ·i is the sum of the elements of the Hadamard product. Thus, ˆ X G and { M i } 4 i =1 can be obtained by minimizing the augmented Lagrangian function (4), alternately . The Lagrangian multipliers are updated as Λ k +1 i = Λ k i + β ( ˆ X k +1 G − M k +1 i ) for i = 1 , . . . , 4 . IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 5 Fig. 3. Illustration for how to exploit the lo w-rankness of the group of similar patches using four unfolding matrices. The mode-1 unfolding is of size w × w bn , the mode-2 unfolding is of size w × w bn , the mode-3 unfolding is of size b × w 2 n , and the mode-4 unfolding is of size n × w 2 b . First, ˆ X G is obtained by solving the following optimization subproblem: min ˆ X G ( 1 ˆ Ω G ˆ Y G ( ˆ X G ) + 4 X i =1 β 2 k ˆ X G − M k i + Λ k i β k 2 F ) . (5) It is obvious that the objecti ve function of (5) is dif ferentiable, thus ˆ X k +1 G has a closed form solution: ˆ X k +1 G = 1 4 β 4 X i =1 ( β M k i − Λ k i ) ! ˆ Ω c G + ˆ Y G ˆ Ω G , (6) where ˆ Ω c G is the complementary set of the indices set ˆ Ω G . Next, { M i } 4 i =1 -subproblems are solved. Note that M i - subproblems are independent, and thus we can solv e them separately . W ithout loss of generality , the typical v ariable M i is solved through the follo wing problem: min M i, ( i ) α i β L ω k ( M i, ( i ) ) + 1 2 k M i, ( i ) − ˆ X k +1 G ( i ) − Λ k i, ( i ) β k 2 F , (7) where ω k = (1 / ( σ 1 ( M k i, ( i ) ) + ε i ) , · · · , 1 / ( σ J i ( M k i, ( i ) ) + ε i )) T , ( J 1 , J 2 , J 3 , J 4 ) = ( w , w , b, n ) . M k +1 i, ( i ) can be obtained using a thresholding operator [56], [61], M k +1 i, ( i ) = U k ( Σ k − τ diag ( ω k )) + ( V k ) T , (8) where U k Σ k ( V k ) T is the SVD of ˆ X k +1 G ( i ) + 1 β Λ k i, ( i ) and ( X ) + = max { X , 0 } . Thus, M k +1 i = fold i ( M k +1 i, ( i ) ) . Based on the previous deriv ation, we develop the lo w-rank method to reconstruct missing information in multitemporal remotely sensed images, as outlined in Algorithm 2 . Then the proposed NL-LR TC method is outlined in Algorithm 1 . I V . E X P E R I M E N T S A N D D I S C U S S I O N A. T est Data The proposed reconstruction method, NL-LR TC, for mul- titemporal remotely sensed images is applied to three data Algorithm 1 NL-LR TC for multitemporal remotely sensed images inpainting. Input: Data Y and inde x set Ω , radius of searching window r , patch size w , and indicator threshold γ 2 , parameters β , α , and ε . 1: Obtain the 3D tensors ˆ Y and ˆ Ω by rearranging the data Y and Ω , respectiv ely; 2: while ˆ Ω c 6 = 0 do 3: Find ( i, j ) subject to ˆ Ω i,j = 0 , that means the pixel ( i, j ) is co vered by clouds; 4: Search the similar patches for patch ˆ Y i,j in the search- ing windo w; 5: Stack these similar patches as a group ˆ Y G , and obtain the corresponding inde x set ˆ Ω G ; 6: Estimate the missing pixel v alues in ˆ Y G using Algo- rithm 2 and set ˆ Ω G = 1 ; 7: Replace the corresponding entries in ˆ Ω and ˆ Y with ˆ Ω G and ˆ X G , respecti vely . 8: end while Output: Recov ered data X . Algorithm 2 Lo w-rank reconstruction of missing information via ADMM. Input: Group ˆ Y G and index set ˆ Ω G , parameters β , α , and ε . 1: Initialize: { M i } 4 i =1 = 0 , { Λ i } 4 i =1 = 0 , tol = 10 − 5 , and K = 100 . 2: for k = 1 to K do 3: Update ˆ X k +1 G by (6); 4: for i = 1 to 4 do 5: Update M i by (8); 6: Update Λ k +1 i by: Λ k +1 i = Λ k i + β ( ˆ X k +1 G − M k +1 i ) ; 7: end f or 8: If    ˆ X k +1 G − ˆ X k G    F /    ˆ X k G    F < tol, stop iteration. 9: end f or Output: Recov ered data ˆ X G for group ˆ Y G . IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 6 sets for simulated and real-data experiments. The first data set was taken over Munich, German y , by Landsat-8. The data set has nine bands, and three bands with 30-m res- olution (red, green, blue) are used. The data set was ac- quired over the Munich sub urbs (which consist of forests, mountains, hills, etc.) and includes six temporal images de- noted as “M102014”, “M012015”, “M022015”, “M032015”, “M042015”, and “M082015”, where “MXXYYYY” means the data is taken over Munich in XX-th month YYYY -year; see Fig. 4. The second data set was taken over Beijing, China, by Sentinel-2 with six spectral bands at a ground sampling distance of 20 meters (bands 5, 6, 7, 8A, 11, 12). The data set w as acquired ov er the Beijing suburbs (which consist of villages, mountains, etc.) and includes fiv e temporal images denoted as “BJ122015”, “BJ032016”, “BJ072016”, “BJ082016”, and “BJ092016”, where “BJXXYYYY” means the data is tak en over Beijing in XX-th month YYYY -year; see Fig. 5. The third data set w as taken ov er Eure, France, by Sentinel-2 and atmospheric correction has been processed by MA Y A [62]. The data set includes four temporal images and four spectral bands at a ground sampling distance of 10 meters (bands 2, 3, 4, and 8), see Fig. 6. In our e xperiments, the observed multitemporal remotely sensed data contain four different temporal data, namely , the observed tensor Y is of size m × n × b × 4 . F or the first data set, “M032015”, “M042015”, and “M082015” are the reference data; four subimages of “M012015” and “M022015” are used for the simulated experiments in that the size of the tested images is 512 × 512 in the spatial domain; “M102014” is used for the real experiment in that the size of the tested images is 1080 × 1920 in the spatial domain. For the Munich area, the surface reflectance is changed with time due to sno w , seasonal change of vegetation, etc. Thus, “M012015”, “M022015”, and “M102014” are greatly different from the other reference data. W e also study how NL-LR TC performs when the temporal difference is not so great with the second data set. For this data set, “BJ092016” and “BJ072016” are used for simulated and real experiment, respectively . The other temporal data are as reference data. The size of tested Beijing images is 256 × 256 in the spatial domain. W e test another real experiment with “EU082017” in the third data set whose size is 400 × 400 in the spatial domain. (a) M102014 (b) M012015 (c) M022015 (d) M032015 (e) M042015 (f) M082015 Fig. 4. Data set taken by Landsat-8. “MXXYYYY” means the image is taken over Munich in XX-th month YYYY -year . (a) BJ122015 (b) BJ032016 (c) BJ072016 (d) BJ082016 (e) BJ092016 Fig. 5. Band 6 of Beijing data. “BJXXYYYY” means the image is taken over Beijing in XX-th month YYYY -year . (a) EU102016 (b) EU012017 (c) EU052017 (d) EU082017 Fig. 6. RGB bands (bands 4, 3, and 2) of Eure data. “EUXXYYYY” denotes the image taken over Eure in XX-th month YYYY -year . B. P erformance Evaluation In the simulated experiments, the performance of multitem- poral remotely sensed images reconstruction is quantitativ ely ev aluated by peak signal-to-noise ratio (PSNR) [21], structural similarity (SSIM) index [63], metric Q [64], a verage gradient (A G) [65], and blind image quality assessment (BIQA) [66]. The PSNR and SSIM assess the recovered image by com- paring it with the original image from the gray-level fidelity and the structure-le vel fidelity aspects, respecti v ely . The metric Q, A G, and BIQA assess the recov ered image without the reference image based on the human vision system. Gi ven a reference image ˜ X ∈ R m × n , the PSNR of a reconstructed image X ∈ R m × n is computed by the standard formula PSNR ( X , ˜ X ) = 10 log 10 N ˜ X 2 max k ˜ X − X k 2 F , (9) where N = mn denotes the number of the pixels in the image, and ˜ X max is the maximum pixel value of the original image. The SSIM of the estimated image X is defined by SSIM ( X , ˜ X ) = (2 µ X µ ˜ X + c 1 )(2 σ X ˜ X + c 2 ) ( µ 2 X + µ 2 ˜ X + c 1 )( σ 2 X + σ 2 ˜ X + c 2 ) , (10) where µ X and µ ˜ X represent the av erage gray values of the recov ered image X and the original clear image ˜ X , respectiv ely . σ X and σ ˜ X represent the standard de viation of X and ˜ X , respecti v ely . σ X ˜ X represents the covariance between X and ˜ X . The metric Q of an image is defined by Q ( X ) = s 1 s 1 − s 2 s 1 + s 2 , (11) where s 1 and s 2 are two singular values of the gradient matrix of the image X . The A G is computed by A G ( X ) = 1 ( m − 1)( n − 1) m − 1 X i =1 n − 1 X j =1 q (∆ 1 x 2 i,j + ∆ 2 x 2 i,j ) / 2 , (12) where ∆ 1 x i,j and ∆ 2 x i,j are the first dif ferences along both directions, respecti vely . The BIQA can be calculated IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 7 according to [66] and its code is available online 1 . In the real experiments, the performance is quantitati vely ev aluated by the metric Q, A G, and BIQA. In our e xperiments, the PSNR, SSIM, Q, A G, and BIQA v alues of a multispectral image are the average v alues of those for all bands. For all the fi ve indicators, the lar ger the value, the better the results. W ithout any special instructions, the parameters are set as following: the number of time series t = 4 , patch size w = 4 (patch size w should be multiples of t due to the rearrangement procedure), indicator thresholding value γ 2 is 0.91, radius of searching windo ws r = 100 ( r /t > 20 is recommended to make the search region large enough), searching step is 2 (half of patch size w ), penalty parameter β = 1 or 10 , and ε = 10 − 4 , 10 − 2 for the images whose value ranges are [0 , 255] and [0 , 1] respecti vely . Three of the most adv anced missing information reconstruction methods, HaLR TC [53], ALM-IPG [46], and PM-MTGSR [50], are compared with NL-LR TC. The parameters for these three compared methods are tuned to maximize the reconstruction PSNR v alue for each data set. All the experiments are performed under W indo ws 10 and Matlab V ersion 9.0.0.341360 (R2016a) running on a desktop with an Inter(R) Core(TM) i7-6700 CPU at 3.40 GHz and 16 GB of memory . C. Simulated Experiments In this section, the simulated experiments are presented to test NL-LR TC. The test data of Munich are “M012015” and “M022015”. T o assess the performance of NL-LR TC fully and efficiently , three subimages of “M012015” and one subimage of “M022015” shown in Fig. 7 are tested in the simulated experiments. These subimages are of size 512 × 512 × 3 . The structure and details of these subimages are dif ferent: the data set “Image 1” is cut from Fig. 4(b) (red square) and corre- sponding areas of Fig. 4(d)-(f) and mostly contains vegetation areas with relativ ely low contrast due to flat terrains; “Image 2” is cut from Fig. 4(b) (green square) and corresponding areas of Fig. 4(d)-(f) and mainly contains mountains, hills, and riv ers; and “Image 3” and “Image 4” are extracted from Fig. 4(b) (blue square) and (c) (white square), respecti vely , and the corresponding areas of Fig. 4(d)-(f). The data sets “Image 3” and “Image 4” contain both characteristics of “Image 1” and “Image 2”. The simulated clouds and stripes removal results are sho wn in Fig. 8, where Exps. 1–4 are for clouds remov al and Exps. 5 and 6 are for stripes remo v al. First, the cloud remov al e xperiments (Exps. 1–3) are un- dertaken using “Image 1”, “Image 2”, and “Image 3” shown in Fig. 7: the first column is the simulated corrupted data, and the other three columns are the supplementary data. The simulated cloud masks are generated by Matlab manually . The recovery results compared with HaLR TC, ALM-IPG, and PM-MTGSR are shown in Fig. 8. Exps. 1–3 indicate that PM-MTGSR and NL-LR TC succeed in estimating the missing entries. Howe ver , the results of HaLR TC and ALM- IPG hav e a noticeable dif ference between the known and the estimated areas. HaLR TC utilizes only the low-rankness of the 1 https://cn.mathworks.com/matlabcentral/fileexchange/30800-blind-image- quality-assessment-through-anisotropy Image 1 Image 2 Image 3 Image 4 Fig. 7. Simulated images cut from Fig. 4. From top to bottom, data set denoted as “Image 1”, “Image 2”, “Image 3”, and “Image 4”, respectiv ely . observed 4-order multitemporal data, in another word, it does not adopt patch strategy . ALM-IPG mainly studies the smooth relationships between the continuous temporal images, while the temporal series of the gi v en data is not continuous, i.e., the image acquired in the adjacent time is different largely . T o compare the results of PM-MTGSR and NL-LR TC in detail, a zoomed region is displayed in Fig. 9. From this figure, we can see that there are obvious edges between the estimated and known areas for the PM-MTGSR results, especially for Exps. 2 and 3. In contrast, the estimated areas of NL-LR TC hav e the similar tint with the kno wn areas, i.e., the estimated areas of NL-LR TC are in harmony with the kno wn areas. T o further contrast the original and reconstructed pixel values, scatter plots between the original and restoration pixels in the missing areas for Exps. 1–3 are sho wn in Fig. 10. The scatter plots results of Exp. 1 are visually better than those of Exp. 2 and Exp. 3 because the difference between the target image and reference images in Exp. 1 are slighter than those of Exp. 2 and Exp. 3. The points on the scatter plot of NL-LR TC show a more compact distribution than those of HaLR TC, ALM-IPG, and PM-MTGSR, b ut the difference between PM- MTGSR and NL-LR TC is not obvious for Exp. 1. For Exps. 2 and 3, the scatter plots of HaLR TC, ALM-IPG, and PM- MTGSR are obviously worse than those of ours, especially for the larger pixel values. This is consistent with the visual results of Exps. 2 and 3 shown in Figs. 8 and 9 where the reconstructed areas by HaLR TC, ALM-IPG, and PM-MTGSR are noticeably dif ferent from the known areas. Furthermore, one more simulated Beijing data shown in Fig. 5 is tested to demonstrate the ef fecti v eness of NL-LR TC (denoted as Exp. 4). In this experiment, we test ho w the four methods perform when the temporal dif ference between cloud- contained data and other temporal cloud-free data is not too large. “BJ092016” is adopted as the cloud-contained data. The IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 8 Exp. 1 Corrupted HaLRTC ALM-IPG PM-MTGSR NL-LRTC Exp. 2 Exp. 3 Exp. 4 Exp. 5 Exp. 6 Fig. 8. Simulated experiments, Exps. 1–6, for clouds and stripes remo val. Exps. 1–3 are cloud remov al for “Images 1–3”, respectively; Exp. 4 is cloud remov al for Beijing data (results for band 6 are sho wn in this figure); Exps. 5 and 6 are stripes remov al for “Image 3” and “Image 4”, respectively . results of band 6 reconstructed by the four methods for Exp. 4 are shown in Fig. 8. The Exp. 4 shows that the results by all the four methods are almost visually similar . One can get a similar conclusion from the scatter plots shown in Fig. 10. The points on the scatter plots of all the four methods are mostly distributed surrounding the blue diagonal, but there are also a few points deviating from the diagonal line. In this e xperiment, we can see that the results obtained by ALM-IPG are improv ed compared to Exp. 1–3 because ALM-IPG mainly studies the smoothness of adjacent temporal data. In conclusion, when the IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 9 Fig. 9. Zoom results of PM-MTGSR and NL-LR TC for Exps. 1–3. From left to right: original data, corrupted data, zoom part of results reconstructed by PM-MTGSR and NL-LR TC, respectiv ely . From top to bottom: zoom results for Exps. 1, 2, 3, respectively . For the corrupted data, the black areas are missing. temporal difference is not large, all the four methods perform in a similarly effecti ve manner . Besides the cloud removal experiments, the destriping ex- periments (Exps. 5 and 6) are also conducted using “Image 3” and “Image 4” shown in Fig. 7: the first column is the simulated corrupted data, and the other three columns are the supplementary data. Some regular diagonal and random vertical stripes are manually added into “Image 3” and “Image 4”, respecti vely , as shown in Fig. 8. It is obvious that the results by NL-LR TC are the best visually . The scatter plots comparing the original and reconstructed pixel v alues in the missing areas for Exps. 5 and 6 are shown in Fig. 10. The points on the scatter plot of HaLR TC and ALM-IPG results obviously de viate from the diagonal. The points for PM- MTGSR and NL-LR TC are better , but the points of PM- MTGSR distributed in the direction orthogonal to the diagonal line are wider than those of NL-LR TC. That means the scatter plots of our method are the best. Overall, the proposed method obtains the best results for the remov al of stripes. The quantitati ve comparison is shown in T ab. II. The table shows that all the four methods can recover better results compared to the corrupted image itself. PSNR and SSIM ev aluate the recovered image by comparing with the ground truth image. By analyzing the PSNR and SSIM results, HaL- R TC obtains the worst results for all the experiments. This is because HaLR TC only takes advantage of the low-rankness of the observed data. ALM-IPG obtains better results than HaLR TC, because it considers the lo w-rankness and temporal continuous property simultaneously . Howe ver , it assumes the smoothness of adjacent temporal images, the results depend on the similarity of the adjacent temporal images. This algorithm is suited to process high-temporal-resolution images, such as videos. PM-MTGSR obtains better results than HaLR TC and ALM-IPG because it makes use of the patch similarity . NL- LR TC also takes the patch similarity into consideration and makes the best use of lo w-rankness of the three dif ferent dimensions. Thus NL-LR TC obtains the best results. For Exp. 4, since the dif ference between the cloud-contained and reference data is not great, the result of ALM-IPG is better than those of Exps. 1–3. Although the PSNR value of PM- MTGSR is higher than that of NL-LR TC, the difference is slight. Moreover , the SSIM v alue of NL-LR TC is higher than that of PM-MTGSR. The Q, A G, and BIQA results also show NL-LR TC obtains the best results for Exps. 1–4. The elapsed times of PM-MTGSR and NL-LR TC are longer than those of HaLR TC and ALM-IPG because they are patch-based methods in which searching similar patches costs much more time. NL- LR TC is much faster than PM-MTGSR because PM-MTGSR processes multispectral images band by band, while NL-LR TC reconstructs the missing areas of all band at the same time. In conclusion, the proposed method performs the best for cloud and stripe remov al when the temporal dif ference is large and can also obtain promising results when the temporal difference is slight. Next, two simulated data containing more than one piece of cloud are tested. The recovered results are shown in Fig. 11, which are for the “Image 3” taken by Landsat-8 (Exp. 7) and Beijing data taken by Sentinel-2 (Exp. 8). Exp. 7 and Exp. 8 perform the similar results with Exps. 1–3. and Exp. 4, respectiv ely . The results recovered by PM-MTGSR and NL- LR TC in Exp. 7 are visually similar , b ut are visually better than those reconstructed by HaLR TC and ALM-IPG. Exp. 8 shows all the four methods obtain the visually similar results. T ab . III shows that, for both Exps. 7 and 8, NL-LR TC obtains the best quantitati v e results. At last, we analyze the impact of the number of time series on the reconstruction performance by changing the number ( t = 2 , 4 , · · · , 16 ). The test data were taken over Munich by Landsat-8 on between December , 2014 and April, 2017 2 . The SSIM and PSNR v alues with respect to the number of time series are displayed in Fig. 12. This figure shows that reconstruction results are becoming better with increasing of the number of time series. When the number of time series reach to a lar ge amount, the PSNR and SSIM values reach to the highest with a little fluctuation. This is because more temporal data not only provide more correlative information but also contain more interference information especially when the acquisition times of the cloud-contained and reference data are far form each other . D. Real Experiments In this section, real-data experiments are undertak en. The experimental data are “M102014”, “BJ072016”, and “EU082017”. The cloud detection is not our focus in this work and is complex for dif ferent kinds of atmospheric conditions. For the Landsat data “M102014”, the cloud is detected via a modified version of the thresholding-based cloud detection method proposed in [46]; see Appendix A for more details. Beijing data “BJ072016” contains shadows that cannot be de- tected by a simple thresholding method. Thus, for “BJ072016”, 2 The sixteen test data were taken during four years: 2014 (December), 2015 (January-April, June, August, October), 2016 (April-September), and 2017 (January , February). IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 10 Exp. 1 Exp. 2 Exp. 3 Exp. 4 Exp. 5 Exp. 6 Fig. 10. Scatter plots between the original and restoration pixels in the missing areas for Exps. 1–6. the mask for clouds and their shado ws is manually drawn. For the Eure data “EU082017”, the cloud detection is processed by MA Y A [62]. The corresponding recovery results are shown in Figs. 13, 14, and 15. For “M102014” (see Fig. 13), the color composite images of reconstruction areas obtained by HaLR TC, ALM-IPG, and PM-MTGSR are obviously different from the known areas. The reconstruction area of NL-LR TC shows a more natural visual effect. For “BJ072016” (see Fig. 14), the recovery results by ALM-IPG hav e obvious stairs in the edge of missing and known areas. HaLR TC and PM- MTGSR fail in reconstructing clear details. NL-LR TC sho ws better results containing more clear details and being more natural compared to the kno wn area. F or “EU082017” (see Fig. 15), the reco very results by HaLR TC, ALM-IPG, and PM- MTGSR are visually worse than that by NL-LR TC, that means the missing areas recovered by NL-LR TC are in harmony with the know areas. Moreover , quantitativ e results for Figs. 13, 14, and 15 are shown in T ab . IV. For the Beijing and Eure real IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 11 T ABLE II Q UA N TI TA TI V E R E S U L T S C O MPA R IS O N F O R E X P S . 1 – 6 . Methods Exp. 1 Exp. 2 SSIM PSNR Q A G BIQA Time(s) SSIM PSNR Q A G BIQA Time(s) Corrupted 0.8107 9.35 - - - - 0.6237 6.54 - - - - HaLR TC 0.9660 36.78 0.0230 0.0250 0.0036 39.53 0.8357 24.36 0.0342 0.0303 0.0075 44.45 ALM-IPG 0.9708 38.05 0.0232 0.0252 0.0036 49.94 0.8575 25.81 0.0345 0.0318 0.0074 41.34 PM-MTGSR 0.9722 38.40 0.0230 0.0260 0.0036 4028.05 0.8594 25.85 0.0373 0.0330 0.0080 4032.23 NL-LR TC 0.9769 38.95 0.0245 0.0266 0.0038 478.73 0.8975 27.32 0.0443 0.0375 0.0090 1216.38 Methods Exp. 3 Exp. 4 SSIM PSNR Q A G BIQA Time(s) SSIM PSNR Q A G BIQA Time(s) Corrupted 0.5002 4.97 - - - - 0.4851 22.02 - - - - HaLR TC 0.8117 25.39 0.0235 0.0271 0.0052 41.16 0.9664 39.09 0.0092 0.0080 0.0014 20.05 ALM-IPG 0.8341 26.26 0.0260 0.0291 0.0051 41.00 0.9677 40.41 0.0095 0.0080 0.0013 18.81 PM-MTGSR 0.8419 26.67 0.0286 0.0305 0.0058 4056.01 0.9704 41.17 0.0086 0.0082 0.0013 2391.16 NL-LR TC 0.8629 27.18 0.0403 0.0351 0.0068 1274.66 0.9707 41.13 0.0085 0.0085 0.0017 345.75 Methods Exp. 5 Exp. 6 SSIM PSNR Q A G BIQA Time(s) SSIM PSNR Q A G BIQA Time(s) Corrupted 0.6848 21.70 - - - - 0.4920 15.67 - HaLR TC 0.9308 29.90 0.0399 0.0363 0.0110 37.83 0.8749 24.23 0.2298 0.0586 0.0151 43.35 ALM-IPG 0.9341 30.75 0.0407 0.0363 0.0106 41.97 0.9285 27.78 0.0873 0.0505 0.0081 41.66 PM-MTGSR 0.9407 30.69 0.0356 0.0367 0.0106 4148.10 0.9381 28.73 0.0694 0.0466 0.0068 4169.75 NL-LR TC 0.9846 36.43 0.0299 0.0355 0.0073 1428.18 0.9679 31.71 0.0523 0.0483 0.0064 1338.02 Exp. 7 Corrupted HaLRTC ALM-IPG PM-MTGSR NL-LRTC Exp. 8 Fig. 11. Cloud remov al results for “Image 3” and Beijing data containing more than one piece of cloud. Exp. 7 is for “Image 3” and Exp. 8 is for Beijing data. T ABLE III Q UA N TI TA TI V E R E S U L T S C O MPA R IS O N F O R E X P S . 7 A N D 8 . HaLR TC ALM-IPG PM-MTGSR NL-LR TC Exp. 7 SSIM 0.8323 0.8439 0.8538 0.8862 PSNR 25.32 26.26 26.29 27.59 Exp. 8 SSIM 0.9618 0.9645 0.9692 0.9711 PSNR 38.25 39.81 40.96 41.61 data, all the three index values (Q, A G, and BIQA) for NL- LR TC are better than those for HaLR TC, ALM-IPG, and PM- MTGSR. For the Landsat-8 real data, the Q and BIQA values for NL-LR TC are worse than those for HaLR TC, ALM-IPG, and PM-MTGSR. Howe v er , the difference is slight. The AG value for NL-LR TC is better than the other three compared methods. The real-data experiments also demonstrate that the 2 4 6 8 10 12 14 16 # of time series 0.85 0.9 0.95 1 SSIM 2 4 6 8 10 12 14 16 # of time series 20 25 30 35 40 45 PSNR Fig. 12. SSIM and PSNR values with respect to the numbers of time series. proposed method is effecti ve. V . C O N C L U S I O N In this paper, a non-local lo w-rank tensor completion (NL- LR TC) method has been proposed to reconstruct the missing IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 12 Original HaLRTC ALM-IPG PM-MTGSR NL-LRTC Fig. 13. Results for real Landsat experiment. Original Cloud mask HaLRTC ALM-IPG PM-MTGSR NL-LRTC Fig. 14. Results for real Sentinel-2 data taken over Beijing. T ABLE IV Q UA N TI TA TI V E R E S U L T S C O MPA R IS O N F O R F I G S . 1 3 , 1 4 , A N D 1 5 . HaLR TC ALM-IPG PM-MTGSR NL-LR TC Fig. 13 Q 0.0384 0.0384 0.0366 0.0371 A G 0.0331 0.0335 0.0337 0.0339 BIQA 0.0058 0.0056 0.0058 0.0056 Fig. 14 Q 0.0132 0.0100 0.0073 0.0144 A G 0.0060 0.0069 0.0050 0.0088 BIQA 0.0012 0.0009 0.0005 0.0012 Fig. 15 Q 0.0264 0.0354 0.0258 0.0404 A G 0.0185 0.0196 0.0158 0.0222 BIQA 0.0053 0.0068 0.0057 0.0070 information in the multitemporal remotely sensed images. By proposing a non-con ve x approximation for tensors rank, all the three domains (spatial, spectral, and temporal) relationships were considered in NL-LR TC. T o take advantage of the spatial correlations, we grouped the 3-order similar patches into a 4- Original Cloud mask HaLRTC ALM-IPG PM-MTGSR NL-LRTC Fig. 15. Results for real Sentinel-2 data taken over Eure, France. In this figure, the images are sho wn in color format using bands 4, 3, and 2. order tensor and considered the tensor lo w-rankness. Because NL-LR TC made use of the global correlations of all the three domains, it is good at processing not only the temporally contiguous data but also the data that hav e large differences between the adjacent temporal images regarding the character- istics and conditions of the Earth’ s surface. In the simulations with various image data sets, NL-LR TC showed comparable or better results than HaLR TC, ALM-IPG, and PM-MTGSR, which are three of the state-of-the-art algorithms. For the real- data experiments, our method obtained visually more natural and quantitati vely better reconstruction results. A P P E N D I X A C L O U D D E T E C T I O N In this section, we present an automatic thresholding method for cloud detection motiv ated by the algorithm of [46]. This method assumes that most cloud values in the remotely sensed images are larger than other cloud free v alues, i.e., clouds are predominantly white. Giv en a 4-order observ ation remotely sensed image Y ∈ R m × n × b × t , where m × n denotes the number of pix els of remotely sensed images, b denotes the number of spectral channels of remote sensors, and t is the number of time series. In this research, we talk about how to restore the image at one time according to the other times. Suppose Y t 1 taken at time t 1 is the cloud contained image, then the other images (denoted as Y ˆ t 1 ) are the references. The cloud detector is to produce a set of indices Ω ∈ R m × n × b × t , where the position ( i, j, k, l ) ∈ Z m × Z n × Z b × Z t is co vered by cloud if Ω i,j,k,l = 0 and is cloud free if Ω i,j,k,l = 1 . Note that all the spectral bands of the practical remotely sensed images taken at the same position and same period will be co vered by the same clouds, i.e, Ω ( i, j, k , t 1 ) = 0 , ∀ k ∈ { 1 , 2 , · · · , b } if exist k 1 ∈ { 1 , · · · , b } subject to Ω ( i, j, k 1 , t 1 ) = 0 . In this research, there are some other cloud free references. The key point of the detection method is to maximize the similarity of cloud contained and free images in the existing IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 13 region Ω . Gi ven a similar function f ( x , y ) , find the indices set Ω by optimizing the follo wing problem: ˜ Ω = max Ω f ( Y t 1 Ω , Y ˆ t 1 Ω ) . (13) The function can be an y similarity functions, such as correla- tion coefficients, cosine coefficients, generalized Dice coeffi- cients, and generalized Jaccard coefficients. Here, we use the correlation coefficients. Besides the mentioned method (13), one can also minimize the distance function such as Euclidean distance, the mean absolute error (MAE), and the mean relati v e error (MRE) to get the indices set Ω . The thresholding produce is detailedly summarized in Algorithm 3. Algorithm 3 Thresholding method. Input: temporal sequence of cloudy images Y , parameter step for increase the thresholding value s . 1: Obtain the initial indices set Ω 0 = Y t 1 > 0 ; 2: Set ˜ Ω = Ω 0 and γ 1 = 0 ; 3: while f ( Y t 1 Ω , Y ˆ t 1 Ω ) increase do 4: Update the thresholding value γ 1 = γ 1 + s ; 5: Calculate the correlation f ( Y t 1 Ω , Y ˆ t 1 Ω ) . 6: end while Output: ˜ Ω : initial guess for index set. The above described thresholding method regards the sta- tionary white background as the cloud, as seen in Fig. 17. It is better that these white objects remain in the reconstructed images. Fortunately , the clouds usually cover a big continuous area, this fact motiv ates us to delete the discrete points. T o this end, we propose a KNN-lik e method. In detail, the pixel in ˜ Ω c is not regarded as cloud if most of its surrounding pix els are not cloud, i.e., most of its neighbor pixels are in ˜ Ω . The procedure for finding the white background is shown in Fig. 16. The two stages of cloud detection is summarized in Algo- rithm 4 belo w . Algorithm 4 Cloud detection. Input: temporal sequence of cloudy images Y , parameter thresholding value γ 1 , and parameter r in the similar k- nearest-neighbors search. 1: Obtain the initial guess Ω via Algorithm 3: Ω = ˜ Ω ; 2: for ( m, n, : , l ) / ∈ ˜ Ω do 3: Extract the patch P with the center ( m, n, : , l ) and a radius of r 1 : P ( i, j, : , l ) = Ω ( i, j, : , l ) , where k ( i, j ) − ( m, n ) k ∞ < r 1 ; 4: Calculate the percent ( p ) of the cloud pix els number; 5: if p < 0 . 5 then 6: Ω = Ω ∪ { ( m, n, : , l ) } . 7: end if 8: end f or Output: Ω : index set of non-cloudy pixels. Fig. 16. Modified cloud detection procedure. The black pixels denote 0 (cloud), others denote 1. In this figure, r =3. For the red target pixel, it should be the white background rather than cloud, because in the search window , most pixels are 1. While for the blue target pixel, it should be cloud. (a) Cloud contained (b) Cloud remov al (c) Alg. 3 (d) Alg. 4 Fig. 17. Illustration of the proposed cloud detection procedure. (a) cloud contained image, (b) removing the detected cloud, (c) cloud detected via Alg. 3, (c) cloud detected via Alg. 4. A C K N O W L E D G M E N T The authors would like to thank Prof. H. Shen and X. Li from W uhan Univ ersity for sharing the codes of the PM- MTGSR method. R E F E R E N C E S [1] M.-D. Iordache, J. M. Bioucas-Dias, and A. Plaza, “Sparse unmixing of hyperspectral data, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 49, no. 6, pp. 2014–2039, 2011. [2] M.-D. Iordache, J. M. Bioucas-Dias and A. Plaza, “T otal v ariation spatial regularization for sparse hyperspectral unmixing, ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 50, no. 11, pp. 4484–4502, 2012. [3] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P . Gader , and J. Chanussot, “Hyperspectral unmixing overvie w: Geometrical, statistical, and sparse regression-based approaches, ” IEEE Journal of Selected T opics in Applied Earth Observations and Remote Sensing , vol. 5, no. 2, pp. 354–379, 2012. [4] X.-L. Zhao, F . W ang, T .-Z. Huang, M. K. Ng, and R. J. Plemmons, “Deblurring and sparse unmixing for hyperspectral images, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 51, no. 7-1, pp. 4045–4058, 2013. [5] N. Y okoya, T . Y airi, and A. Iwasaki, “Coupled nonne gativ e matrix factorization unmixing for hyperspectral and multispectral data fusion, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 50, no. 2, pp. 528–537, Feb 2012. [6] N. Y okoya, J. Chanussot, and A. Iw asaki, “Nonlinear unmixing of hyperspectral data using semi-nonnegativ e matrix factorization, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 52, no. 2, pp. 1430–1437, Feb 2014. IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 14 [7] J. Bieniarz, E. Aguilera, X. X. Zhu, R. Muller, and P . Reinartz, “Joint Sparsity Model for Multilook Hyperspectral Image Unmixing, ” IEEE Geoscience and Remote Sensing Letters , vol. 12, no. 4, pp.696–700, 2014. [8] D. Lunga, S. Prasad, M. M. Crawford, and O. Ersoy , “Manifold-learning- based feature extraction for classification of h yperspectral data: A revie w of advances in manifold learning, ” IEEE Signal Pr ocessing Magazine , vol. 31, no. 1, pp. 55–66, Jan 2014. [9] Y . T arabalka, M. Fauvel, J. Chanussot, and J. A. Benediktsson, “SVM- and MRF-based method for accurate classification of hyperspectral images, ” IEEE Geoscience and Remote Sensing Letter s , vol. 7, no. 4, pp. 736–740, Oct 2010. [10] J. C. Harsanyi and C.-I. Chang, “Hyperspectral image classification and dimensionality reduction: an orthogonal subspace projection approach, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 32, no. 4, pp. 779–785, Jul 1994. [11] T . Matsuki, N. Y okoya, and A. Iwasaki, “Hyperspectral tree species classification of japanese complex mixed forest with the aid of lidar data, ” IEEE Journal of Selected T opics in Applied Earth Observations and Remote Sensing , vol. 8, no. 5, pp. 2177–2187, May 2015. [12] P . Ghamisi, R. Souza, J. A. Benediktsson, X. X. Zhu, L. Rittner , and R. A. Lotufo, “Extinction profiles for the classification of remote sensing data, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 54, no. 10, pp. 5631–5645, Oct 2016. [13] L. Mou, P . Ghamisi, X. Zhu, “Deep Recurrent Neural Networks for Hyperspectral Image Classification, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 55, no. 7, pp. 3639–3655,2017. [14] L. Mou, P . Ghamisi, X. Zhu, “Unsupervised Spectral-Spatial Feature Learning via Deep Residual Conv-Decon v Network for Hyperspectral Image Classification, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 56, no. 1, pp. 391 – 406,2018. [15] D. Manolakis, D. Marden, and G. A. Shaw , “Hyperspectral image pro- cessing for automatic target detection applications, ” Lincoln Laboratory Journal , vol. 14, no. 1, pp. 79–116, 2003. [16] R. M. Willett, M. F . Duarte, M. A. Davenport, and R. G. Baraniuk, “Sparsity and structure in hyperspectral imaging : Sensing, reconstruc- tion, and target detection, ” IEEE Signal Pr ocessing Magazine , vol. 31, no. 1, pp. 116–126, Jan 2014. [17] N. M. Nasrabadi, “Hyperspectral target detection : An overvie w of current and future challenges, ” IEEE Signal Pr ocessing Magazine , vol. 31, no. 1, pp. 34–44, Jan 2014. [18] N. Y okoya, N. Miyamura, and A. Iwasaki, “Detection and correction of spectral and spatial misregistrations for hyperspectral data using phase correlation method, ” Applied Optics , vol. 49, no. 24, pp. 4568–4575, Aug 2010. [19] N. Y okoya and A. Iw asaki, “Object detection based on sparse repre- sentation and hough v oting for optical remote sensing imagery , ” IEEE Journal of Selected T opics in Applied Earth Observations and Remote Sensing , vol. 8, no. 5, pp. 2053–2062, May 2015. [20] M. Shahzad and X. X. Zhu, “ Automatic detection and reconstruction of 2-D/3-D building shapes from spaceborne TomoSAR point clouds, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 54, no. 3, pp. 1292–1310, March 2016. [21] H. Shen, X. Li, L. Zhang, D. T ao, and C. Zeng, “Compressed sensing- based inpainting of Aqua moderate resolution imaging spectroradiometer band 6 using adaptiv e spectrum-weighted sparse bayesian dictionary learning, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 52, no. 2, pp. 894–906, 2014. [22] L. W ang, J. J. Qu, X. Xiong, X. Hao, Y . Xie, and N. Che, “ A new method for retrieving band 6 of Aqua MODIS, ” IEEE Geoscience and Remote Sensing Letters , vol. 3, no. 2, pp. 267–270, April 2006. [23] C. Zeng, H. Shen, and L. Zhang, “Recovering missing pixels for Landsat ETM+ SLC-off imagery using multi-temporal regression analysis and a regularization method, ” Remote Sensing of Envir onment , vol. 131, pp. 182 – 194, 2013. [24] H. Shen, X. Li, Q. Cheng, C. Zeng, G. Y ang, H. Li, and L. Zhang, “Missing information reconstruction of remote sensing data: A technical revie w , ” IEEE Geoscience and Remote Sensing Magazine , vol. 3, no. 3, pp. 61–85, 2015. [25] C.-H. Lin, K.-H. Lai, Z.-B. Chen, and J.-Y . Chen, “Patch-based informa- tion reconstruction of cloud-contaminated multitemporal images, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 52, no. 1, pp. 163–174, Jan 2014. [26] C. Y u, L. Chen, L. Su, M. Fan, and S. Li, “Kriging interpolation method and its application in retrieval of MODIS aerosol optical depth, ” in 19th International Conference on Geoinformatics , June 2011, pp. 1–6. [27] A. Maalouf, P . Carre, B. Augereau, and C. Fernandez-Maloigne, “ A bandelet-based inpainting technique for clouds removal from remotely sensed images, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 47, no. 7, pp. 2363–2371, July 2009. [28] C. Ballester, M. Bertalmio, V . Caselles, G. Sapiro, and J. V erdera, “Filling-in by joint interpolation of vector fields and gray levels, ” IEEE T r ansactions on Image Processing , vol. 10, no. 8, pp. 1200–1211, Aug 2001. [29] A. Bugeau, M. Bertalmio, V . Caselles, and G. Sapiro, “ A comprehensiv e framew ork for image inpainting, ” IEEE T ransactions on Image Pr ocess- ing , vol. 19, no. 10, pp. 2634–2645, Oct 2010. [30] W . He, H. Zhang, L. Zhang, and H. Shen, “T otal-variation-re gularized low-rank matrix factorization for hyperspectral image restoration, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 54, no. 1, pp. 178–188, 2016. [31] Q. Cheng, H. Shen, L. Zhang, and P . Li, “Inpainting for remotely sensed images with a multichannel nonlocal total v ariation model, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 52, no. 1, pp. 175–187, 2014. [32] H. Shen and L. Zhang, “ A MAP-based algorithm for destriping and in- painting of remotely sensed images, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 47, no. 5, pp. 1492–1502, May 2009. [33] H. Shen, Y . Liu, T . Ai, Y . W ang, and B. W u, “Universal reconstruction method for radiometric quality improvement of remote sensing images, ” International J ournal of Applied Earth Observation and Geoinforma- tion , vol. 12, no. 4, pp. 278 – 286, 2010. [34] A. Criminisi, P . P ´ erez, and K. T oyama, “Region filling and object remov al by exemplar-based image inpainting, ” IEEE T ransactions on Image Processing , vol. 13, no. 9, pp. 1200–1212, Sept 2004. [35] A. A. Efros and T . K. Leung, “T exture synthesis by non-parametric sampling, ” in Pr oceedings of the IEEE International Conference on Computer V ision (ICCV) , 1999. [36] P . Rakwatin, W . T akeuchi, and Y . Y asuoka, “Restoration of Aqua MODIS band 6 using histogram matching and local least squares fitting, ” IEEE Tr ansactions on Geoscience and Remote Sensing , vol. 47, no. 2, pp. 613–627, Feb 2009. [37] H. Shen, C. Zeng, and L. Zhang, “Recovering reflectance of A QU A MODIS band 6 based on within-class local fitting, ” IEEE Journal of Selected T opics in Applied Earth Observations and Remote Sensing , vol. 4, no. 1, pp. 185–192, March 2011. [38] I. Gladkov a, M. D. Grossberg, F . Shahriar , G. Bonev , and P . Romanov , “Quantitativ e restoration for MODIS band 6 on Aqua, ” IEEE Tr ansac- tions on Geoscience and Remote Sensing , vol. 50, no. 6, pp. 2409–2416, June 2012. [39] X. Li, H. Shen, L. Zhang, H. Zhang, and Q. Y uan, “Dead pixel comple- tion of aqua MODIS band 6 using a robust M-estimator multire gression, ” IEEE Geoscience and Remote Sensing Letters , vol. 11, no. 4, pp. 768– 772, April 2014. [40] C. Zeng, H. Shen, M. Zhong, L. Zhang, and P . W u, “Reconstructing MODIS LST based on multitemporal classification and rob ust regres- sion, ” IEEE Geoscience and Remote Sensing Letters , vol. 12, no. 3, pp. 512–516, March 2015. [41] Y . Julien and J. A. Sobrino, “Comparison of cloud-reconstruction methods for time series of composite NDVI data, ” Remote Sensing of En vir onment , vol. 114, no. 3, pp. 618 – 625, 2010. [42] W . Zhu, Y . Pan, H. He, L. W ang, M. Mou, and J. Liu, “ A changing- weight filter method for reconstructing a high-quality NDVI time series to preserve the integrity of vegetation phenology , ” IEEE T ransactions on Geoscience and Remote Sensing , vol. 50, no. 4, pp. 1085–1094, April 2012. [43] A. Savitzky and M. J. E. Golay , “Smoothing and differentiation of data by simplified least squares procedures. ” Analytical Chemistry , vol. 36, no. 8, pp. 1627–1639, 1964. [44] L. Lorenzi, F . Melgani, and G. Mercier, “Missing-area reconstruction in multispectral images under a compressiv e sensing perspective, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 51, no. 7, pp. 3998–4008, July 2013. [45] X. Li, H. Shen, L. Zhang, H. Zhang, Q. Y uan, and G. Y ang, “Recovering quantitativ e remote sensing products contaminated by thick clouds and shadows using multitemporal dictionary learning, ” IEEE T r ansactions on Geoscience and Remote Sensing , vol. 52, no. 11, pp. 7086–7098, Nov 2014. [46] J. W ang, P . A. Olsen, A. R. Conn, and A. C. Lozano, “Removing clouds and recovering ground observations in satellite image sequences via temporally contiguous robust matrix completion, ” in Pr oceedings of the IEEE Conference on Computer V ision and P attern Recognition (CVPR) , 2016. IEEE TRANSA CTIONS ON GEOSCIENCE AND REMO TE SENSING, IN PRESS 15 [47] Q. Cheng, H. Shen, L. Zhang, Q. Y uan, and C. Zeng, “Cloud remov al for remotely sensed images by similar pixel replacement guided with a spatio-temporal MRF model, ” ISPRS Journal of Photogrammetry and Remote Sensing , vol. 92, pp. 54 – 68, 2014. [48] S. Benabdelkader and F . Melgani, “Contextual spatiospectral postrecon- struction of cloud-contaminated images, ” IEEE Geoscience and Remote Sensing Letters , vol. 5, no. 2, pp. 204–208, April 2008. [49] X. Li, H. Shen, L. Zhang, and H. Li, “Sparse-based reconstruction of missing information in remote sensing images from spectral/temporal complementary information, ” ISPRS Journal of Photogrammetry and Remote Sensing , vol. 106, pp. 1–15, 2015. [50] X. Li, H. Shen, H. Li, and L. Zhang, “Patch matching-based mul- titemporal group sparse representation for the missing information reconstruction of remote-sensing images, ” IEEE Journal of Selected T opics in Applied Earth Observations and Remote Sensing , vol. 9, no. 8, pp. 3629–3641, 2016. [51] J. Kang, Y . W ang, M. Schmitt, X. Zhu, “Object-based Multipass InSAR via Robust Low Rank T ensor Decomposition, ” IEEE T ransactions on Geoscience and Remote Sensing , in press. [52] T . G. Kolda and B. W . Bader, “T ensor decompositions and applications, ” SIAM Review , vol. 51, no. 3, pp. 455–500, 2009. [53] J. Liu, P . Musialski, P . W onka, and J. Y e, “T ensor completion for estimating missing values in visual data, ” IEEE T ransactions on P attern Analysis and Machine Intelligence , vol. 35, no. 1, pp. 208–220, 2013. [54] T .-Y . Ji, T .-Z. Huang, X.-L. Zhao, T .-H. Ma, and G. Liu, “T ensor completion using total variation and low-rank matrix factorization, ” Information Sciences , vol. 326, pp. 243–257, 2016. [55] M. M. Deza and E. Deza, in Encyclopedia of Distances . Springer, 2009, pp. 1–583. [56] T .-Y . Ji, T .-Z. Huang, X.-L. Zhao, T .-H. Ma, and L.-J. Deng, “ A non-con ve x tensor rank approximation for tensor completion, ” Applied Mathematical Modelling , vol. 48, pp. 410–422, 2017. [57] S. Gu, L. Zhang, W . Zuo, and X. Feng, “W eighted nuclear norm minimization with application to image denoising, ” in Pr oceedings of the IEEE Conference on Computer V ision and P attern Recognition (CVPR) , 2014. [58] Z. Lin, M. Chen, and Y . Ma, “The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices, ” arXiv preprint arXiv:1009.5055 , 2010. [59] B. He, M. T ao, and X. Y uan, “ Alternating direction method with gaussian back substitution for separable con vex programming, ” SIAM Journal on Optimization , vol. 22, no. 2, pp. 313–340, 2012. [60] X.-L. Zhao, F . W ang, and M. K. Ng, “ A new conv ex optimization model for multiplicati ve noise and blur remov al, ” SIAM Journal on Imaging Sciences , vol. 7, no. 1, pp. 456–475, 2014. [61] W . Dong, G. Shi, X. Li, Y . Ma, and F . Huang, “Compressive sensing via nonlocal lo w-rank regularization, ” IEEE T r ansactions on Imag e Pr ocessing , vol. 23, no. 8, pp. 3618–3632, 2014. [62] V . Lonjou, C. Desjardins, O. Hagolle, B. Petrucci, T . Tremas, M. Dejus, A. Makarau, and S. Auer, “Maccs-atcor joint algorithm (MAJ A), ” in Pr oceedings of SPIE Remote Sensing , vol. 10001, 2016, pp. 1–13. [63] Z. W ang, A. C. Bovik, H. R. Sheikh, and E. P . Simoncelli, “Image quality assessment: from error visibility to structural similarity , ” IEEE T r ansactions on Image Pr ocessing , vol. 13, no. 4, pp. 600–612, 2004. [64] X. Zhu and P . Milanfar, “ Automatic parameter selection for denoising algorithms using a no-reference measure of image content, ” IEEE T r ansactions on Image Processing , vol. 19, no. 12, pp. 3116–3132, Dec 2010. [65] Z. Li, Z. Jing, X. Y ang, and S. Sun, “Color transfer based remote sensing image fusion using non-separable wavelet frame transform, ” P attern Recognition Letters , vol. 26, no. 13, pp. 2006–2014, 2005. [66] S. Gabarda and G. Crist ´ obal, “Blind image quality assessment through anisotropy , ” Journal of the Optical Society of America A , vol. 24, no. 12, pp. B42–B51, Dec 2007. T eng-Y u Ji receiv ed the B.S. de gree from the School of Mathematical Sciences, Uni v ersity of Electronic Science and T echnology of China, Chengdu, China, in 2012, where he is currently pursuing the Ph.D. degree with the School of Mathematical Sciences. His current research interests include tensor de- composition and applications, including tensor com- pletion and remotely sensed image reconstruction. Naoto Y okoya (S’10–M’13) receiv ed the M.Sc. and Ph.D. degrees in aerospace engineering from the Univ ersity of T ok yo, T okyo, Japan, in 2010 and 2013, respectively . From 2012 to 2013, he was a Research Fellow with Japan Society for the Promotion of Science, T okyo, Japan. Since 2013, he is an Assistant Pro- fessor with the Uni versity of T okyo. From 2015 to 2017, he was also an Alexander v on Humboldt Research Fellow with the German Aerospace Center (DLR), Oberpfaffenhofen, and T echnical University of Munich (TUM), Munich, Germany . His research interests include image analysis and data fusion in remote sensing. Since 2017, he is a Co-chair of IEEE Geoscience and Remote Sensing Image Analysis and Data Fusion T echnical Committee. Xiao Xiang Zhu (S’10-M’12-SM’14) received the Master (M.Sc.) degree, her doctor of engineering (Dr .-Ing.) degree and her “Habilitation” in the field of signal processing from T echnical University of Munich (TUM), Munich, Germany , in 2008, 2011 and 2013, respectively . She is currently the Professor for Signal Process- ing in Earth Observation (www .sipeo.bgu.tum.de) at T echnical University of Munich (TUM) and German Aerospace Center (DLR); the head of the T eam Sig- nal Analysis at DLR; and the head of the Helmholtz Y oung Inv estigator Group ”SiPEO” at DLR and TUM. Prof. Zhu was a guest scientist or visiting professor at the Italian National Research Council (CNR- IREA), Naples, Italy , Fudan Univ ersity , Shanghai, China, the Uni versity of T okyo, T okyo, Japan and Univ ersity of California, Los Angeles, United States in 2009, 2014, 2015 and 2016, respectiv ely . Her main research interests are remote sensing and Earth observation, signal processing, machine learning and data science, with a special application focus on global urban mapping. Dr . Zhu is a member of young academy (Junge Akademie/Junges Kolle g) at the Berlin-Brandenbur g Academy of Sciences and Humanities and the German National Academy of Sciences Leopoldina and the Bavarian Academy of Sciences and Humanities. She is an associate Editor of IEEE Transactions on Geoscience and Remote Sensing. Ting-Zhu Huang recei v ed the B. S., M. S., and Ph. D. degrees in Computational Mathematics from the Department of Mathematics, Xian Jiaotong Uni- versity , Xian, China. He is currently a professor in the School of Mathematical Sciences, UESTC. He is currently an editor of The Scientific W orld Journal, Advances in Numerical Analysis, J. Appl. Math., J. Pure and Appl. Math.: Adv . and Appl., J. Electronic Sci. and T ech. of China, etc. His current research interests include scientific computation and appli- cations, numerical algorithms for image processing, numerical linear algebra, preconditioning technologies, and matrix analysis with applications, etc.