2.1.

Basic of Tensor Algebra

A tensor is defined as a multidimensional array, and its number of dimensions is referred to as mode or order. The tensor definition unifies the framework for depicting the data. Scalar, vector, and matrix are special cases of the tensor, with the modes of 0, 1, and 2, respectively. Scalars are denoted by lowercase letters, e.g., $x$; vectors are denoted by boldface lowercase letters, e.g., $\mathbf{x}$; matrices by boldface capital letters, e.g., $\mathbf{X}$. For simplicity, only the third-order tensor is considered throughout this paper, which is denoted by Euler script letters, i.e., $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$, whose single entry is a real scalar and expressed as $x_{i_1{i}_2{i}_3}$, where $1\le i_n\le I_n(1\le n\le 3)$ are indexes along each mode. Similar to the Frobenius norm of a matrix, the norm of a tensor $\mathcal{X}$ is the inner product with itself, i.e., ${\Vert \mathcal{X}\Vert }_F=\sqrt{\sum _{i_1=1}^{I_1}\sum _{i_2=1}^{I_2}\sum _{i_3=1}^{I_3}{x}_{i_1{i}_2{i}_3}^2}$.

In terms of a tensor, fibers are the higher-order analogue of matrix rows and columns. The mode-$n$ fibers are the $I_n$-dimensional vectors obtained by varying index $i_n$, while keeping other indexes fixed. Matricization or unfolding is defined as the process of reshaping a tensor into a matrix. The mode-$n$ matricization of a tensor $\mathcal{X}$ arranges the mode-$n$ fibers as the columns of the resulting matrix, which is denoted by ${\mathbf{X}}_{(n)}$. For instance, the mode-1 matricization process transforms the original tensor $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ into a matrix ${\mathbf{X}}_{(1)}\in {\mathbb{R}}^{I_1\times I_2{I}_3}$, whose columns consist of $I_2{I}_3$ mode-1 fibers, with each of size $I_1\times 1$. Figure 1 illustrates the mode-$n$ matricization process of a third-order tensor.

Fig. 1

Illustration of the mode-$n$ fibers and its corresponding matricization process.

The mode- $n$ product of a tensor $\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ with a real matrix $\mathbf{A}\in {\mathbb{R}}^{J\times I_n}$ is defined as the multiplication of mode- $n$ unfolding matrix ${\mathbf{X}}_{(n)}$ and $\mathbf{A}$, which is expressed as Ref. 15

The above definitions can be easily extended to higher-mode tensor and more mathematical details could be found in Ref. 15.

2.2.

Tensor Representation of PolSAR Image Pixels

Under the assumption of monostatic case and medium reciprocity, PolSAR data are expressed as multilook complex data by spatial averaging¹

Based on the covariance matrix, several parameters can be directly calculated such as its diagonal elements and correlation coefficients, which provide a good understanding of the relationships among different polarimetric channels. In addition, many effective target decomposition methods have been developed over the past few decades,¹ including the Pauli,¹⁸ Freeman,^19,20 Krogager,²¹ Van Zyl,²² and Yamaguchi²³ methods. These decomposition theorems provide different polarimetric descriptors that could give an interpretation of canonical scattering mechanisms. Each descriptor has its own strength and weaknesses for discriminating different terrain types. Some literatures have employed different combinations of multiple features, whose experimental results indicate the significant improvement of the classification accuracy under certain experimental conditions.^7–10

Therefore, multifeatures combination is investigated in this paper. The considered features in Table 1 are selected as Qi et al. did in Ref. 24, which are generally selected in a rather random manner, without further careful consideration. Totally $D=21$ polarimetric features are calculated as listed in Table 1. All these features are calculated using PolSARpro software and the notations are consistent for better understanding.²⁵ By stacking all the extracted features, the PolSAR image could be characterized by a third-order PolSAR feature tensor $\mathcal{F}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$, as shown in Fig. 2. This representation provides a global description for the PolSAR data. By employing a sliding window with the size of $W\times W$ on the feature tensor $\mathcal{F}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$, each pixel could be expressed as a third-mode pixel subtensor $\mathcal{P}\in {\mathbb{R}}^{W\times W\times D}$, where the first two modes represent the neighboring pixels, and the third represents the polarimetric features. This tensor representation provides a local neighborhood description for PolSAR data.

Table 1

Extracted polarimetric features.

Fig. 2

Global and local tensor representation of a pixel in PolSAR data.

Under traditional matrix algebra framework, the PolSAR feature tensor $\mathcal{F}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ is reshaped into a feature matrix ${\mathbf{F}}_{(3)}\in {\mathbb{R}}^{I_3\times I_1{I}_2}$, with each column represents the feature vector of a pixel $f\in {\mathbb{R}}^{I_3}$. It is worth noting that this process is equivalent to the mode-3 matricization of the PolSAR feature tensor. However, this matrix representation uses only the 21 polarimetric attributes of one pixel, which ignores the spatial correlation of scattering coefficients among neighboring pixels. By employing a sliding window with the size of $W\times W$ on the feature tensor $\mathcal{F}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$, each pixel could be expressed as a third-mode pixel subtensor $\mathcal{P}\in {\mathbb{R}}^{W\times W\times I_3}$, where the first two modes represent the neighborhood pixels, and the third represents the polarimetric features. Inclusion of polarization diversity captures abundant structural and textural information of the medium. On the contrary, if we only consider spatial correlation without considering polarimetric information, the classification accuracy is rather poor. This is consistent with the fact that single-polarization SAR data are not good candidates for land cover classification.²⁶

Unlike traditional matrix-based representation, the proposed single-pixel tensor representation $\mathcal{P}\in {\mathbb{R}}^{W\times W\times I_3}$ not only combines multiple features, but also takes the spatial homogeneity into consideration. This local tensor representation conserves the natural data structure and would be beneficial for analysis of local scattering mechanisms, and is the model foundation for subsequent processing.

However, this tensor representation with high dimensionality requires larger memory storage, and thus increases computation complexity. Meanwhile, the redundancy existed in both spatial and feature dimension may pose a hindrance for accurate classification. Therefore, it is necessary to find a tensor subspace that captures most of the variation in the original tensor object, and thus improves the accuracy and efficiency of the classification scheme. In the following, we will provide solutions for this issue.

3.1.

Multilinear Principal Component Analysis

The goal of MPCA is to reduce the dimensionality of a tensor consisting of a large number of interrelated variables, while retaining as much as possible the variation present in the original tensor. Given a training set $\{{\mathcal{P}}_m\in {\mathbb{R}}^{I_1\times I_2\times I_3},\hspace{0.17em}m=1,\dots ,M\}$ consists of $M$ tensors of size $I_1\times I_2\times I_3$, and a set of projections matrices $\{{\mathbf{U}}^{(n)}\in {\mathbb{R}}^{I_n\times J_n},n=\mathrm{1,2},3\}$, the projected subtensors $\{{\mathcal{Y}}_m\in {\mathbb{R}}^{J_1\times J_2\times J_3},\hspace{0.17em}m=1,\dots ,M\}$ are expressed as

Equivalently, Eq. (4) can be rewritten in the form of mode-$n$ matricization

MPCA seeks to find the best projection matrices ${\mathbf{U}}^{(n)}$ such that the projected subtensors have maximum energy, which is essentially an optimization problem¹⁶

Based on Eq. (4), the energy of projected tensors ${\mathrm{\Psi }}_{\mathcal{Y}}$ can be expanded as

Similarly, based on Eqs. (6) and (7), ${\mathrm{\Psi }}_{\mathcal{Y}}$ in Eq. (9) also can be rewritten in terms of the mode-2 and mode-3 matricizations:

Therefore, the original optimization problem can be formulated as three-equivalent optimization problems

The formulation of the optimization problem in Eq. (13) is the same with that of the matrix PCA.²⁷ Therefore, the optimal solution to ${\mathbf{U}}^{(n)}$ is obtained by applying the eigenvalue decomposition on ${\mathbf{\varphi }}^{(n)}$ and assigning the eigenvectors corresponding to the largest $J_n$ eigenvalues as the columns of ${\mathbf{U}}^{(n)}$. However, the optimal solution to ${\mathbf{U}}^{(n)}$ depends on other projection matrices, and it is rather difficult to solve all the projection matrices simultaneously. Therefore, an alternative least square (ALS) scheme is applied to iteratively solve the projection matrices. The projection matrices ${\mathbf{U}}^{(n)}$ define the mapping from the original high-dimensional training tensor into an intrinsic low-dimensional tensor subspace. This tensor subspace is assumed to capture most of the variation in the training tensors set and would benefit for the classification process.

3.2.

Comparison with Conventional Matrix-Based Methods

Under the matrix algebra framework, traditional matrix-based DR methods such as PCA,²⁷ independent component analysis (ICA),²⁸ factor analysis (FA),²⁹ and linear discriminate analysis (LDA)²⁹ require to reshape the third-order PolSAR feature tensor $\mathcal{F}\in {\mathbb{R}}^{I_1\times I_2\times I_3}$ to the two-dimensional feature matrix ${\mathbf{F}}_{(3)}\in {\mathbb{R}}^{I_3\times I_1{I}_2}$, as shown in Fig. 2. Then the feature matrix is multiplied by a transformation matrix so that a reduced feature matrix is obtained. The generative model of the linear matrix-based DR method can be expressed as

However, the projection in Eq. (14) does not take the spatial relations among neighboring pixels into consideration. It is worth noting that the reduced feature matrix $\mathbf{Y}$ can also be expressed as the mode-3 matricization of a third-order tensor. Comparing Eq. (4) with Eq. (14), we can easily tell the difference between the matrix-based and tensor-based learning techniques. For matrix-based DR methods, only mode-3 feature information is utilized, while ignores the spatial information in mode-1 and mode-2. For the proposed tensor-based scheme, the pixel subtensor $\mathcal{P}\in {\mathbb{R}}^{W\times W\times I_3}$ is processed directly without further reshaping process. Therefore, a dimensionality reduction algorithm operating directly on a tensor object rather than its vectorized version is desirable. Owing to the iterative ALS procedure, the proposed scheme takes into account the cross dependence between each mode, which means any projection along a given mode depends on the projections along all other modes. The projection matrix for feature extraction is estimated based on the update of the spatial-mode information. In summary, the main differences between the proposed tensor-based method with other matrix-based methods reside in two aspects: data representation and projection matrix optimization. This explains how the spatial information is used in the proposed tensor representation and learning scheme.

3.3.

Flowchart of Classification Scheme

Figure 3 summarizes the flowchart of the proposed supervised classification scheme. After obtaining the tensor representation of all the pixels as shown in Sec. 2, a set of training samples and a set of test samples are determined according to the prior knowledge of the ground truth. Then the tensor learning technique, i.e., MPCA, is applied on the training tensors set to explore the data structure and find projection matrices along each mode. Further, all the test tensor samples undergo the same mapping with the training samples via the projection matrices. The projected training samples are utilized to train a specific classifier, such as the KNN classifier.¹⁷ Subsequently, the projected test samples are fed to the trained classifier, and the classification result is obtained.

Fig. 3

Flowchart of supervised classification scheme based on tensor representation and learning.

Indeed, advanced classifiers could improve the classification performance to some extent. However, without well-selected discriminative features, it would be difficult to obtain high classification accuracy, even if using a rather complex classifier. Otherwise, a simple classifier still can obtain a good classification result if well-separable features are given. Since our work mainly focuses on the DR process, the rationality behind the choice of KNN is to show that the extracted features by the proposed method are separable so that even a simple classifier can achieve a satisfactory result.

4.1.

Results of Simulated Data

For PolSAR data, it is common to assume that the multilook covariance matrix is complex Wishart distributed, whose probability density function is expressed as Ref. 2

In order to well-characterize heterogeneity of the scattering medium, a more refined non-Gaussian product model that incorporating the texture variation and speckle is defined as Refs. 5 and 6:

Based on Eq. (16), if given the mean covariance matrices $\mathrm{\Sigma }$, the number of looks $L$ and the texture parameter $t$, one can simulate different realizations of the PolSAR covariance matrix.¹ In this experiment, four-look fully polarimetric data are simulated for illustration. It is $800\times 1000$ pixels in size and has seven types of land covers. The mean covariance matrices $\mathrm{\Sigma }$ of each type are extracted from the real datasets of AIRSAR Flevoland dataset.³⁰

Figure 4(a) shows the color-coded Pauli image of the simulated data with different land-cover types and their corresponding class labels. The number of samples for each class is listed in Table 2. All the classes are supposed to follow Gaussian distributed, except that the class 3 is more non-Gaussian distributed considering the graininess-like texture variation. This can be verified from Fig. 4(b), which illustrates the non-Gaussianity of the simulated data using the relative kurtosis as defined in Ref. 5. A larger relative kurtosis value indicates a relatively larger degree of non-Gaussianity. The presence of several straight lines results from the heterogeneous boundary between different classes.

Fig. 4

Simulated data with corresponding class labels: (a) color-coded Pauli image of the original simulated data and (b) illustration of non-Gaussianity using relative kurtosis.

Table 2

Samples for each class in simulated data.

According to Table 1, a total of 21 polarimetric features are extracted. To balance the contributions of each feature, all the features are scaled in the range of 0 to 1 with respect to their maximum values. Then each pixel is formulated as a third-order tensor $\mathcal{P}\in {\mathbb{R}}^{W\times W\times 21}$ as illustrated in Fig. 2. In this simulation, the sliding neighboring window size is set as $W=5$, and the dimension of projected features is set as $p=3$ for illustration. Only 1% of the samples is selected as training samples as shown in Table 2.

First, several classical matrix-based linear DR methods, i.e., PCA,²⁷ ICA,²⁸ FA,²⁹ and LDA,²⁹ are applied to obtain a reduced feature set, respectively. It is worth noting that nonlinear methods including ISOMAP, LPP, and Laplacian Eigenmaps require constructing neighborhood graph, which has extensive computation and is very time consuming. Therefore, these nonlinear techniques are not considered for comparison here.

The scatter plots of the reduced feature set by different matrix-based DR methods are illustrated in Figs. 5(a)–5(d). Each dot represents a training sample after DR, and different colors indicate different land-cover types. From Figs. 5(a)–5(c), it is shown that the unsupervised DR methods, i.e., PCA, ICA, and FA, fail to distinguish the samples belonging to different classes, and almost all the samples are mixed together in the feature space. Unsupervised learning cannot properly model underlying structures and characteristics of different classes. LDA is the supervised approach to learn discriminant subspaces by utilizing the a priori label information, where the between-class scatter of samples is maximized and the within-class scatter is minimized at the same time. Hence the samples are more discriminative in Fig. 5(d) than that in Figs. 5(a)–5(c). However, the scatter plots in Figs. 5(a)–5(d) indicate that the projected features are not distinguishable enough, which would pose hindrance for accurate classification even with a rather complex classifier.

Fig. 5

Scatter plots of extracted intrinsic features obtained by different methods: (a) PCA, (b) ICA, (c) FA, (d) LDA, and (e) MPCA.

Further, the proposed tensor-learning technique, i.e., MPCA, is applied and its resulting scatter plot is presented in Fig. 5(e). It is shown that each class is much more concentrated and different classes are more distinguishable, i.e., the samples have low within-class variation and large between-class discriminability. The proposed method extracts the most salient information embedded in the local tensor object by exploiting the mutual correlation between the spatial and feature dimension. This indicates that the redundancy among features is minimized and a simpler classifier should have rather high classification accuracy.

First, the complex Wishart classifier (CWC)² is applied as an evaluation benchmark. The CWC only utilizes the covariance matrix, which classifies the PolSAR data by exploiting statistical properties under the Bayesian maximum likelihood rule.² A particular result of CWC is shown in Fig. 6(a). It is shown that the CWC suffers from the adverse impact of speckle, and misclassifies many samples in a homogeneous area, especially in class 3.

Fig. 6

A particular classification result using different methods: (a) CWC, (b) PCA + KNN, (c) ICA + KNN, (d) FA + KNN, (e) LDA + KNN, and (f) proposed: MPCA + KNN.

Then after applying various DR techniques, i.e., PCA, ICA, FA, LDA, and proposed MPCA, one can obtain a low-dimensional intrinsic feature set, respectively. Further, these samples are fed to the KNN classifier as illustrated in Sec. 3.3. For the KNN classifier, the number of nearest neighbors $k$ is chosen by employing fivefold cross validation.

Figures 6(b)–6(e) shows a particular classification result for PCA + KNN, ICA + KNN, FA + KNN, and LDA + KNN, respectively. It is shown that all these matrix-based DR methods lead to very messy classification results because they fail to obtain distinguishable feature sets and suffers from the influence of speckle noise. Under the situation with extremely limited training samples, i.e., only 1% of all the samples, the conventional matrix-based mapping has large generalization error. Although the spatial information is considered, but without well extracted, the redundancy may pose an adverse impact on the classification accuracy.

Figure 6(f) shows a particular classification result for MPCA + KNN. It is shown that many more pixels are correctly classified as a class in the homogeneous area. The tensor representation and learning copes with the problem of insufficient training samples, which is more realistic in practical applications.

Figure 7 compares the overall accuracy (OA) values for different combinations of DR method and KNN classifier, which is calculated over 100 independent Monte Carlo runs. The dimension of projected features is set as $p=3$, neighborhood window size $W=5$, and ratio of the training samples $\rho =1\%$. These matrix-based methods suffer from the curse of dimensionality, i.e., it is difficult to obtain high OA accuracy with a small number of training samples available. PCA + KNN has the worst performance. The very large deviations of OAs for ICA and LDA indicate their poor generalization performance to the test samples. This phenomenon attributes to the fact that they are not deterministic methods and the projections are highly dependent on the training samples.²⁸ The proposed MPCA has the highest OA values, which indicate its projected features are more discriminative. Even compared with the CWC, the accuracy is improved as much as $\sim 12\%$. in addition, its robustness is proved by rather small deviation of OAs.

Fig. 7

Comparison of OA values for different combinations of DR method and KNN classifier.

The dimensionality of the projected feature $p$, the ratio of training samples $\rho$, and the window size $W$ are critical parameters for the DR process. To provide a thorough analysis on the performance of the proposed tensor learning technique, some results on different experimental conditions are presented in the following.

4.1.1.

Influence of projected feature dimension

The dimensionality of the projected feature defines a feature space that characterizes the data samples. High-dimensional features would increase the computation burden, which leads to an inefficient classification process. Moreover, for a classifier, a high-dimensional space requires much more training samples to obtain a better classification performance. To provide a comprehensive evaluation, Fig. 8 presents the variation of mean OAs under different dimensionalities of projected features. It is shown that the mean OA increases with the dimensionality of projected features at first and becomes almost stable subsequently. This proves the validity of the assumption that the significant information of the original higher dimensional is embedded in a low-dimensional feature space. The proposed MPCA-based method has a superior performance.

Fig. 8

Simulated data: variation of OA with various dimensions of projected features.

4.1.2.

Influence of training samples size

The training samples are representatives of the data and are used to determine the structure of the classifier. The variability of the training data number has an influence on the classification accuracy. In practical application, the size of the training samples may be limited. Therefore, the performance of the proposed scheme is investigated for different ratios of training samples ranging from 0.1% to 20%. In this comparison, the projected feature dimension is set as $p=3$ and the neighboring window is set as $W=5$.

Figures 9(a)–9(d) shows the performance of the conventional matrix-based methods. It is shown that these methods perform poorly when the raining samples are very limited. This attributes to the fact that they suffer from the curse of dimensionality because of the high-dimensional nature of the feature vector. As more training samples are available, these methods have significant performance improvement. This plot highlights the importance of labeled training samples.

Fig. 9

Simulated data: variation of OA with various percentages of training samples: (a) PCA, (b) ICA, (c) FA, (d) FA, and (e) proposed: MPCA.

As shown in Fig. 9(e), with the number of training points increase, the mean value curve of OAs has a trend of increasing with no $>1\%$. Meanwhile, the deviation decreases because more training samples would mitigate the ambiguity of the classifier. Even when the number of training samples is rather small, i.e., 0.1%, the OA of proposed method is close to 98%, which is superior to the conventional matrix-based methods with 90% training samples available. This demonstrates that the features learned by the proposed scheme are of high generality. The increase of training samples would complicate the training process and slow down the entire classification step. For the proposed scheme, a small training set is able to obtain rather accurate classification result, which could accelerate the training phase and improve efficiency. The proposed method could cope with the curse of dimensionality and could exploit the available training dataset more effectively.

4.1.3.

Influence of neighboring window size

Another important adjustable parameter is the window size that used for tensor representation. In this comparison, the projected feature dimension is set as $p=3$ and the ratio of training samples is set as $\rho =1\%$. Figure 10 shows the OAs under various window sizes of the proposed tensor-based methods. It is shown that OA increases with the enlargement of neighboring window size. This phenomenon demonstrates that the tensor representation plays a significant role in improving the classification performance by considering more spatial information. However, a larger window requires more memory storage and also more computation time, which definitely reduces the efficiency. Moreover, a larger window may not be suitable for analyzing scattering areas with high heterogeneity. Therefore, there is a trade-off between efficiency and accuracy.

Fig. 10

Simulated data: variation of OA with various window sizes.

From the above analysis, it is shown that the conventional matrix-based methods are very sensitive to the ratio of training samples, and the projected feature dimension and neighboring window size have little impact. Under the situation without sufficient training samples, the proposed method is capable of coping with the curse of dimensionality and still has superior generalization performance. The proposed tensor-based classification scheme is less sensitive to the tuning parameters than other methods and has a superior performance under the same parameter configuration.

4.2.

Results of Real EMISAR Data

In this part, a real dataset is provided to test the validity of the proposed scheme. This dataset is acquired over the Foulum area by Danish EMISAR.³¹ A subarea is extracted from this dataset with a size of $286\times 337$ pixels. Figure 11(a) shows the color-coded Pauli image of the area. Based on Ref. 5, we manually specify nine types of terrain in this area and the corresponding ground-truth map is shown in Fig. 11(b).

Fig. 11

EMISAR L-band fully polarimetric data of Foulum, Netherlands: (a) color-coded Pauli image and (b) ground-truth map based on Ref. 5.

Next, the real Foulum dataset undergoes the same processing flow with the simulated data, as illustrated in Sec. 3.3. Each pixel is represented as a third-mode tensor $\mathcal{P}\in {\mathbb{R}}^{5\times 5\times 21}$, and a training set is determined based on a priori ground-truth map. Then different DR methods are applied on the training set to obtain an intrinsic low-dimensional feature set. Furthermore, the projected intrinsic features are proceeding to the KNN classifier for supervised classification. The experimental settings are the same with the simulated data in Sec.4.1.

Figure 12(a) presents the classification result of CWC. Many misclassification points occur in the homogeneous areas such as type 1 (i.e., right bottom corner) due to the influence of speckle noise. Figures 12(b)–12(e) illustrate the results after applying the PCA + KNN, ICA + KNN, FA + KNN, and LDA + KNN, respectively. It is shown that PCA + KNN has the worst performance, even type 6 and type 7 are mixed together. The LDA + KNN has a comparatively better performance than other methods, because it is a supervised DR method. All these methods including CWC seem to suffer from the adverse impact of speckle noise. Figure 12(f) shows the result after realizing the proposed MPCA + KNN. It is shown that the classification result is much smoother in homogeneous areas such as type 1. It implies that the introduction of spatial neighboring information could alleviate the speckle noise to some extent, which facilitates the classification and improves the accuracy. Further, the OA values are calculated over 100 independent Monte Carlo runs, as shown in Fig. 13. It is worth noting that only 10% of samples in the ground-truth map are used for training, and all other samples are used to evaluate the OA. Since the total samples of simulated data are almost 9 times than the real data, and thus here the training ratio is set as 10%, not set 1% as in simulated data. According to Fig. 13, it is shown that the proposed scheme has the best performance with a rather large lead in the OA values.

Fig. 12

Real EMISAR data: a particular classification result using different methods: (a) CWC, (b) PCA + KNN, (c) ICA + KNN, (d) FA + KNN, (e) LDA + KNN, and (f) proposed: MPCA + KNN.

Fig. 13

Real EMISAR data: comparison of OA values for KNN classifier.

Fig. 14

Real EMISAR data: variation of OA with various dimensions of projected features for KNN classifier.

Fig. 15

Real EMISAR data: variation of OA with various percentages of training samples: (a) LDA, (b) ICA, (c) PCA, (d) FA, and (e) proposed: MPCA.

Fig. 16

Real EMISAR data: variation of OA with various window sizes for the proposed method.

Similarly, as in Sec. 4.1, some further performance discussions are provided regarding with the influence of the projected dimension, ratio of training samples, and window size. Figure 14 plots the variation of OAs with various projected dimensions. Figure 15 shows the variation of OAs for different methods with various percentages of training samples. Figure 16 shows the variation of proposed method with various window sizes. Similar observations with the simulated dataset can be drawn from these figures.

4.3.

Results of Real AIRSAR Data

In order to further evaluate the effectiveness of the proposed data, other real data are also tested. This dataset is acquired over Flevoland in Netherlands by the National Aeronautics and Space Administration (NASA)/Jet Propulsion Laboratory (JPL) AIRSAR on August 16, 1989. The PolSAR image used here has a size of $270\times 250$ pixels. Figure 17(a) shows the color-coded Pauli image of the area. The PolSAR image is impaired with speckle noise, which poses a hindrance for image interpretation. Based on Ref. 32, there are totally six terrain types in this area, and the corresponding ground-truth map is shown in Fig. 17 (b).

Fig. 17

AIRSAR L-band fully polarimetric data of Flevoland, Netherlands: (a) Color-coded Pauli image and (b) ground-truth map.

In Secs. 4.1 and 4.2, the experiments of simulated data and real EMISAR data are conducted for nonfiltered PolSAR data. In this part, in order to evaluate the effects of filtering on the classification performance, the real AIRSAR data are filtered using Lee refined filter.³³ The projected feature dimension is set as $p=3$, and the ratio of training samples is set as 10%.

Fig. 18

Real AIRSAR data: a particular classification result using different methods: (a) CWC, (b) PCA + KNN, (c) ICA + KNN, (d) FA + KNN, (e) LDA + KNN, and (f) proposed: MPCA + KNN.

Fig. 19

Real AIRSAR data: variation of OA with various dimensions of projected features for KNN classifier.

Fig. 20

Real AIRSAR data: variation of OA with various percentages of training samples: (a) LDA, (b) ICA, (c) PCA, (d) FA, and (e) proposed: MPCA.

Fig. 21

Real AIRSAR data: variation of OA with various window sizes for the proposed method.

Figure 18 compares the classification results of using existing matrix-based methods and proposed methods. Figures 19 and 20 provide detailed quantitative performance comparisons in terms of the projected feature dimension and training samples size. Meanwhile, Fig. 21 plots the variation of overall classification accuracy with the neighboring window size. From these figures, it is shown that after filtering preprocessing, the OA of proposed method is still the best, even better than the benchmark CWC for nearly 4%. This also indicates the superiority of the proposed tensor-based processing scheme.

4.4.

Discussion on Computation Time

Table 3 compares the computational time for each DR method for the simulated data and real data. In this comparison, the projected dimension is fixed as $p=3$, the ratio of training samples is $\rho =10\%$, and the neighboring window size is $W=5$. This comparison is realized in a personal laptop with Intel^® Core i7-6820 using MATLAB R2013a software. The time is averaged under 100 independent runs. It is shown that matrix-based methods are quite efficient, the proposed tensor-based method costs much more time than these matrix-based methods. This is reasonable because the proposed method involves alternative optimizing along the spatial and feature dimension, whereas the matrix-based methods only processing along the feature dimension. The number of samples in the simulated data is larger than that of the real data, and thus its average time is much higher. The solving process of FA involves an iterative expectation–maximization–optimization under maximum likelihood criterion. As the number of samples increases, more iterations are needed for convergence. Therefore, the cost time of FA value is much higher than normal trend. Although the samples of simulated data are $\sim 9$ times the real data, the matrix-based method still is quite efficient. The computation time of the proposed technique grows significantly with the increase of processing samples. However, the extra time cost is worthy considering the great improvement of classification accuracy.

Table 3

Comparison of computation time for each DR method.