Multitemporal Polarimetric SAR Change Detection for Crop Monitoring and Crop Type Classiﬁcation

—The interpretation of multidimensional synthetic aperture radar (SAR) data often requires expert knowledge. In fact, it requires to simultaneously consider several time series of polarimetric features to understand the physical changes of a target and its temporal evolution. In an effort to characterize the changes over time, multitemporal polarimetric SAR (MTPolSAR)

has become a key tool in remote sensing. The day-and-night and all-weather capabilities of SAR satellites offer significant opportunities for continuous and remote monitoring regardless of cloud cover. Additionally, the information extracted from the acquired SAR data can be augmented by means of polarimetric analysis, exploiting the covariance between the different polarization channels measured by the satellite [1], [2]. However, the process of deriving information is often not intuitive due to the large amount of polarimetric features that change over time driven by natural or human-made processes. This analysis is essential for evaluating changes in earth processes, land cover and use, including applications such as forest [3], agriculture [4], [5], glacier changes [6], and urban settlement monitoring, among others. Agricultural fields are just one of the examples of such dynamical behavior since the morphological parameters of the plants evolve along the season varying the PolSAR response accordingly. The most widely used polarimetric analyses consider the simultaneous interpretation of several polarimetric time series, resulting from multitemporal backscatter measurements and PolSAR decompositions [4], [7]. However, the process of identifying the adequate PolSAR decomposition to use the set of features that are best suited for the problem at hand and how to combine them to identify informative patterns is not straightforward.
Other methods focus on identifying the scattering mechanism (SM) that are suffering the largest change by means of PolSAR change detection. Test statistic methods such as the likelihood ratio test focus on evaluating the equality of two covariance matrices that follow the complex Wishart distribution [8] [9]. Likewise, the complex Hotelling-Lawley trace test statistic is used in [10] to evaluate the similarity between two complex covariance matrices which is subsequently thresholded to determine changes between images. Statistical information theory has also been applied to measure divergence between complex covariance matrices [11] as well as distance-based methods which evaluate the polarimetric similarities between images to identify changed regions [12], [13]. Note that although these methods exploit the PolSAR information to quantify the intensity of changes, they do not focus on identifying the type of changes between the images. Therefore, they are not suitable for analyzing the SM dynamics and its temporal evolution. In [14], a method was presented in order to detect polarimetric changes between two images. However, a way to combine this information with the changes in the backscatter amplitude is not directly embedded in the detection process. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ In terms of multitemporal change detection, the likelihood ratio test statistic was extended [15], [16] in order to handle multiple images and test the equality of several covariance matrices. In this case, the authors focused on detecting whether change occurred in any of the images in the stack and try to estimate when it happened. Note, however, that the interpretation regarding the type of change which is vital for change analysis is not considered in this approach. In [17], the same authors propose to improve this method by adding the capability of identifying the change direction. This is done by computing the definiteness (positive or negative definiteness) of the difference of covariance matrices that model the PolSAR response in the two dates considered for the change detection. This way the method detects not only that there was a change and its associated intensity but also the direction of the change, e.g., whether the change was due to a target(s) that appeared or disappeared between the two dates. A faster computation of the same method is presented in [18]. Note that while knowing the direction of change provides more information, it does not explain in physical sense, i.e., as an SM, the type of change. This may be sufficient for some applications, specifically for change analysis of man-made targets, however, understanding the types of changes is crucial, e.g., for environmental applications, such as understanding the physical changes in the vegetation development or degradation of a crop or a forest.
Other studies introduced the concept of change detection matrix to denote changes not only between consecutive images but also between all the images in a stack. In [19], the authors used it for speckle filtering by identifying the changed and unchanged areas over a time series; however, the analysis was limited to single polarization data and focused exclusively on the image filtering strategy. In [20], the concept was extended for PolSAR data and included multitemporal change detection. However, besides being a supervised change detection methodology, the change matrix (CM) and the change maps only provide information about how dynamic or stable the resolution cells are over the time series without giving information about the types of changes and the specific dates of occurrences.
In [21], a PolSAR change detector inspired by the polarimetric matched filter [22] was introduced, which proposes to find changes between two images based on the optimization of the power ratio between them. In this case, focus is given to the intensity of the change for anomaly detection, although the SMs are not analyzed. Based on the same concept, a method to analyze and characterize multitemporal PolSAR changes was proposed recently [23], also based on the maximization of the power ratio between acquisitions, with the aim of providing not only the intensity of changes between images but also its type. A novel visualization of the PolSAR changes, previously introduced in [24] and [25] is also formalized.
The authors of [23] extend the method to not only analyze changes between consecutive images but also among all images in the stack and apply it to visually inspect changes in agricultural fields. Note, however, that the projection vectors of the optimization (i.e., eigenvectors of the optimization) that maximize the power ratio are not guaranteed to correspond to SMs but are the vectors that provide the maximum contrast between the covariance matrices of the two dates. This implies that the type of changes detected may not correspond to physical changes in the target and can only be interpreted with the aid of transformation models which significantly restricts the potential of the method.
In this article, a modified version of the change detector is proposed in order to improve the type of change interpretation so that the obtained changes represent PolSAR SMs directly. This allows one to gain insight into the scene evolution without depending on transformation models that consider matrix multiplications [26]. Note that part of this method was introduced in [26]- [28]. However, in this article, we also present the analysis for different crop types and test sites, confirming the transferability of the methodology proposed, and present a novel application of the same method in which we test the potential for image classification. We exploit the fact that analyzing the intensity and type of changes during the whole growing season and encoding them into a CM allows the characterization of a target, in this case a crop type, in terms of the increase and decrease of SMs over time. In other words, we propose a classifier that is not only based on states (at different dates) but also on the dynamics (or changes) of such states directly embedded in the change matrix. Unlike other crop-type classification algorithms, the method proposed here provides the advantage of having a single method for in-season phenology monitoring as shown in [23] and for crop-type classification as will be presented here.

II. TEST SITE AND DATASETS
The methods presented in this article are tested on two separated test sites where ground truth has been collected. The sites are geographically distant and include several crop types in order to include crop and region diversity. In the first site, the focus is on monitoring crop growth while in the second site, additional to crop growth the interest is to test the potential of the methods presented here for crop classification. The following is a short description of the two locations.

A. South of Spain
The first site corresponds to rice fields growing in parcels with an average size of 11 Hectares, located near Seville, in the south of Spain as shown in Fig. 1. The ground truth includes phenology data for validation collected weekly from May to September 2014 covering one complete agricultural season. In this case, the ground data are used for comparing and validating the results obtained from the proposed methodology, particularly for monitoring crop growth along the season as it will be shown in Section VI. With regard to the satellite data, as shown in Table I, five quad polarimetric images from the C-band RADARSAT-2 satellite are considered. The images are acquired with a temporal resolution of approximately 24 d and an incidence angle of 39 degrees.

B. Indian Head, Canada
This dataset includes several other crops types including barley, oats, wheat, canola, flax, lentils, and field peas, in parcels  with an average size of 70 Hectares. These crops are located in Indian Head, Canada, where ground truth was gathered as part of the ESA-funded AgriSAR 2009 campaign, which covers phenology among other biophysical variables. Images between June and September of the same year are selected, ensuring that the cultivation period is covered and the images affected by severe rainfalls are removed. Fig. 2 presents the location of the fields and Table II shows the seven quad polarimetric images from the C-band RADARSAT-2 satellite with the associated incidence angles that are considered for crop growth monitoring and crop type mapping (Sections VI-A2 and VII).

C. SAR Data Preprocessing
In both test sites, the original spatial resolution of the images is 5 × 5 m, however, the final resolution may be affected by the following SAR data pre-processing steps. The covariance matrices utilized as input for the analysis are obtained after applying radiometric calibration and speckle filter. A 9 × 9 boxcar filter is used considering that parcels are sufficiently big and several pixels fall within a parcel. Then, terrain correction and geo-coding is performed to ensure that all the images are geo-coded over the same grid even if they are acquired by different beams or pass directions.

III. METHODOLOGY
The scattering process measured by a quad-polarimetric radar imaging system is conventionally represented for each resolution The elements S ij of the matrix correspond to the complex scattering coefficients associated with the amplitude and phase of the backscattered signals, where H and V stand for linear horizontal and vertical and the double letter is for transmitterreceiver. The scattering matrix can also be represented in a vectorized form derived from a chosen set of 2 × 2 complex basis matrices, as k = 1 2 [T r(SΨ)], where T r is the matrix trace operator, and Ψ is a set of 2 × 2 complex basis matrices [1], [2].
If a monostatic backscattering system is used and assuming reciprocity (HV = V H), the polarimetric target vector in Pauli basis is given by where the superscript T denotes the vector transpose, and the √ 2 is required to preserve the total scattered power. The polarimetric target vector presented in (2) corresponds to the case of a single look complex image. The central limit theorem can be applied assuming that the number of scatterers inside a resolution cell is large, and thus, k p follows a multivariate complex circular Gaussian distribution with zero mean and probability density function as follows [1], [2]: where * represents the complex conjugate, . represents the matrix determinant, and Σ is the covariance matrix of the Pol-SAR target vector in Pauli basis (also known as the coherency matrix), obtained as Σ = k p k * T p , where . is the statistical expectation operator.
Equation (3) corresponds to the distribution k p ∼ N (0, Σ), where Σ is a Hermitian positive semidefinite matrix and contains the necessary information to characterize a target [1], [2].
Since the amplitude and phase of a resolution cell are the coherent and linear combination of backscattered signals from individual scatterers within it, the measurements are affected by a random variation denoted as speckle [29]. To reduce the randomness of the acquired signals, L number of iid resolution cells can be averaged (or speckle filtered), and the matrix of (3) becomes the L-looked coherency matrix where i=1,2,..., L is the number of averaged samples or realizations.

A. Covariance Matrix Eigenvector/Eigenvalue Decomposition
Several theorems have been proposed to decompose the target covariance matrix of (4) into simpler objects to aid its physical interpretation. These can be coherent decompositions to characterize the so-called coherent or pure targets, or incoherent decompositions for partial or distributed targets [2]. In this article, a brief introduction to the eigenvector/eigenvalues decomposition in the context of PolSAR data analysis [2], [7] is given since as it will be seen in Section IV, it is linked to the optimization required in the change detector proposed. The coherency matrix T of (4) can be decomposed using the eigendecomposition theorem as where D is a 3 × 3 diagonal matrix that contains the nonnegative and real eigenvalues diag(λ 1 , λ 2 , λ 3 ), and U = [u 1 u 2 u 3 ] is a 3 × 3 unitary matrix in which the column vectors u 1 , u 2 , and u 3 are the three orthogonal eigenvectors.
In order to provide a physical interpretation of this decomposition, the eigenvectors u 1 , u 2 , and u 3 or more generally, v i with i = 1, 2, 3 of T can be rewritten as The average or dominant SM of a resolution cell can then be obtained by combination of the three orthogonal eigenvectors as follows [2]: In the three later definitions, P i corresponds to the pseudo probabilities Using (7), it is possible to visualize the average or dominant SM of each resolution cell of a PolSAR image as an RGB composite as shown in Fig. 1 using the following colour representation:

IV. POLSAR CHANGE DETECTION BASED ON DIFFERENCE OF COVARIANCE MATRICES
This section introduces the proposed PolSAR change detector which allows us to extract polarimetric changes between two acquisitions, identifying the occurrence of changes and providing interpretation about the type of change that a target suffered.

A. Pairwise Change Detection
For the case of bidate change detection, we denote the pair of coregistered acquisitions to be compared as date1 and date2. The coherency matrices represented by (4) and denoted as T 1 and T 2 are used to characterize a resolution cell at date1 and date2, respectively.
It is possible to write T 2 as the sum of T 1 plus a separate matrix including the added and subtracted components as where the matrix T C calculated as T C = T 2 − T 1 is Hermitian symmetric, since it is produced by the linear combination of two Hermitian matrices, but may not be positive semidefinite. This means that the diagonal elements of T C are real and the upper triangular part is the complex conjugate of the lower triangular part [27]. However, there is a difference with ordinary coherency matrix, since it is not bound to be positive semidefinite. This means that its trace and its quadratic form P = ω * T T C ω (using a generic projection vector ω) can be negative. This is an expected feature since we need a matrix that is able to communicate if the change in the partial target has brought an increase or reduction of a) a specific SM and b) the overall power of the final partial target (trace).
By representing T C with its quadratic form P = ω * T T C ω, we can investigate the amount of change that a scattering vector represented by the projection vector ω suffers. In this sense, by optimizing over all the possible projection vectors ω, it is possible to find the one that experiences the largest or smallest change. Note that the optima ω also correspond to the SM that was added/subtracted to the partial target in date1.

B. Optimization
To find the maximum and minimum projection vectors, we can apply the well-known Lagrangian optimization for the quadratic form P = ω * T T C ω, i.e., By constraining ω to be unitary we can obtain the Lagrangian as Differentiating with respect to ω * T and setting the derivative equal to zero, we obtain the equation Note that (13) corresponds to the eigendecomposition theorem described in (5), and therefore, the optimization can be completed by a diagonalization of the matrix T C . Since T C is Hermitian, the eigenvalues will be real but as T C is not positive semidefinite, the eigenvalues are not bound to be positive. This is because a change can increase or decrease the resulting power of a SM (e.g., we could have that the surface SM increases or decreases).

C. Change Visualization
After the optimization procedure described in the previous section, the eigenvalues and eigenvectors of T C are obtained. In this case, the eigenvectors represent SMs that change between a pair of acquisitions and the eigenvalues represent the intensity of such changes. Note that since the eigenvalues are not bound to be positive, the negative eigenvalues represent a SM that has been removed from the scene.
The visualization and interpretation can be performed in the same way as described in Section III-A. The dominant or average SM for both cases, i.e., when is added or removed from a scene, can be represented with (9). However, (8) requires special consideration given that the eigenvalues can be negative. It is adjusted by separating the positive and negative eigenvalues corresponding to added or removed SMs, respectively, which allows for their independent visualization as follows: Equation (14) indicates that if an eigenvalue is negative, the pseudoprobability that the associated eigenvector or SM was added to the scene is zero and if the eigenvalue is positive, the pseudoprobability that the associated eigenvector or SM was removed to the scene is zero. We can then compute the average parameters λ, α, β separately for added or removed SMs which enable us to employ the RGB representation of (9) as follows: Using (15) and (16) into (9), an RGB representation for the added and removed SMs between a couple of acquisitions can be obtained. Fig. 3 shows the change detection results applied over four different pairs of images along the season in the test site 1. To describe the changes between images, we will focus on the rice fields which are represented by the pink color in Fig. 1. Although other interesting behavior is observed in the smaller fields nearby, no description is given since ground truth for these fields is not available in this dataset. A single incidence angle is considered to provide an interpretation of the observed PolSAR changes between the following set of pairs of images: Pair 1: Fig. 3(a) shows the change in SMs from the beginning to the end of June when the phenological information confirms that on 2014-06-05 most of the parcels have been already flooded while on 2014-06-29, the majority of the parcels reached the tillering stage. In the top plot of Fig. 3(a), we can see a strong increase in double bounce (red) mostly due to the SAR signal bouncing in the water surface, then to the emerging plants in tillering and back to the satellite. Few parcels also depict blue color since not all parcels share the same starting date or might be shortly delayed in the growth process. Note, at the bottom plot of Fig. 3(a) that mostly dark areas are present meaning that the SMs that are removed from the scene are not strong. This can be associated with the fact that the parcels transitioned from flooded fields to the fields with short vegetation and, therefore, the power of the acquisition is stronger at the end of June than at the beginning, causing the majority of the SMs to be added rather than removed.

V. RESULTS AND INTERPRETATION OF POLSAR CHANGE DETECTION
Pair 2: A similar behavior can be observed for change in SMs from the June 29 to July 23 when the fields reach the vegetative stage and plants continue to get more vigorous. However, the intensity of the SMs that are added to the scene is not as strong as the change between the previous pair since the overall change in the backscatter power is not as high. Note that the color depicted in Fig. 3(b) resembles more of a pink color rather than only red, indicating a combination of added double bounce and surface scattering. This might be explained since stronger and more vegetative structures are added to the scattering scene due to the increase in the plants height and the number of tillers. Note also that this does not necessarily mean that surface scattering is predominant but that the surface scattering in July was stronger than in June. This may be a consequence of the ability of the detector to identify an increase of surface scattering even when a surface scattering was already present in the scene. This can happen when, for instance, the roughness or the moisture increased between images, and it is flagging that one of these parameters has changed. In this case, the surface scattering increases because between these stages, the soil goes from a smoother water surface to a rougher surface. The detector is telling us that the smooth water surface is no longer present. With regard to the SMs that are removed between June and July, a combination of dark areas and weak red and pink colors can be seen in the bottom plot of Fig. 3. The parcels in which the dark areas dominate indicate that no SMs are being removed from these parcels while the areas in red indicate the removal of double bounce. Note that the latter removal occurs in those parcels that are ahead in the development process and the combination of the dark areas with red/pink areas indicates the differences in development rate between these two types of parcels. The meaning of removing the pink color is better explained observing the pair number 3 (next pair), since it is caused by the same effect.
Pair 3: The pair from the July 23 to August 16 shows mostly dark areas indicating no added SMs to the scene [ Fig. 3(c)]. This can be explained since the plants are reaching the advanced vegetative state on the August 16th image, therefore, starting to lose the vigor and vertical structure and producing a less strong backscatter response than on the July 23rd image. On the contrary, the bottom plot shows that the combined double bounce and surface scattering that had been added from June to July are now being removed from the scene given the more mature condition of the canopy structure and the lower stem and leave water content of the plants, once they reach late vegetative stages. This is the only pair of images in which the removed SMs are stronger than the added ones.
Pair 4: Fig. 3(d) shows an increase in green color in the top plot, corresponding to the addition of volume scattering from August 16 to September 19. By the latter date, plants have reached maturation, lost water content, and the vigoros vertical structure allowing the SAR signal to interact with plants, soil, and simultaneously with both in a combination of the three SMs. In the bottom plot of Fig. 3(d), it is possible to see that it is mostly characterized by dark areas with again weak pink color being removed from some parcels, in this case, due to these parcels being delayed in the growth process and showing the same effect than the removed SMs in the previous pair.

VI. CM AND TARGET DYNAMICS
Since the information about the crop evolution is obtained from the analysis of pairs of acquisitions, it is possible to evaluate not only changes from consecutive dates but also to see how any given acquisition differs with respect to all the other acquisitions in the stack of coregistered images. This process identifies the evolution of the SMs throughout the season.
Inspired by [23], we organized the PolSAR changes not only between consecutive images but also between all the images in the stack as a N × N ×3 square matrix, with N being the number of images in the stack, and three given the RGB representation of (9). The OFF-diagonal elements correspond to the evaluation of polarimetric changes between the available dates. The upper triangular part is designed to represent the SMs that are added to the scene between a given pair of images and the lower triangular part to represent the ones that have been removed. Note that the diagonal elements of this matrix can also be employed if desired by obtaining the dominant or average SM in the scene for a given date as introduced in Section III-A.
It is worth noticing that the color representation of added or removed SMs in the square matrix is in the same way as described in (9). Therefore, the CM is an RGB composite that depicts with the color the type of change and with the contrast the intensity of the change between a given pair of dates.

1) Target Dynamic Analysis:
To illustrate the evolution of a target SM over time and relate them to physical changes, Fig. 4 shows the SMs throughout the season of a rice parcel in the test site 1, where the ground truth is known. The colors and intensities of the CM can be interpreted as a conventional Pauli RGB composite, i.e., red for double bounce, green for oriented dihedral or volume, and blue for surface scattering.
Due to the rice crop cultivation practices in which the season starts with seeding over flooded soils with low backscatter returned to the satellite, the multitemporal change detection process mostly identifies increase in SMs when compared to this initial stage. This is reflected by the upper triangular elements of the CM of Fig. 4, having marginally stronger intensity than the lower triangular part. In Fig. 4, the squares and circles with the numbers 1 to 5 correspond to the main growth stages of rice. The SMs added between two growth stages can be seen as the intersection between the squares of both stages in the upper triangular part. The intersection of these stages in the lower triangular part represents the SMs that are removed from the scene. It is possible to see, for instance, that double bounce is added to the scene between stages 1 and 2 represented by the red color, which is assumed to be caused by the interaction of the SAR signal with the plant stems emerging from the water at tillering stage. Although it is almost imperceptible, dark blue appears in the left most column of the CM indicating that surface scattering is removed from all stages when compared to square 1 (flooded soil).
When the crop reaches the advanced vegetative stage in July, more vegetative structures are added to the scattering scene due to the increase in the plants' height and the number of tillers. As a consequence, there is an increase in the power of the returned backscatter and addition of both surface scattering and double bounce is observed. This combination is reflected in the CM with purple color in the intersection between squares 1 to 3 and 2 to 3. On the contrary, the biggest decrease in the intensity of change occurs in the transition between vegetative to reproductive stages (end of July to end of August). This is observed due to the plants starting to lose the vigoros vertical structure when entering into the reproductive stage. The CM shows, at this point, decrease in surface and double bounce in the intersection between squares 3 and 4 in the lower triangular. The intersections between 3 and 4 and between 4 and 5 in the upper triangular part remain dark indicating that the intensity of the SMs added between these dates is not strong. Note that rather than showing dark squares, these are shown as a smooth transition between green to blue and again back from blue to green (from 4 to 5) due to the linear interpolation used to smooth up the transition between dates.
At the end of the season (rightmost column), an increase in green color corresponding to volume scattering being added to the scene when the crop reaches maturation can be seen. This is communicating that comparing stage 5 with any stage from 1 to 4, volume scattering is added. This is as a consequence of the plants' random orientation and also due to the plants getting drier and, therefore, more transparent to the SAR signal causing multiple scattering to be present in the scene, and considering that the soil is not flooded as it was at the beginning of the season. Traditional polarimetric analysis has reported similar results with high values of entropy at this stage [4], [30].
2) Other Crop and Land Types: To generalize the usefulness of the CM for target dynamics analysis, we obtained the CM for the crops included in the test site 2. The obtained results are shown in Fig. 5, and a detailed interpretation of the CM for field pea crops together with the corresponding field photos is presented in the supplementary material accompanying this article.
It is possible to see by visual inspection how different the resulting matrices are for each of these three crop types and to that of the rice crop presented in Fig. 4. This shows how the matrix encodes differently the temporal PolSAR evolution for each crop type characterizing it in a specific way. The analysis correlating the crop morphology at the different growth stages can also be used with these change matrices to understand the interaction of the SAR signal with each crop throughout the season.
It is important to notice that because the matrix contains separately each SM (in each channel), the analysis of the parcel evolution can be done for each mechanism separately. For instance, we can compare if the double bounce in a parcel due to the plants' emergence above the water level is delayed in time relative to other parcels or a reference parcel.
On the other hand, similar effects were observed for other land types, such as rivers, forests, and cities where each of these land types created a CM with specific behavior. These differences described by the colors that represent the evolution of SMs, the intensity of the CM, and the times when specific events occur can be used to characterize a crop type or more generally a land class and, consequently, this information can be exploited as input for image classification.

VII. CROP-TYPE CLASSIFICATION
It was shown in Section VI how the CM can be used to visualize, understand, and characterize the interaction and response This dataset is split according to the percentages described in Section VII-A1 of the radar signal with a target that evolves over time. This section presents the application of the CM for multitemporal image classification. It considers the proven fact that multiple acquisitions provide better performance results than single image classification while also considering the fully polarimetric information measured by the satellite [31], [32]. In order to investigate the CM potential for image classification, change matrices are used as inputs for a fully connected neural network (NN) classifier. Its performance is then compared with a classification based on time series of PolSAR features, where the features are derived from covariance matrix decompositions [2].
A. Dataset, Metrics, and Classifier 1) Dataset Splits: In all cases, the parcels, where ground truth is available, are divided into three groups: training (65%), validation (15%), and test parcels (20%) preserving in all cases the same proportion of samples per class as in the original dataset. Note that we split the data ensuring that the resolution cells belonging to a parcel are only present in one of the sets (i.e., training parcels independent from validation and test parcels) considering that neighboring pixels would have correlation due to spatial filters and other preprocessing steps, such as multilooking, speckle filter, and coregistration. Without doing this, the test set would contain samples very similar to those in the training set and the model would give overoptimistic performance after predicting on test set.
After splitting the available ground truth and removing the crop types where the data is not enough for training and testing a model, the dataset is reduced to ten crop types and results in a highly imbalanced distribution of pixels, as shown in Table III. Consequently, both the classifier and the metrics used are adjusted to account for this effect as introduced in next section.
2) NN Classifiers: A fully connected NN is used for classification in all the cases presented in the following sections. Hyperparameter search was performed to select network depth, width and need for dropout layers, and to select the optimizer and learning rate. This resulted in an NN with three hidden layers of 256 neurons each, using batch normalization before every hidden layer and before the output layer, and Nadam optimizer with learning a rate of 0.0005. To deal with the class imbalance as previously stated, the model is trained using weighted cross entropy loss function, in which the weights correspond to ratios obtained from the class imbalance. Note that randomly subsampling the majority classes before training was also tested, delivering similar or slightly lower performances. On the other hand, Random Forest based classifiers were also tested with lower performance.
3) Metrics: The accuracy metric is traditionally reported to evaluate classifiers performance. However, since it relies on the frequency with which the model predictions on the test set match the ground truth, in imbalanced classification problems the results are dominated by the majority classes and with small contribution by the minority classes [33], [34]. To account for this, we report the balanced accuracy metric which is the accuracy adjusted by their corresponding class weight (ratio of samples per class over the total number of samples). We also report the per-class F1 score, recommended for imbalanced classification [33], [34], and compute the macro average F1 measure to average across classes.
To be able to compare the obtained results with other works, we additionally report the traditional overall accuracy which we call imbalanced accuracy, since it does not account for the class imbalance present in the data.

B. CM-Based Classification
The CM of size N × N × 3 derived in Section VI, where N is the number of dates, is flattened as a 1-D array and received by the input layer of the NN. For the present analysis, seven images corresponding to the dates shown in Fig. 5 are used, thus resulting in a 147-input feature vector. Fig. 6 shows the confusion matrix obtained after predicting on the test set. It can be noticed that the CM-based classifier performs well predicting crops such as canola, field peas, lentil, and mixed pasture and that it underperforms predicting crops such as barley, durum wheat, oats, and spring wheat, mostly due to confusion in the prediction between these four classes. This is in part explained considering that these crop types have similar biophysical characteristics being all of them part of the family of cereal crops, therefore, producing similar PolSAR responses. Table IV shows the corresponding metrics obtained on the test set after training the models. It can be seen in the change matrix column that both the overall balanced accuracy and the F1 macro metrics are relatively low mainly due to the confusions made by the classifier within the cereal crops predictions. Note that the overall imbalanced accuracy gives the impression of a better classifier performance, although in this particular case, it can be attributed to the class imbalance in the dataset. If instead of classifying these cereal crops separately but in a single class, the classification performance increases significantly as presented in Fig. 7. In this case, most classes are predicted correctly, being flax the crop less accurately predicted with some confusions with the cereal crops and pastures. Also, notice that the performance of the classifier predicting canola in the test set is very high since this crop structure and its evolution over time is considerably different with respect to the other crops under analysis.

C. PolSAR-Features Based Classification
Time series of backscatter and PolSAR features derived from the H/A/alpha decomposition of the coherency matrix [7], [35], were obtained and used as input for an additional Neural network to compare its performance with that of the CM-based classification. A total of 16 features time series were considered, including the VV, VH, HH backscatter, the ratios VH/VV, HH/VH, and HH/VV, dominant and average eigenvalues, dominant and average alpha angles, dominant beta angle, entropy, anistropy and the magnitude of the SMs computed as √ λ * cos(α), √ λ * sin(α) * cos(β), √ λ * sin(α) * sin(β). Fig. 8 shows the confusion matrix for this case, with most classes having a similar performance than the CM-based classification, except for Canary seed and field peas where the CM-based classifier seems to have slightly better performance, and showing that using PolSAR time series features the classifier also has confusions identifying the cereal crops.
In terms of the overall performance, Table IV shows that the classifiers using both input data types (i.e., change matrices or PolSAR time series features) seem to have very similar resulting test metrics, when predicting ten crop types. The balanced accuracy as well as the F1 macro, which are particularly important for imbalanced classification, are slightly better for the CM-based classifier. Note, however, that because the differences are not large enough (around 1%), a definite conclusion about which input data type performs best is not possible.
After grouping the cereal crops in a single class and training an NN classifier to predict the reduced set of classes, the results obtained are shown in Fig. 9. Again significant improvements are achieved compared to a classifier predicting the cereal crops separately and improvements in all crop types apart from mixed pasture are also achieved with respect to the CM-based classification. With regard to overall performance for the case of the reduced set of classes, Table V shows that the classifiers using both input data types have similar resulting test metrics, however, note that both the balanced accuracy and the F1 macro are approximately 3% better for the PolSAR time series based classifier than for the CM-based equivalent. This is opposite to the results of Table IV where the CM-based outperforms the PolSAR-features based classification. Since neither of the  two tables show large differences between metrics, a definite conclusion about which input data type performs better is not possible but rather we see that the performances are comparable.

D. Prediction Maps
The trained models can now be used for predicting on the whole test site as shown in Figs. 10 and 11 for the 10 classes case, i.e., the cereals crops as independent classes, and for the 6 classes by combining the cereal crop types, respectively.
In Fig 10, we notice similar results to those obtained in the confusion matrices, in which canola, lentil, field peas, and mixed pasture have strong agreement with the ground truth, whereas the cereal crops are confused between them. Fig. 11 also shows the confusion of the model to predict flax as a cereal crop.

VIII. DISCUSSION
This article addresses a method to process and analyze multitemporal and quad polarimetric SAR data. This aims at developing a methodology that helps in reducing the complexity when analyzing time series of polarimetric features. As shown in Section VI, inspired by [23], we use quad-pol change detection to encode the evolution of SMs over time in a CM. The proposed CM can be interpreted intuitively by looking at its colors and intensities. In this article, in order to have a straightforward interpretation of the SMs, a change detector based on the difference of the covariance matrices is used as opposed to the power ratio used in [23]. Note that a detailed analysis between  these two approaches will be addressed in subsequent work. Compared to previous approaches that propose a multitemporal change detection [15], [16] and change matrices [19], [20], the methodology proposed in this article also introduces an improvement, since it not only exploits intensity of changes but also the types of changes.
In Section VII, we addressed the performance for image classification using the CM as input. This was compared to more traditional classifiers that use multitemporal quad-pol features as inputs. Experimental results indicated that comparable levels can be obtained for image classification purposes when using the CM as input for a machine learning classifier. However, the CM allows one to have an interpretable idea of the input data being used, while simultaneously providing a single method for crop monitoring over time and crop type classification.
A current limitation of the proposed method is that the crop stage monitoring was addressed by visually linking phenological stages rather than a more granular index of phenology or other crop biophysical variable. This aspect is currently being addressed and will be approached in future work, where a comparison with traditional methods will be presented. In addition, the influence of the incidence angle will be further investigated. This is because combining different acquisition geometries currently results in undesired artifacts due to combined detection of changes in both the crop growth and the incidence angle. Using the CM for crop growth, a simple solution would be to restrict the detection of changes to the same orbit as in Fig. 4, resulting in increasing the time between images. However, for current constellations, this might be acceptable since it can reach times of six days, e.g., for the case of Sentinel-1. A different approach could be to perform incidence angle normalization before the change detection step as, for instance, in [36]. Note that the influence of the incidence angle when using the CM for crop type classification is expected to be less significant. This is because both the training and test pixels (or pixels in the region of interest) contain the effect of the incidence angle. Therefore, a classifier could predict crop types since it receives the same information at training and testing times. The acquisition geometry, however, requires more detailed consideration and will be addressed in future work. Also for forthcoming investigation, other methods that create multitemporal matrices can be integrated to the present methodology. These methods, recently introduced in the literature, use interferometric coherence between images [37], [38] instead of the polarimetric information. Therefore, combining the PolSAR and InSAR data may result in complementary information being used as input for the classifier in order to boost performance. Similarly, a next step developing the change matrices methodology is to adapt it to dual pol systems such as the Sentinel-1 satellite. This is done by replacing the 3 × 3 pixel covariance matrices derived for monostatic Quad-pol systems of equation (10), for the 2 × 2 matrices formed using the two polarisations provided by Sentinel-1. Note, however, that the visualization and interpretation need to be adapted accordingly.
On the other hand, we have noticed that differences in color and intensity and deviations from a typical CM can be identified visually. This can also be used as a tool for in-season crop growth anomaly detection and will be the focus of future research. Regarding the classification, the use of quad-pol features (focused on the system states) and change detection matrix (focused on the system change of states) provide similar performance, but they are quite complementary. As a future work, we want to study an optimal way to merge the two information sources together in a classifier.

IX. CONCLUSION
This article presents a change detection based method to visualize and analyze multitemporal quad polarimetric images. The change detector is derived as the pixelwise difference of covariance matrices between all the coregistered images in a stack, which organized in a matrix form, allows the interpretation of the evolution of SMs over time. We showed how growth stages of rice fields crops can be related to the SMs observed and how it differs from other crop types. Based on this property of encoding the SMs for each crop type, we tested the performance of the proposed method for image classification to identify several crop types. Results show that this method yields similar performance to traditional classifiers, which use time series of polarimetric features as input with differences in the overall balanced accuracies and F1-macro metrics of around 1% or 2% in favor of one or another based on the number of classes used. The method presented here achieves similar classification performances while providing additional advantages in terms of interpretability and insight into the physical changes of a target over time. It also provides a different view of the change of crop states that could be incorporated to provide a broader view of the whole dynamic system.