Despeckling Multitemporal Polarimetric SAR Data Based on Tensor Decomposition

Despeckling is an essential task in polarimetric synthetic aperture radar (PolSAR) image processing. Most of the existing filters developed for multitemporal synthetic aperture radar (SAR) images make use of either real or complex information. Real information refers to amplitude or intensity values, whereas complex information refers to complex covariance matrix (CCM) derived from either interferometric SAR (InSAR) or PolSAR data. The InSAR CCM is formed using images of the same polarimetric channel but acquired at different dates, and the PolSAR CCM contains information acquired simultaneously in different polarimetric channels. Therefore, these despeckling methods may present good performance in some applications and scenes but fail in other cases, due to differences in input sources. In order to achieve a more robust result in all cases, we develop a method for multitemporal polarimetric SAR data filtering based on tensor decomposition, which aims at improving the identification of homogeneous pixels for spatially adaptive filtering. The key element of this approach consists of exploiting tensor theory to construct a new CCM that contains both polarimetric and interferometric information, as well as multitemporal information for each pixel. The effectiveness of the proposed method and its performance are evaluated with simulated and real SAR data in comparison with several established methods.


I. INTRODUCTION
M ODERN synthetic aperture radar (SAR) systems provide an unprecedented ability to acquire polarimetric images over the same area with high revisit frequency.Polarimetric SAR (PolSAR) is an advancement over the conventional singlepolarization SAR because of the additional information content of the various polarimetric channels measured.Polarization diversity enhances discrimination capabilities, which can benefit in, for example, forest, agriculture, cryosphere, urban, and ocean applications [1].Unfortunately, SAR images are intrinsically contaminated by speckle, which complicates the aforementioned applications.
In the past decades, single-polarization SAR images have been vastly used for applications such as ground surface deformation monitoring [2], [3], [4] and land cover mapping [5].As a result, speckle filters are mainly focused on intensity or interferometric complex information.The simplest despeckling method is the spatial multilooking process, which consists in averaging neighbor pixels within a spatial window.If the window is predefined, e.g., rectangular, it achieves good smoothing at the cost of blurring spatial details and fine structures.To avoid selecting pixels that belong to different targets, adaptive techniques were proposed.Since 1980s, filters based on the minimum mean square error (MMSE) [6], [7] have been widely applied, as they consider local statistics.Several improved algorithms have been derived from MMSE filters, such as Kuan [8], Frost [9], and the recently improved iterative MMSE filter [10] (and its PolSAR version [11]).Later, the intensitydriven adaptive-neighborhood filter [12] was proposed, which requires good connectivity between pixels.Based on the concept of patch matching, single-temporal SAR filtering methods, such as probabilistic patch-based method [13], SAR block matching 3-D method [14], weighted sparse representation despeckle method [15], and also multitemporal SAR block-matching 3-D (MSAR-BM3D) method [16], lead to a good performance, especially in detail preservation.Innovative methods based on multichannel logarithm with Gaussian denoising were proposed in [17], [18], and [19], and deep learning methods were also adapted to SAR speckle filtering [20], [21], [22].To reduce speckle in multitemporal sets of images more efficiently, some fast strategies were introduced in [23].
Another hot topic in this domain is the use of adaptive filters based on identifying statistically homogeneous pixels (SHPs) through a region growing method, which exploits the similarity of backscatter intensity [24], [25], [26], [27], [28].In the context of differential interferometry (DInSAR) based on time series, the DespecKS method was designed for the SqueeSAR T M technique [24].Similarity between a central pixel and its surroundings is estimated by the Kolmogorov-Smirnov test [29].To overcome the computational burden of DespecKS, the fast SHP selection (FaSHPs) [25] and its variant [30] were developed.They construct a confidence interval for each pixel using the amplitude images and select SHPs according to that interval.However, all the methods mentioned here only consider the real values of the images (amplitude or intensity).In order to exploit interferometric complex information, the complex-covariancematrix-based multitemporal filter (CCM-MTF) [26] and the covariance matrix patch filter [31] were proposed.They identify SHPs by using the generalized likelihood ratio test [32] to compare InSAR CCMs of the reference pixel and its surrounding pixels.
Most of these algorithms were developed for PolSAR data or InSAR data acquired at a single date or two dates at most.In addition, a ratio-based speckle reduction approach for multitemporal PolSAR images was introduced in [45].Moreover, the GLR test on the averaged PolSAR CCM was exploited in a multitemporal polarimetric filter (MPF) [28] for the extension of differential interferometry with PolSAR data.
In general, available filters based on complex information only use the InSAR CCM and ignore the PolSAR CCM, or vice versa.To the best of the authors' knowledge, limited research exists on using interferometric and polarimetric information jointly for despeckling.One study uses the sum of Kronecker product decomposition to extract polarimetric and interferometric coherence information simultaneously [46], but it is only applicable to fully polarimetric data.The most recent study exploiting both interferometric and polarimetric information is based on BPTs [40].
In order to utilize the multitemporal InSAR and PolSAR information in a simple and efficient way, in this article, we propose a novel tensor-decomposition-based multitemporal polarimetric filter (TD-MPF).The key aspect of this filter is the generation of a new CCM, which jointly contains the information of the InSAR CCM and the PolSAR CCM for the given multichannel images.This new CCM is constructed using tensor formulation.The reason for employing tensors is that they can jointly process multi-dimensional information (spatial, spectral, and temporal) [47], [48], [49], [50].Once the tensor-based CCM is constructed for every pixel, the set of SHPs of each pixel can be identified by applying the GLR test between their CCMs, and finally, the estimated CCM is computed by averaging over the SHPs.
Both simulated and real PolSAR datasets acquired for a short period are used to verify the effectiveness of the proposed method.The filtering results are compared against MSAR-BM3D, FaSHPS, CCM-MTF, and MPF.These methods were chosen because their strategies are different and their codes are available.

II. METHODOLOGY
Throughout this article, scalars are denoted by italic lowercase and uppercase letters (e.g., m and N ), vectors are denoted by boldface lowercase letters (e.g., k), matrices are represented by bold uppercase letters (e.g., C and U), and tensors are denoted by calligraphy letters (e.g., X , Y, and C).The overall processing flow is shown in Fig. 1, in which p indicates the number of images used and m is the number of polarization channels.

A. Statistical Descriptors of Multitemporal Data 1) Quad-Pol Data:
To construct a tensor that contains the information of multitemporal quad-pol data, it is necessary to obtain two types of covariance matrix: a polarimetric covariance matrix and a single-polarization interferometric covariance matrix for each polarimetric channel.
For each pixel, a single-look complex PolSAR image acquired at date t a is represented by a scattering vector k FP,t a , which is formed under the assumption of reciprocity as follows: k FP,t a = [S HH,t a , S HV,t a , S VV,t a ] T (1) where the superscript T denotes the matrix transpose operator, FP indicates fully polarimetric, and S HV,t a is the complex backscattering coefficient measured by transmitting horizontal polarization and receiving vertical polarization, acquired at time t a .The other elements are defined similarly.The multilook of the polarimetric CCM is performed over samples at the same position along the temporal dimension.For p-dimensional temporal images, i.e., when there are p acquisition times, the PolSAR CCM can be obtained as where H denotes transpose and conjugation.The size of the resultant matrix is 3 × 3.
The formation of the interferometric CCM is inspired by the CCM-MTF method [26].A single polarization scattering vector of dimension p, denoted as k SP (where SP represents the single polarimetric channel, e.g., k VV , k HH , k HV , or k VH ), can be grouped into p/m scattering vectors k SP of m-dimension.For quad-pol data, we set m = 3; hence Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The procedure to estimate C SP for each channel in the quad-pol case is illustrated in Fig. 2. The combination of images into p/3 scattering vectors is motivated by the goal of forming polarimetric and interferometric matrices of equal dimensions, so that they can be later stacked as a tensor.2) Dual-Pol Data: The much wider availability of dual-pol data, compared to quad-pol data, from satellites such as Sentinel-1 makes it a suitable option for testing the performance of the proposed method.
For dual-pol data with VV and VH channels (which can also be altered to include other polarization channels, i.e., HH and VV, or HH and HV), the scattering vector at one pixel becomes Then, the multilook of the polarimetric CCM of size 2 × 2 is defined as As there are only two channels in dual-pol data, only two interferometric CCM can participate in the construction of the tensor, for example, C VV and C VH .To obtain these matrices, it is also necessary to group the available p-dimensional temporal images in (3) into p/m scattering vectors of m-dimension, where m equals 2 for dual-pol data: The procedure to estimate C SP for each channel in the dual-pol case is illustrated in Fig. 3.
Finally, the single-polarization interferometric matrix of each channel, with size 2 × 2, is computed according to (9).

B. Tensor Modeling of Multitemporal PolSAR Images
According to the assumption of reciprocity and using ( 2), (5), and ( 7), a total of four (C FP , C HH , C VV , C HV ) and three covariance matrices (C DP , C VV , C VH ) (for our datasets) are collected to provide a comprehensive description of one pixel with quad-pol and dual-pol multitemporal data, respectively.The detail of the construction of the tensor is illustrated in Fig. 4 for quad-pol data and for two combinations of dual-pol data (HH and VV, VV and VH).
It is important to note that, if prior knowledge of the scene is available, the weights assigned to the polarimetric and interferometric matrices can be adjusted to obtain specific results.In this study, equal weights are set for the PolSAR CCM and the InSAR CCM.According to the notation in Fig. 4, this means that w 1 = 0.5, and w 2 + w 3 + w 4 = 0.5 or w 2 + w 3 = 0.5, with w 2 = w 3 = w 4 .Since the intensity of the cross-polar channel (VH or HV) is significantly lower than that of the copolar channels (HH and VV), it is essential to compensate its weight before proceeding with the tensor modeling.For this purpose, parameter x is employed for C VH (see Fig. 4), obtained as the maximum ratio between the median intensity value of the VV (or HH) channel and the VH channel for all acquisition dates, i.e., where INT VV,t a is the intensity value of the VH channel at acquisition time t a , and similar definitions apply for the other parameters.In case of using only the copolar channels, VV and HH, the calculation of x is not required.At this stage, these matrices can be combined to create a fiveorder tensor X ∈ R F1 ×F 2 ×F 3 ×F 4 ×F 5 , where F 1 and F 2 indicate the spatial dimensions of the image, F 3 = F 4 = m, and F 5 = 4 for quad-pol data or F 5 = 3 for dual-pol data, reflecting the number of covariance matrices previously generated.

C. Alternating-Least-Squares-Based Tucker decomposition (Tucker-ALS) Algorithm for Covariance Matrix Reconstruction
The Tucker-ALS method [51] is composed of two parts: Tucker decomposition [52] and alternating least squares (ALS).
1) Tucker Decomposition: Tucker decomposition is considered as a tool for spatial compression of the modeling tensor (5)  (11) where × n is the mode-n product expressed in [53], Y ∈ R R 1 ×R 2 ×R 3 ×R 4 ×R 5 is the core tensor, and U ∈ R F n ×R n is a series of projection matrices with orthogonal columns, where n = 1, 2, . . .,5. R is the rank of the core tensor Y.This Y tensor must satisfy 2) Alternating Least Squares: The ALS part is aimed at minimizing the following least-squares cost function [54]: (13) Combining ( 12) and ( 13), the process of ALS algorithm can be described as [54,Th. 4 (14) U (n) can be obtained by using the alternating scheme.First, we optimize one of the matrix estimates in U (1) , U (2) , . . ., U (5) , while keeping the other estimates of matrices unchanged.Second, we repeat the first step until the value in ( 13) converges.The steps of Tucker-ALS 1 are summarized in Appendix A.
One advantage of Tucker-ALS is that we can predefine the rank of the output tensor, making it suitable for a subsequent GLR test.In this regard, the following rank values were preset: After dimensional reduction, the 5-D tensor X can be projected to C ∈ R F 1 ×F 2 ×m×m×1 .Therefore, when considering all pixels, we can convert

D. SHP Detection Based on the GLR Test
Under the condition of p ≥ m, the multilooked covariance matrix derived from (2), (5), and ( 7) follows a complex Wishart distribution with p degrees of freedom, where p is the number of multilooking samples, and m is the dimension of the matrix.
Thanks to the linear transform property of Tucker-ALS, the decomposed matrix C can be assumed to follow a complex Wishart distribution [32] where ), and Γ() is the function defined in [55].
If matrices for different pixels are independent and identically with the same covariance matrix Σ and the same multilooking samples p, the problem of evaluating the similarity between , p, Σ) can be rephrased as the following hypothesis test: By referring to a previous study [32], a criterion Q designed specifically for the above hypothesis test can be expressed as Large values of Q indicate that C 1 and C 2 are likely to follow a common distribution.The probability of finding a smaller value of −2ρlnQ can be simplified as The similarity between two pixels can be determined by comparing their matrices.The SHP group for the reference pixel is determined by where w is the similarity between the reference pixel and other pixels within the same search window, and R thr is the threshold.
Once the degrees of freedom m 2 and false alarm rate α are determined, R thr can be obtained by looking up the chi-square distribution table.

E. Refined Estimation of Covariance Matrices by Using SHP Group
The refined matrix Σ for a reference pixel is estimated through multilooking scattering vectors k over its SHP set Ω, as determined by (19): where N is the number of pixels in the SHP set Ω. Different kinds of matrices can be generated through utilizing the obtained SHP set.By employing the vector k in (1) or ( 6), one can obtain p polarization despeckled matrices with a size of 3 × 3 or 2 × 2. In addition, by using the vector k from (3), a filtered interferometric matrix with a size of p × p can be produced.

F. Selected Indicators for Filter Performance Evaluation
The performance of the proposed speckle filter is evaluated and compared with other filters adapted to the same type of input data.The assessment of the filter's performance is based on a set of well-defined indices.
1) Indices for Simulated Images: First, we use the wellknown equivalent number of looks [56] (ENL) index to measure the speckle strength in homogeneous areas in which INT f is the intensity value of a selected homogeneous area after filtering.
In addition, the signal-to-noise ratio (SNR) is calculated over the full image to give a global evaluation of the filtering performance SNR = 10 log 10 var(Data t ) where Data t is the true data, Data f is the data after filtering, and Data o is the noisy data.N is the number of pixels in the selected area.In this index, Data can correspond to intensity, coherence, phase, etc.A higher value of ENL and SNR indicates a better speckle reduction ability.
2) Indices for Real Images: In order to quantitatively measure the quality of filtered amplitude images, we adopted the speckle suppression index (SSI) [57] and also the more sophisticated speckle suppression and mean preservation index (SMPI) [58] for homogeneous areas, which are calculated as where subscripts f and o denote the filtered images and original noisy images, respectively.Lower values of SSI and SMPI indicate a better capability of speckle suppression and mean value preservation.
In addition, an edge enhancing index (EEI) is introduced to measure the edge preserving ability [59] where the subscript 1f or 2f denotes the filtered intensity value on either side of the edge, and 1o or 2o represents the noisy intensity value.A higher value of EEI means better edge preservation.

III. EXPERIMENTAL RESULTS
To demonstrate the effectiveness of the proposed TD-MPF, experiments were conducted using both simulated and real PolSAR datasets.The performance of TD-MPF was compared with four mainstream algorithms in terms of their capability of suppressing speckle while preserving details.These methods are FaSHPS [25], CCM-MTF [26], MPF [28], and MSAR-BM3D filter [16].
MSAR-BM3D and FaSHPS only rely on the intensity or amplitude information of multitemporal SAR images.CCM-MTF utilizes multitemporal single-polarization SAR images to build the InSAR CCM, while MPF employs multitemporal PolSAR images and computes the multilook PolSAR CCM by averaging over the time dimension.The input of TD-MPF is a combination of CCM-MTF and MPF.Table I summarizes the characteristics of the filters used in the comparison.
To ensure comparability among the filtered results produced by the five methods, the window size for FaSHPS, CCM-MTF, MPF, and TD-MPF was fixed as 15 × 15.For MSAR-BM3D, the window size was set to the default recommended values (block size = 8, stack size = 32, and search window = 39).

A. Experiment With Simulated SAR Images
For conducting a reliability analysis, the ground truth must be known.With this purpose, in the first experiment, we generated a multitemporal quad-pol dataset composed of nine single-look complex PolSAR images.The simulated scene consists of four squares with different scattering characteristics, as listed in Table II.The parameters used to describe the scattering properties include: σ, which is the intensity value of HH channel; γ and , which are the scale factor of amplitude value for VV and VH in comparison to HH respectively; ρ t , which is the complex correlation between images acquired at different time points; and ρ p , which is the complex correlation between different polarimetric    The behavior of all different methods can be better interpreted by jointly analyzing all filters' results, presented in Fig. 5, and the SHP detection maps shown in Fig. 6.Four areas of interest are marked as 1, 2, 3, and 4 in Fig. 5(a) for the subsequent analysis.Since MSAR-BM3D is not based on a similarity evaluation between pixels, its SHP map is not shown in Fig. 6.
In terms of edge preservation, VV and HH channels are sensitive to the border of Area-1, while VH channel is only good at distinguish edges of Area-4.As a result, these lines are kept well in Fig. 5, while other edges are blurred.The PolSAR filters MPF and MSAR-BM3D appear to be the best methods for retaining structures, followed by TD-MPF.Due to the changed amplitude value of VH (compared with VV and HH), it reduces the differences between new CCMs, making it difficult to preserve pixels between Area-1 and Area-2, and between Area-1 and Area-3.As a result, like in the outputs of VH channel by FaSHPS and CCM-MTF, TD-MPF also slightly blurs the edge in Area-1, Area-2, and Area-3.This phenomenon could be alleviated by modifying the weight of the VH InSAR CCM (i.e., to make it lower), which participates in the construction of the new CCM.
To evaluate the results obtained from the simulated images, quantitative measurements of SNR and ENL are used (see Table III).Besides, the correlation between different polarimetric channels ρ p (see Fig. 7) and the correlation between images acquired at different time instants ρ t (see Fig. 8) are also estimated and shown in Table III for assessment.Their definitions are given in Appendix B. In Area-1, Area-3, and Area-4, ENL values for TD-MPF are better than other methods, while in Area-2, the value is comparable to MPF.Instead of estimating it independently over the four areas, the SNR value is calculated for the entire image.The TD-MPF method shows the highest SNR values, indicating the highest level of speckle reduction.To sum up, TD-MPF offers the best compromise between detail preservation and speckle reduction among all the selected filters.

B. Experiments With Real SAR Images 1) Agricultural Area With Dual-Pol TerraSAR-X Data:
In the first experiment with real SAR data, eight dual-polarimetric (HH and VV) TerraSAR-X images covering an agricultural area in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

(c) FaSHPS(HH). (d) FaSHPS(VH). (e) CCM-MTF(VV). (f) CCM-MTF(HH). (g) CCM-MTF(VH). (h) MPF. (i) TD-MPF.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Sevilla, Spain, were used to evaluate the performance of the proposed method.These images have a spatial resolution of 6.6 m in azimuth and 2 m in range and were acquired in stripmap mode between May 25, 2017 and August 10, 2017.
The Pauli RGB composites of an image (original and filtered with the evaluated methods) are shown in Fig. 9.By visual inspection, MSAR-BM3D performs better in preserving the borders between adjacent parcels.However, the filtering performance of MSAR-BM3D inside each parcel is insufficient, as speckle is still obvious compared with other methods.On the other band, methods based on single-pol data perform well inside the parcels, but they tend to spread the edges.Apparently, the best filtering is obtained by using polarization information, i.e., by MPF and TD-MPF.
Two polarimetric features, alpha angle and entropy (H), are used to further evaluate the impact of despeckling on the discrimination among various crop types based on their scattering characteristics [60].Fig. 10 proves the edge preservation ability by plotting the H-alpha values measured in the boundary between wheat and cotton (the reference data of 2017 is provided by the official land parcel identification system).The two crop types are almost mixed together in FaSHPS(VV), FaSHPS(HH) and CCM-MTF(HH).Contrarily, TD-MPF and MPF show a comparable result in keeping edges, followed by CCM-MTF(VV).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The SHP detection map is given in Fig. 11.FaSHPS and CCM-MTF provide similar results on the border of the parcels, but CCM-MTF performs better than FaSHPS in terms of noise reduction.Fig. 11 also proves that TD-MPF and MPF are more effective in maintaining edges by selecting less pixels in areas close to different crop types.Moreover, in the homogeneous areas, the result of TD-MPF is less noisy than MPF.Therefore, TD-MPF can also achieve a better compromise between speckle suppression and detail preservation with real SAR data measured at VV and HH channels.
2) Airport Zone With Quad-Pol Radarsat-2 Data: In the second experiment, three quad-pol Radarsat-2 images acquired on January 20, 2010, February 13, 2010, and March 9, 2010 over Barcelona, Spain, were used.The beam mode is FQ9, with azimuth and range resolutions around 7.6 and 5.2 m, respectively.
The filtered results over Barcelona airport are displayed in Fig. 12. MSAR-BM3D, MPF, and TD-MPF provide clearer details at the runway and the terminal building.CCM-MTF and FaSHPS decrease the resolution to some extent.By inspecting the enlarged maps in Fig. 12, it is evident that only MPF and TD-MPF offer results comparable to the powerful MSAR-BM3D method, as the roads (black lines) are retrieved from noise.Conversely, CCM-MTF and FaSHPS lose structural information after filtering, especially in the VH channel.
The corresponding SHP identification results are shown in Fig. 13.The SHP maps in Fig. 13, obtained with only three images, demonstrate that a significant number of points in heterogeneous areas are considered as homogeneous points by FaSHPS and CCM-MTF.As a consequence, their performance is poor in heterogeneous regions.Compared with FaSHPS, CCM-MTF has better results in detect buildings since it takes advantage of  the coherence value of two adjacent images, while homogeneous areas are oversmoothed with this filter.The strong ability of detail detection in Fig. 13(g) and (h) proves the effectiveness of MPF and TD-MPF, since they can effectively exploit the information gathered in multitemporal PolSAR images.Finally, TD-MPF has less residual speckle than MPF in homogeneous areas.
The quantitative assessment is presented in Table IV.Lower SSI and SMPI indicate better results in speckle reduction and mean value preservation, respectively.MSAR-BM3D achieves the highest values for SSI and SMPI, indicating its poor performance in homogeneous areas, leading to a significant difference between the filtered image and the original one.The SSI values of CCM-MTF and FaSHPS are similar, proving their comparable noise reduction ability over the selected area.However, the SMPI value in CCM-MTF is higher than that in FaSHPS, especially in the HH channel.TD-MPF and MPF outperform other algorithms in the tested indexes.As for edge preservation, MSAR-BM3D, MPF, and TD-MPF provide similar values of EEI, indicating their comparable ability in edge preservation for real SAR data.Overall, TD-MPF produces the lowest values of SSI and SMPI but the highest value of EEI, showing its excellent performance in speckle suppression, even when processing a few images.
3) Dam Area With Dual-Pol Sentinel-1 Data: In the last experiment, 16 dual-pol (VV,VH) Sentinel-1 images acquired from May 2017 to November 2017 over the Yellow River Xiaolangdi dam of Henan province, China, were used.The observation mode is Interferometric Wide-Swath, with pixel dimensions of 13.9 and 2.3 m in azimuth and range, respectively.
The filtered amplitude images are given in Fig. 14.By visual inspection, we can see that the white line parts in FaSHPS(VH) and CCM-MTF(VH) are obviously oversmoothed.In contrast, the results from VV channel [MSAR-BM3D, FaSHPS(VV), CCM-MTF(VV)] are able to emphasize spatial structures, but still have noticeable noise in the body of the dam.The major improvement of the proposed method for this scene consists of reducing strong fluctuations while keeping important features.
By comparing the amplitude images with the map of number of identified SHP for each pixel in Fig. 15, we can see that the denoising performance of each filter comes from its SHP detection ability.FaSHPS(VV) and CCM-MTF(VV) preserve the main edge of the dam, but they cannot provide a satisfactory result in homogeneous areas due to their selection of less SHPs in these two parts.Conversely, FaSHPS(VH) and CCM-MTF(VH) select more SHPs for each pixel, leading to the spread of structures.In Fig. 15(e), obtained with the MPF method, a tradeoff performance for using VV or VH channel at the homogeneous area of the dam can be achieved, while the structure of the dam is well maintained.Finally, by combining InSAR CCM and PolSAR CCM, a better noise suppression can be achieved, while the edge of dam is still well preserved.Therefore, the proposed method is also applicable to a dataset with VV and VH channels.

IV. DISCUSSION
Experiments with simulated and real multitemporal PolSAR data have shown that InSAR CCM and PolSAR CCM exhibit different sensitivity to different targets.Therefore, using only one kind of input data, one cannot get the best result.From both visual inspection and quantitative evaluation, we concluded that the proposed TD-MPF achieves a high filtering efficiency Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

A. Construction of Tensors
The performance of the proposed method could be adjusted by setting the weight of the CCMs that participate in the construction of the tensors.The strategy in this work is to compensate the InSAR CCM of the VH channel before using it and then set InSAR CCM and PolSAR CCM in an equal weight.
We use here the Radarsat-2 data over Barcelona airport and the Sentinel-       the SHP maps presented in Fig. 16 were based on only three images with a temporal baseline of approximately one month.As a result, polarimetric information dominates in this case, and retaining more polarimetric information leads to preserve more details.However, in most cases, the good compromise for detail Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.preservation and speckle reduction is obtained when the weight for PolSAR CCM equals 0.5.

B. Suitability for Short Time Series
If only a few images are available during a study period, the lack of sufficient observations reduces the efficacy of FaSHPS and CCM-MTF.In such a case, TD-MPF still can get an accurate SHP detection map by combining time series of observations from each channel, thus benefitting edge preservation and speckle reduction.Many potential applications may benefit from the improved estimation of the identified SHP sets and the estimated CCM, for example, classification with short time series, DInSAR with short time series, etc.Since the proposed algorithm entails a multilook along the temporal dimension, it is not suitable for applications with long-time series, in which specific events and changes of the scene located at particular instants need to be preserved.e), (f), and (g)] failed to preserve the bright scatters that appeared in the first date.Based on the strategy of averaging the CCMs, CCM-MTF, MPF, and TD-MPF can capture these changes, which suggest that these three methods can detect bright targets that partially exist in the data stack along temporal dimension.

D. Computation Efficiency
Finally, we compare the execution time for the filters by using eight dual-pol images and nine quad-pol images of size 256 × 256.Details are given in Table V.The computation time The matrices Σ HH , Σ VV , and Σ VH in the diagonal are the interference matrices of size p × p.The relations of the three channels

Fig. 1 .
Fig. 1.Flowchart of the main steps performed by the proposed tensor-based method.

Fig. 2 .
Fig. 2. Procedure to estimate the initial InSAR CCM for the quad-pol case.

Fig. 3 .
Fig. 3. Procedure to estimate the initial InSAR CCM for the dual-pol case.

Fig. 4 .
Fig. 4. Scheme of the composition of tensors for one pixel and their associated weights.Left: quad-pol case.Right: dual-pol case with two different combinations of channels.
1 data over Xiaolangdi dam to demonstrate the importance of the VH compensation, as well as the equal weighting of InSAR and PolSAR CCMs.The corresponding maps of number of SHPs for different cases are shown in Figs.16 and 17, respectively.Results show that without VH compensation [see Fig. 16(b), (d), and (f) and Fig. 17 (b), (d), and (f)], the SHP maps exhibit differences compared to Fig.16(a), (c), and (e) and Fig.17 (a), (c), and (e).This disparity is particularly evident when the weight assigned to PolSAR CCM (w 1 ) is low, for example, in the main part of the dam (see Fig.

17
) and the homogeneous area around the airport (see Fig.16), the results with VH compensation step are more smooth.It is also important to choose a reasonable weight for PolSAR CCM.Once the weight of PolSAR CCM is fixed, for example w 1 = 0.2, then the weight for each InSAR CCM equals (1 − 0.2)/3 for the quad-pol case and (1 − 0.2)/2 for the dual-pol case.When the weight for PolSAR CCM equals 0.8, the SHP detection maps will be similar to those for MPF [see Figs.13(g) and 15(e)].If the weight for PolSAR CCM equals 0.2, the result is close to that of CCM-MTF [see Fig. 13(d), (e), and (f) and Fig. 15(c) and (d)].It is worth mentioning that

Fig. 18
Fig. 18 is provided to analyze the filters' performance when target distribution changes.In the first acquired image (January 20, 2010), bright targets were observed near the terminal, while in the second (February 13, 2010) and third (March 9, 2010) images, the number of bright targets decreased.The results generated by the different filters for the first image are shown in Fig. 18(d)-(h) and (l).FaSHPs [see Fig.18(e), (f), and (g)] failed to preserve the bright scatters that appeared in the first date.Based on the strategy of averaging the CCMs, CCM-MTF, MPF, and TD-MPF can capture these changes, which suggest that these three methods can detect bright targets that partially exist in the data stack along temporal dimension.

18
Fig. 18 is provided to analyze the filters' performance when target distribution changes.In the first acquired image (January 20, 2010), bright targets were observed near the terminal, while in the second (February 13, 2010) and third (March 9, 2010) images, the number of bright targets decreased.The results generated by the different filters for the first image are shown in Fig. 18(d)-(h) and (l).FaSHPs [see Fig.18(e), (f), and (g)] failed to preserve the bright scatters that appeared in the first date.Based on the strategy of averaging the CCMs, CCM-MTF, MPF, and TD-MPF can capture these changes, which suggest that these three methods can detect bright targets that partially exist in the data stack along temporal dimension.
Fig. 18 is provided to analyze the filters' performance when target distribution changes.In the first acquired image (January 20, 2010), bright targets were observed near the terminal, while in the second (February 13, 2010) and third (March 9, 2010) images, the number of bright targets decreased.The results generated by the different filters for the first image are shown in Fig. 18(d)-(h) and (l).FaSHPs [see Fig.18(e), (f), and (g)] failed to preserve the bright scatters that appeared in the first date.Based on the strategy of averaging the CCMs, CCM-MTF, MPF, and TD-MPF can capture these changes, which suggest that these three methods can detect bright targets that partially exist in the data stack along temporal dimension.
of MSAR-BM3D and FaSHPS depends on both image size and temporal dimension.Only dual-pol case is shown in MSAR-BM3D because it can only deal with dataset with temporal dimension 2 n .The time consumption for CCM-MTF, MPF and TD-MPF relays on image size and matrix size; as a result, quadpol data 256 × 256 × 3 × 3 require longer processing time than dual-pol data 256 × 256 × 2 × 2. Because of the special design of FaSHPS, it results as the fastest method.For the proposed methods, as only one tensor decomposition step is added before the similarity test, the computational cost of TD-MPF is similar to that of other filters based on the GLR test.The execution time is measured on a workstation equipped with an Intel Xeon CPU at 3.50 GHz and 16 GB of RAM.V. CONCLUSIONIn this article, a novel multitemporal PolSAR filtering method named TD-MPF was proposed by integrating InSAR CCM and PolSAR CCM through tensor decomposition.The proposed method was compared with different filters, including MSAR-BM3D, FaSHPS, CCM-MTF, and MPF.Experiments using multitemporal simulated and real PolSAR images demonstrated the robust performance of TD-MPF at different scenarios for edge preservation and speckle suppression.The proposed TD-MPF is suitable for processing both quad-pol and dual-pol SAR data stacks.As a result, the improved estimation of InSAR CCM or PolSAR CCM may benefit a number of potential applications, such as deformation or crop monitoring, classification, and change detection.APPENDIXA PSEUDOCODE OF TUCKER-ALS APPENDIX B SIMULATION OF MULTITEMPORAL POLSAR IMAGESThe simulated dataset has been generated according to the complex Gaussian polarimetric model, for each pixel in the p-dimensional multitemporal PolSAR dataset, assuming a reflection symmetric target, which has a covariance matrix C of size 3p × 3p.From a matrix C, we can obtain p SLC SAR images for each polarimetric channel using the procedure introduced in[61]

TABLE I CHARACTERISTICS
OF THE USED FILTERS

TABLE II DIFFERENT
SIMULATED VALUES FOR FOUR AREAS channels.More details about the simulation procedure are given in Appendix B.

TABLE III QUANTITATIVE
EVALUATION RESULTS FOR SIMULATED IMAGES

TABLE IV STATISTICAL
INDICES OF THE SELECTED HOMOGENEOUS REGION IN THE FILTERED IMAGES CORRESPONDING TO THE AIRPORT AREA WITH QUAD-POL RADARSAT-2 DATA

TABLE V COMPARISON
OF THE COMPUTATION TIME FOR WINDOW SIZE 15 × 15 OVER 256 × 256 PIXELS