Unsupervised Classification of Crop Growth Stages with Scattering Parameters from Dual-Pol Sentinel-1 SAR Data

Abstract

Global crop mapping and monitoring requires high-resolution spatio-temporal information. In this regard, dual polarimetric Synthetic Aperture Radar (SAR) sensors provide high temporal and high spatial resolutions with large swath width. Generally, crop phenological development studies utilized SAR backscatter intensity-based descriptors. However, these descriptors are derived either from the covariance matrix elements or from the eigendecomposition. Therefore, this approach fails to utilize the complete polarization information of the scattered wave. In this study, we propose a target characterization parameter, ${\theta }_{\mathrm{xP}}$ that utilizes the 2D Barakat degree of polarization and the elements of the covariance matrix. We also propose an unsupervised clustering scheme using ${\theta }_{\mathrm{xP}}$ and the scattering entropy, $H_{\mathrm{xP}}$. We utilize time-series Sentinel-1 data of canola and wheat fields over a Canadian test site to show the sensitivity of ${\theta }_{\mathrm{xP}}$ to the development of crop morphology at different phenological stages. During the initial growth stages, ${\theta }_{\mathrm{xP}}$ values are low due to the low vegetation density. In contrast, at advanced phenological stages, we observe decreased values of ${\theta }_{\mathrm{xP}}$ due to the appearance of complex canopy structure. Similarly, the effectiveness of the unsupervised $H_{\mathrm{xP}}/{\theta }_{\mathrm{xP}}$ clustering plane is also evident from the temporal clustering plots. This innovative clustering framework is beneficial for the operational use of Sentinel-1 SAR data for agricultural applications.

1. Introduction

Identification and monitoring of crop phenological stages are essential factors in agriculture for estimating crop production. In this context, Synthetic Aperture Radar (SAR) data have been extensively used, especially for crop classification [1,2], yield estimation, and biophysical parameter retrieval [3,4,5], formulation of vegetation indices [6,7] and descriptors [8]. Moreover, SAR sensors have high spatial resolution and all-weather imaging capability. The scattered SAR signal is significantly affected by the geometry of crop canopy, underlying soil roughness, and dielectric properties of both crop and soil. In addition to the acquisition wavelength, polarization and angle of incidence have a significant influence on the received information [9,10]. Over agricultural fields, crop biophysical parameters impact SAR backscatter, particularly at advanced vegetative stages. Soil roughness and moisture remain the governing factors for SAR backscatter response during initial growth periods [11].

In particular, Sentinel-1 dual-pol images provide opportunities for crop classification and biophysical retrieval-related applications. Sentinel-1 also provides large swaths and lower data volumes due to reduced polarimetric dimensionality [12]. In the literature, dual-pol intensity parameters are widely used for crop type identification [1,13,14,15,16,17,18,19]. The existing studies also extend to retrieval of crop biophysical parameters [4,20,21] and crop phenology monitoring [22,23].

For dual polarimetric HH-HV or VV-VH data, Cloude et al. [24] proposed an eigen-based decomposition technique. In this technique, the 2 × 2 covariance matrix is decomposed into two orthogonal eigenvectors that are further utilized to derive a target characterization parameter. This average scattering angle $\overline{\alpha }$ is obtained from the two orthogonal polarization states weighted by their respective pseudo probabilities obtained from the eigenvalues. The target scattering entropy $H_{\mathrm{xP}}$ is obtained from the pseudo probabilities. This target characterization parameter and $H_{\mathrm{xP}}$ are further utilized to propose an unsupervised clustering scheme. Similarly, Ainsworth et al. [25] proposed another characterization parameter, $\theta$, utilizing an eigen-based technique. The parameter is related to the cross, ${\sigma }_{\mathrm{XY}}^{\circ }$ and co-pol, ${\sigma }_{\mathrm{XX}}^{\circ }$ ratios. In addition, $\theta$ and $H_{\mathrm{xP}}$ are used to propose a clustering technique for scattering target identification. In their study, the clustering plane consists of eight zones that represent different scattering characteristics.

Sugimoto et al. [26] provided a comparison between polarimetric parameters obtained from dual and full polarimetric SAR data. They reported a high 2D correlation between the parameters derived from the full-pol $H/\alpha$ decomposition and the dual-pol (HH+VV) $H/\alpha$ decomposition. For the full polarimetric data, the analysis was performed using model-based three [27], and four [28] component decompositions. Several other descriptors were introduced for diverse applications: clustering [29], maritime applications [30,31], etc.

To monitor the growth of vegetation, different vegetation indices, such as the Radar Vegetation Index (RVI) [32], Dual-Pol SAR Vegetation Index (DPSVI) [33], and Dual-pol Radar Vegetation Index (DpRVI) [34], can be derived. In general, these descriptors capture the randomness from the complex crop canopy to describe their phenological stages. Apart from these, machine learning models have shown promising results for crop biophysical parameter estimation [22,35,36]. It was observed from these studies that the cross-pol ratio is the most important parameter for vegetation studies. Nasrallah et al. [37] utilized the Gaussian function on the time series SAR data to find the date of significant phenology stages for wheat. Wali et al. [38] analyzed the sensitivity of the temporal backscatter for rice crops using the line regression technique. Similarly, the interferometric phase for identifying wheat phenological stages is observed in a study by Schlund et al. [39].

Usually, polarimetric parameters have been suitably attributed to the physical properties of the crop canopy [40,41,42], and have therefore helped monitor crop phenology. Recently, Dey et al. [43,44] proposed a model-free target characterization parameter using full, compact and dual co-pol data [22]. This parameter utilizes the Barakat degree of polarization [45] and elements of the coherency matrix to characterize diverse target types. In their study, it was shown that this parameter enhanced the target characterization capability as compared to the $\overline{\alpha }$ parameter.

Motivated by this work, we derive a target characterization parameter for dualpol SAR data in this study. This parameter jointly utilizes the 2D Barakat degree of polarization and the elements of the covariance matrix. Therefore, this parameter, ${\theta }_{\mathrm{xP}}$, captures the additional scattering information which the other existing descriptors might miss. Alongside this, we also propose an unsupervised clustering technique using ${\theta }_{\mathrm{xP}}$ and entropy, $H_{\mathrm{xP}}$ parameters. We utilize ${\theta }_{\mathrm{xP}}$ and the clustering scheme to analyze the temporal variation of canola and wheat over a Canadian test site derived from Sentinel-1 C-band SAR data.

2. Methodology

2.1. Study Area and Dataset

The test site in Canada is located in southern Manitoba (49°40′32″N, 97°59′57″S). This is one of the Canadian Joint Experiment for Crop Assessment and Monitoring (JECAM) test sites. Nominal field sizes range from 20–30 hectares to 50–60 hectares, with wheat, soybean, canola, and corn as major annual crops. The land under permanent grassland is $<5$. The in-situ measurements were conducted during Soil Moisture Active Passive Validation Experiment 2016 Manitoba (SMAPVEX16-MB) campaign. Figure 1 shows the location of the test site and the distribution of canola and wheat fields.

The in-situ measurement campaign was carried out from June to July 2016 over 50 different fields. The measurement period mainly consisted of two temporal windows: 8 June to 22 June and 8 July to 22 July 2016. Within this temporal extent, most crops advanced from early to advanced phenological stages achieving peak biomass. The nominal size of each field is approximately $800 \mathrm{m}\times 800 \mathrm{m}$. Vegetation sampling, which includes biomass, Plant Area Index (PAI), and plant height, was performed in three points for each field. Crop biomass was collected via destructive sampling. A 0.5 m × 0.5 m square was placed over the canopy for canola and wheat. All above-ground biomass was collected by cutting all vegetation at the soil level within the square. The PAI was measured with hemispherical digital photos. In this technique, a camera with a fisheye lens captures photos of the crop canopy with the camera positioned at least 50 cm above or below the canopy. Details of the in-situ biophysical parameters (viz. PAI, Dry biomass and VWC) of Canola and Wheat are provided in Appendix A.3. One can find further details on the sampling schema in the SMPAVEX16-MB field report [46].

2.2. SAR Data Pre-Processing

We acquired the Terrain Observation with Progressive Scans SAR (TOPSAR) mode Level-1 Sentinel-1 Interferometric Wide (IW) swath Single Look Complex (SLC) data over the test site. To this SLC configuration, the swath length is ≈250 km and the spatial resolution is $5 \mathrm{m}\times 20 \mathrm{m}$. The whole IW swath is divided into three sub-swaths (IW1, IW2, and IW3) with nine bursts in the azimuth direction. Hence, prior to application, these images are pre-processed with standard correction steps [47].

The Sentinel-1 images are ingested in the SNAP 8.0 platform. Following this, the sub-swaths and bursts are selected based on the test site location, and the state vectors are updated with the Sentinel-1 precise orbit file. The images are then radiometrically calibrated. As we are interested in generating the complex covariance ${\mathbf{C}}_2$ matrix information, and hence, we need to extract both amplitude and phase values. Due to this reason, the calibration output is set to complex output values.

All Sentinel-1 images are then back-geocoded using the Sentinel-1 Back Geocoding operator. This operation essentially co-registers the temporal images with sub-pixel accuracy. We utilize the SRTM 1Sec Grid as the Digital Elevation Model (DEM). Subsequently, the Sentinel-1 TOPS deburst and merge operations are performed to produce a single SLC image. The stack of images is then clipped within the in-situ measurement location to reduce the data volume and increase the computation speed.

These image subsets are multi-looked with a factor of 4 in range and 1 in azimuth directions to generate a square pixel. The final pixel resolution is ≈ 15 $\mathrm{m}$. Finally, these multi-looked images are used to generate covariance elements images. A $5\times 5$ boxcar filter is applied to further reduce the speckle information within the images. Next, the baseline information is deleted from the metadata, and the covariance images are exported in the PolSARpro format. Here, please note that the covariance elements have ensembled information that helps in applying the second-order statistics.

Then, we compute the target characterization parameter, ${\theta }_{\mathrm{xP}}$ using the PolSARtools software in the QGIS platform [48]. Subsequently, we geocode the images with UTM projected coordinate system. We further analyzed the in-situ measurement locations and extracted target characterization parameters from the geocoded products. From several Sentinel-1 images acquired during the campaign, four dual-polarization (VV and VH) C-band Sentinel-1 Single Look Complex (SLC) data were selected for use in this study (

Table 1. Specification for Sentinel-1 data acquired for the Canadian test site.

Acquisition Date	Beam Mode	Incidence Angle Range (Deg.)	Orbit	$\mathbf{az}(\mathit{m})\times \mathbf{rg}(\mathit{m})$
13 June 2016	IW	30.22–32.47	Ascending	$15\times 15$
07 July 2016	IW	30.22–32.47	Ascending	$15\times 15$
19 July 2016	IW	30.22–32.47	Ascending	$15\times 15$
24 August 2016	IW	30.22–32.47	Ascending	$15\times 15$

). The selection of Sentinel-1 datasets was based on acquisition dates that were near coincident with in-situ measurement periods.

2.3. Target Characterization Parameter

The scattering vector for a dual-polarized SAR data are represented as, $k_d={\left[S_{XX}{S}_{XY}\right]}^T$, where X and Y correspond to either horizontal (H) or vertical (V) polarization states. Using these scattering vectors, we define the $2\times 2$ covariance matrix as

$${\mathbf{C}}_2=\langle k_d·k_d^{*T}\rangle =\left[\begin{array}{cc}\langle |S_{XX}{|}^2\rangle & \langle S_{XX}{S}_{XY}^{*}\rangle \\ \langle S_{XY}{S}_{XX}^{*}\rangle & \langle |S_{XY}{|}^2\rangle \end{array}\right]$$

(1)

where $\langle ·\rangle$ denotes ensemble average, and T denotes vector transpose. We define the elements of the $2\times 2$ matrix as, ${\mathbf{C}}_{11}=\langle |S_{XX}{|}^2\rangle ,{\mathbf{C}}_{22}=\langle |S_{XY}{|}^2\rangle$ and ${\mathbf{C}}_{12}=\langle S_{XX}{S}_{XY}^{*}\rangle$.

Similar to the conventional degree of polarization, the 2D Barakat degree of polarization [49] given in Equation (2) ($0\le m_{\mathrm{xP}}\le 1$) also characterizes the state of polarization (or purity) of an EM wave. For a completely polarized EM wave, $m_{\mathrm{xP}}=1$ and for a completely unpolarized EM wave, $m_{\mathrm{xP}}=0$. In between these two extreme cases, the EM wave is said to be partially polarized, $0<m_{\mathrm{xP}}<1$:

$$m_{\mathrm{xP}}=\sqrt{1-\frac{4|{\mathbf{C}}_2|}{{\left(\mathrm{Tr}\left({\mathbf{C}}_2\right)\right)}^2}},$$

(2)

where $\left|·\right|$ is the determinant of a matrix and $\mathrm{Tr}$ is the trace of a matrix. Let us consider two auxiliary quantities defined as

$$tan{\eta }_1=\frac{C_{11}}{m_{\mathrm{xP}}\mathrm{Span}}\mathrm{and}tan{\eta }_2=\frac{C_{22}}{m_{\mathrm{xP}}\mathrm{Span}},$$

(3)

where ${\eta }_1$ and ${\eta }_2$ are two free variables, and $\mathrm{Span}=C_{11}+C_{22}$. The quantity ${\displaystyle \frac{C_{11}}{m_{\mathrm{xP}}\mathrm{Span}}}$ indicates the co-polarized scattering power with respect to the total polarized power, while ${\displaystyle \frac{C_{22}}{m_{\mathrm{xP}}\mathrm{Span}}}$ indicates the cross-polarized scattering power with respect to the total polarized power. Furthermore, using a simple relationship, we obtain:

$$\begin{array}{cc}\hfill tan{\theta }_{\mathrm{xP}}& =tan\left({\eta }_1-{\eta }_2\right)\hfill \\ \hfill & =\frac{m_{\mathrm{xP}}\mathrm{Span}\left(C_{11}-C_{22}\right)}{C_{11}{C}_{22}+m_{\mathrm{xP}}^2{\mathrm{Span}}^2}.\hfill \end{array}$$

(4)

where ${\theta }_{\mathrm{xP}}\in [-45,45]$.

The eigen-decomposition of ${\mathbf{C}}_2$ can be expressed as

$${\mathbf{C}}_2={\mathbf{U}}_2\Sigma {\mathbf{U}}_2^{-1},$$

(5)

where $\Sigma$ is a $2\times 2$ diagonal matrix with non-negative real elements, ${\lambda }_1\ge {\lambda }_2\ge 0$, which are the eigenvalues of ${\mathbf{C}}_2$, and ${\mathbf{U}}_2$ is $2\times 2$ unitary matrix, where $u_{\mathrm{i}}$’s are the unit orthogonal eigenvectors. We define the pseudo probabilities, $p_i$, in terms of the eigenvalues as

$$p_i=\frac{{\lambda }_i}{{\sum }_{k=1}^2{\lambda }_k},$$

(6)

which we then use to define the scattering entropy as

$$H_{\mathrm{xP}}=-\sum _{k=1}^2{p}_k{log}_2\left(p_k\right).$$

(7)

The degree of polarization $m_{\mathrm{xP}}$ characterizes the degree of coherence of partially polarized waves. In contrast, the entropy $H_{\mathrm{xP}}$ describes the degree of statistical disorder associated with partially polarized waves. The parameter ${\theta }_{\mathrm{xP}}$ describes the scattering characteristics of a target. Therefore, in this work, these parameters provide information about the growth stage or phenology of crops.

2.4. $H_{\mathrm{xP}}/{\theta }_{\mathrm{xP}}$ Bound

The feasible regions in the $H_{\mathrm{xP}}/{\theta }_{\mathrm{xP}}$ plot is represented by the bounding curves, Curve I (${\mathbf{C}}_{2\left(I\right)}$) and Curve II (${\mathbf{C}}_{2\left(II\right)}$) as shown in Figure 2 and is given as

$${\mathbf{C}}_{2\left(I\right)}={\displaystyle \frac{1}{1+k}}\left[\begin{array}{cc}1& 0\\ 0& k\end{array}\right],{\theta }_{\mathrm{xP}}^{\max }={tan}^{-1}\left({\displaystyle \frac{k^2-2k+1}{k^2-k+1}}\right)$$

(8)

$${\mathbf{C}}_{2\left(II\right)}={\displaystyle \frac{1}{k+1}}\left[\begin{array}{cc}k& 0\\ 0& 1\end{array}\right],{\theta }_{\mathrm{xP}}^{\min }=-{tan}^{-1}\left({\displaystyle \frac{k^2-2k+1}{k^2-k+1}}\right)$$

(9)

where $k\in [0,1]$ denotes the scattering amplitude ratio, and therefore, ${\theta }_{\mathrm{xP}}\in [-45,45]$. One can observe from Equation (4) that when:

• $m_{\mathrm{xP}}=0$ (i.e., when there exists no polarization structure in the scattered EM wave), then ${\theta }_{\mathrm{xP}}=0{}^{\circ }$ characterize random scattering from targets.
• $m_{\mathrm{xP}}=1$, and ${\theta }_{\mathrm{xP}}=45{}^{\circ }$, characterize coherent scattering from deterministic targets (i.e., trihedral or dihedral).
• $m_{\mathrm{xP}}=1$, and ${\theta }_{\mathrm{xP}}=-45{}^{\circ }$, characterize cross-polarized scattering from complex targets.

Therefore, ${\theta }_{\mathrm{xP}}\in [-45{}^{\circ },45{}^{\circ }]$ characterizes diverse scattering-type information in between these extreme cases, and hence it is suitable to characterize scattering-type information from various targets.

2.5. Example of Variation of $\widehat{\alpha }$ and ${\theta }_{\mathrm{xP}}$

In dual polarimetry, Cloude proposed a target characterization parameter, $\overline{\alpha }\in [0{}^{\circ },90{}^{\circ }]$ using the eigendecomposition approach [24]. In order to compare this with the proposed ${\theta }_{\mathrm{xP}}\in [-45{}^{\circ },45{}^{\circ }]$, we have rescaled $\overline{\alpha }$ to $\widehat{\alpha }=45°-\overline{\alpha }$, where both $\widehat{\alpha }$ and ${\theta }_{\mathrm{xP}}$ vary from −45${}^{°}$to 45${}^{°}$.

To show the efficacy of ${\theta }_{\mathrm{xP}}$ over $\widehat{\alpha }$, we have considered three different land cover types, i.e., urban, oriented urban, and vegetation from the VV-VH data extracted from a full-pol C-band RADARSAT-2 data. The variations of ${\theta }_{\mathrm{xP}}$ and $\widehat{\alpha }$ over these targets are shown in Figure 3. One can note from Figure 3 that the median value of $\widehat{\alpha }$ over orthogonal urban area is approximately 35°. In contrast, the median value of ${\theta }_{\mathrm{xP}}$ is approximately 42°. The degree of polarization over the area is approximately $0.91$. This high degree of polarization essentially confirms the nearly coherent scattering characteristics of the target. In this regard, ${\theta }_{\mathrm{xP}}$ better characterizes the target to be coherent than $\widehat{\alpha }$.

However, over oriented urban area, we observe a large variation of both ${\theta }_{\mathrm{xP}}$ and $\widehat{\alpha }$. This large variation might be due to the diversity in the polarization randomness due to the orientation of buildings to the radar line of sight and orientation in the azimuth direction. Over vegetation surfaces, we observe a partially polarized scattered wave that consists of a major coherent part and a significant incoherent part. The coherent part might be due to the direct scattering from the nearly smooth soil surface, while the incoherent part is due to the scattering from the volume media. In this case, ${\theta }_{\mathrm{xP}}$ also correctly characterizes the coherent contribution of the target during scattering as compared to $\widehat{\alpha }$. A $\widehat{\alpha }-{\theta }_{\mathrm{xP}}$ plot for RADARSAT-2 data over San-Francisco in VV-VH mode is shown in Appendix A.1. Therefore, from these analyses, we can observe the enhancement of target scattering characterization using ${\theta }_{\mathrm{xP}}$. Thus, ${\theta }_{\mathrm{xP}}$ can be further utilized to characterize different phenological stages of diverse crop types.

2.6. Unsupervised Clustering Zones over Vegetative Surface

In this study, we are interested in observing the changes in the scattering mechanisms over different phenological stages of crops. One should note that, for natural targets, $C_{11}\ge C_{22}$ almost always, and hence, ${\theta }_{\mathrm{xP}}\in [0{}^{\circ },45{}^{\circ }]$. Let us now characterize $m_{\mathrm{xP}}$ and ${\theta }_{\mathrm{xP}}$ for a few particular scattering scenarios:

• For a pure diffused target, $m_{\mathrm{xP}}=0$, implies, ${\theta }_{\mathrm{xP}}=0$.
• For pure or point scatterer, $m_{\mathrm{xP}}=1$ and ${\eta }_2=0$, implies ${\theta }_{\mathrm{xP}}=45$.
• Infeasible scattering: $C_{11}=C_{22}$ and $m_{\mathrm{xP}}=1$.

Therefore, in terms of vegetation development, for a highly random vegetative structure, $m_{\mathrm{xP}}\approx 0$ and ${\theta }_{\mathrm{xP}}\approx 0$. However, a slightly rough soil surface acts like a Bragg scatterer, and hence $m_{\mathrm{xP}}\approx 1$ and ${\theta }_{\mathrm{xP}}\approx 45$.

The variations of ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ over bare field and vegetative field types are shown in Figure 4. It can be noted that ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ are highly sensitive with the canopy structure. Interestingly, we observed that the sensitivity is much higher for the highly erectophile crop structures such as canola, corn, and wheat. In contrast, the sensitivity is low for the soybean crop. This high sensitivity of ${\theta }_{\mathrm{xP}}$ might be due to the volume decorrelation with an increased canopy structure.

In this study, we utilize the ${\theta }_{\mathrm{xP}}$ and $H_{\mathrm{xP}}$ parameters to propose a new unsupervised clustering scheme to describe the target scattering behavior. The overall clustering scheme consists of 12 zones: Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8, Z9, Z10, Z11, and Z12. These divisions of the $H_{\mathrm{xP}}/{\theta }_{\mathrm{xP}}$ plane are based on certain scattering symmetry assumptions. In particular, the scattered Stokes vector, ${\overrightarrow{\mathbf{S}}}_s$ is a function of the target property, represented in terms of the Kennaugh matrix, $\mathbf{K}$ and the transmitted Stokes vector, ${\overrightarrow{\mathbf{S}}}_t$. Moreover, the two Stokes vectors are related as ${\overrightarrow{\mathbf{S}}}_s=\mathbf{K}{\overrightarrow{\mathbf{S}}}_t$. Furthermore, one can note that ${\overrightarrow{\mathbf{S}}}_s$ is also a function of the scattering order (n) [50]. Thus, the definition of scattering entropy, H [51,52] can be expressed in terms of n. Following this, we can observe that H increases with increasing n. For example, $H=0$ for $n=0$; $H\approx 0.3$ for $n=1$; $H\approx 0.5$ for $n=2$, and, $H\approx 0.7$ for $n=3$. However, for $n>3$, change in H is insignificant. Therefore, following this trend in the variation of the entropy, we divided the radial axis $H_{\mathrm{xP}}$ into four sub-sections by considering the values: 0.3, 0.5, and 0.7. An illustration of the $H_{\mathrm{xP}}/{\theta }_{\mathrm{xP}}$ clustering plane with the zones is shown in Figure 5.

3. Results and Discussion

In this section, we analyze the temporal dynamics of crops using the proposed dualpolarimetric descriptors. Furthermore, we utilize the proposed clustering framework to assess the phenological stages of the two crops, i.e., canola and wheat, from the C-band Sentinel-1 dual-pol SLC SAR data.

3.1. Canola

The spatio-temporal changes of ${\theta }_{\mathrm{xP}}$ are shown in Figure 6. The temporal variations of ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ are shown in Figure 7 and Figure 10. For temporal analysis, we considered 48 sampling points from three fields (Field No. 206, 208 and 224). We plot the variations of the proposed SAR descriptors with respect to the phenological changes of canola at different dates.

According to the in-situ information, canola was sown at the end of May. Hence, during June, the phenological stage of canola was bounded within the early to mid vegetative stage. In this regard, it is worth mentioning that canola is a broad-leaf crop with comparable leaf size to the wavelength of C-band SAR (≈5.6 cm). Therefore, due to the formation of dense rosette near the soil surface, the leaf structure considerably affects the SAR backscatter values. Thus, due to the high dynamic variations of crop morphology, canola is of particular interest in SAR research.

From Figure 7, we can observe that, on 13 June, the median value of $m_{\mathrm{xP}}\approx 0.78$. This high value of $m_{\mathrm{xP}}$ typically suggests high polarized scattered components from canola fields. This trend might be due to the low vegetation density in the fields. From the in-situ information, we observe that the canola was in the stem elongation stage during this time. Consequently, low biomass and low Plant Area Index (PAI) are evident at this period. Moreover, the vertical crop structure slightly attenuated the transmitted vertical (V) polarization. Hence, we observe a median value of ${\theta }_{\mathrm{xP}}\approx 40°$. Apart from this, the low depolarization in the scattered wave might be due to the scattering from leaves and canopy structures combined.

Furthermore, Figure 8 provides information about the different morphological characteristics of canola at this particular period. We observe that the points from these canola fields are mainly distributed in two separate zones: Z4 and Z7 in Figure 8a. This information typically infers that, during 13 June, the crop density in some canola fields was low. As a result, low entropy pure scattering was evident from those fields (Z4). In contrast, other fields had comparatively dense foliage structures, which increased scattering randomness increased (Z7). From in-situ measurement, we observe that fields 206 and 208 had higher PAI and biomass values than field 224. In field 224, most canola crops were at the leaf development stage. In addition, we observe a point in the Z1 region and points in the Z10 zone, indicating more advanced phenological stages. These variations from Z1 to Z10 reflect the differences in sowing dates for canola, and, consequently, differences in crop emergence and development among fields early in the season.

On 7 July, canola crops reached their flowering stage. During this time, we can observe a sudden drop in the ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ values. The median value of $m_{\mathrm{xP}}\approx 0.44$. This sudden decrease in the scattered polarization structure is due to the development of flowers, stems, and branches. During this time, the high attenuation of the V-polarized transmitted wave also lowered the ${\theta }_{\mathrm{xP}}$ value. Furthermore, the generation of high cross-polarization due to the complex crop canopy structure also aided in lowering ${\theta }_{\mathrm{xP}}$ values. During this time, the median value of ${\theta }_{\mathrm{xP}}$ is approximately 20°.

Similar changes in the clustering zones are evident from Figure 8b. We observe a shift of clustering zone from Z4 and Z7 to Z11. The formation of the dense cluster in Z11 indicates high randomness in the scattered wave due to the complex flower and branch structures. Moreover, during the flowering and early pod development stages, the leaf density drastically reduces. Canola crops drop their leaves during this period of pod and seed development. However, this dense cluster in Z11 also suggests that the number of flowers per plant was high, and the health of the crop was good. On the other hand, we observe some points in the comparatively low entropy region, Z10.

On 19 July, small changes in ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ are observed as compared to 7 July. This effect might be due to the existence of a similar complex canopy structure due to which the depolarization is evident in the scattered wave. During this early to the mid-pod development stage, ramified stems and seeds appear. In the early stage, seeds remain watery and translucent, while seeds reach maximum dry weight at the advanced stage. Therefore, depending upon this particular phenomenon, changes in the scattering mechanisms might be observed. However, during the satellite acquisition, most fields reached the early to mid-stage of pod development. As a result, the variance in ${\theta }_{\mathrm{xP}}$ is low during this time. The median value of ${\theta }_{\mathrm{xP}}$ is ≈18°, while the median of $m_{\mathrm{xP}}$ is approximately $0.4$.

We also observe the accumulation of $H_{\mathrm{xP}}/{\theta }_{\mathrm{xP}}$ cluster in Z11 region in Figure 8c. During this time, the complex vegetation canopy structure increased the randomness in the scattered wave. Consequently, we observe a marginal increase in $H_{\mathrm{xP}}$ compared to 7 July. However, as stated earlier, most of the fields reached the early to mid-stage of pod development, and thus we observe the majority of the point cloud in the Z11 region.

Subsequently, on 24 August, we observe an increase in ${\theta }_{\mathrm{xP}}$ values. Similarly, the increase in $m_{\mathrm{xP}}$ suggests that the polarized component in the scattered wave has increased during this period. This increase in the values of ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ might be due to the changes in the canopy moisture content. However, during this period, high variance in the data is evident. During this period, the median value of ${\theta }_{\mathrm{xP}}$ is approximately 34°, while the median value of $m_{\mathrm{xP}}$ is approximately $0.58$. In Figure 8d, we observe three distinct zones; Z7, Z10, and Z11. This phenomenon is due to the harvest of canola crops. The post-harvest residue might have generated marginal depolarization in the scattered wave. The percentages of pixels at each date are shown in

Table 2. Temporal variation in the percentage of data points in each zone for different phenology stages of canola. The zones with high percentages of points at a particular phenology stage are highlighted in bold.

Dates	Z1	Z2	Z3	Z4	Z5	Z6	Z7	Z8	Z9	Z10	Z11	Z12
13 June	2.0	0.0	0.0	58.0	0.0	0.0	35.4	0.0	0.0	4.6	0.0	0.0
7 July	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	6.2	91.6	2.2
19 July	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	6.2	93.8	0.0
24 August	0.0	0.0	0.0	0.0	0.0	0.0	12.6	0.0	0.0	64.6	22.8	0.0

.

3.2. Wheat

In this section, we analyse the temporal variation of ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ over wheat fields. Wheat is an edible cereal grass with an erectophile canopy structure. Over the test site, sowing of wheat was performed during the beginning of May. This study has considered three wheat fields (Field No. 233, 220, 62) for the temporal evaluation. The temporal patterns of ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ are shown in Figure 10. Alongside this, we represent the dynamic behaviour of unsupervised clusters in Figure 9.

On 13 June, wheat crops reached the tillering stage. This stage starts with the appearance of tillers and ends with mature leaf sheaths. During this time, the crop density remains low. As a result, the scattered wave remains significantly polarized. It can be observed from Figure 10 that the median value of $m_{\mathrm{xP}}$ during this period is $0.78$. This essentially indicates that approximately 78% of the scattered wave is completely polarized, while approximately 22% is depolarized. This marginal depolarization effect might be due to the canopy interaction. Alongside this, the median value of ${\theta }_{\mathrm{xP}}$ is ≈39°, indicating that the target is a closely coherent scatterer. However, the high standard deviations in both ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ represent the variations in the uneven phenological advancement within the fields due to variations in seeding dates.

The clustering zones on 13 June are shown in Figure 9a. Like canola, we can observe that the clusters are distributed in four distinct zones; Z4, Z7, Z10, and Z11. The existence of these distinct zones is due to the differences in the morphological characteristics of wheat crops. From in-situ measurements, it is observed that the plant density of Field No. 220 was low (≈$100 {\mathrm{m}}^{-2}$). In contrast, the plant density over Field Nos. 233 and 62 was ≈$125 {\mathrm{m}}^{-2}$ and ≈$190 {\mathrm{m}}^{-2}$, respectively. Similarly, high PAI and biomass are evident over these two fields. As a result, a high attenuation of the transmitted wave is observed over Field No. 233 and 62. Hence, we observe clouds of clusters at Z10 and Z11. On the other hand, the attenuation of the scattered wave is comparatively lower over Field No. 220. Hence, cluster formation at Z4 and Z7 is observed in Figure 9a.

On 7 July, we observe a decrease in both ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ values due to the appearance of more branches and leaves. The median value of ${\theta }_{\mathrm{xP}}$ is around 27°, while the median value of $m_{\mathrm{xP}}$ is around $0.49$. At this stage, wheat has advanced to the early flowering stage. At this stage, flowers appear on the upper portion of the canopy layer. As a result, the dense wheat structure and the flowers increase the multiple interactions of the EM wave. Hence, a high proportion of depolarized components in the scattered wave is evident during this period. Similar changes in the clusters are evident from Figure 9b. At this period, the clusters shift to Z11 from Z4 and Z7 due to the increased scattering entropy due to the randomness in the scattered wave by the targets. However, a cloud of clusters is seen at Z10, which might be due to some wheat crops having late phenological development.

On 19 July, wheat crops advanced to the early dough stage. At this phenological period, wheat grains appear, and it remains milky. Therefore, the dielectric property of grains affects the scattering mechanism. Moreover, the denser wheat canopy increases multiple scattering. As a result, the values of $m_{\mathrm{xP}}$ reduce further as compared to 7 July. Similarly, the coherent component in the scattered EM wave also reduces. The median value of $m_{\mathrm{xP}}$ is $\approx 0.37$. The median value of ${\theta }_{\mathrm{xP}}$ is $\approx 20°$. The high variation in the values of ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ might be due to the randomly oriented wheat stems and heads. Similar scattering characteristics from the wheat during the dough stage are reported by Wu et al. [53]. In particular, during this period, the scattering phenomenon occurs primarily from two parts of the wheat canopy, i.e., thick upper canopy layer and relatively sparse lower canopy layer.

Further from Figure 9c, we can observe the formation of clusters in the Z11 region. However, we can also observe some points in Z12 due to additional complex scattering mechanisms from the wheat canopy subjected to partial lodging or related canopy anomalies. It is important to note that we observe much spread of points in the Z11 zone, unlike other phenological stages of canola and wheat. This suggests that, during this particular phenology period, the appearance of grains and their orientations might be responsible for breaking the reflection asymmetry condition during scattering from the wheat canopy. Moreover, as stated earlier, the complex canopy geometry, the spatial distribution of grains, and their random orientations increased the overall randomness in the scattered EM wave. As a result, for the majority of the points, we observe high entropy values.

During 24 August, the harvest of wheat started. At this stage, the interaction of EM waves mainly takes place with the soil surface and post-harvest residue. Consequently, a reverse trend is observed in ${\theta }_{\mathrm{xP}}$ and $m_{\mathrm{xP}}$ values. From Figure 10, we can observe that the median value of $m_{\mathrm{xP}}$ is approximately $0.58$, while the median value of ${\theta }_{\mathrm{xP}}$ is approximately 34. These values indicate the coherent component increase in the scattered EM wave from the wheat canopy. Additionally, Figure 9d shows three different clustering zones, Z7, Z10, and Z11, during this period. The percentages of temporal variations of clusters at different zones are shown in

Table 3. Temporal variation in the percentage of data points in each zone for different phenology stages of wheat. The zones with high percentages of points at a particular phenology stage are highlighted in bold.

Dates	Z1	Z2	Z3	Z4	Z5	Z6	Z7	Z8	Z9	Z10	Z11	Z12
13 June	0.0	0.0	0.0	27.0	0.0	0.0	47.9	0.0	0.0	10.4	14.7	0.0
7 July	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	18.8	81.2	0.0
19 July	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	83.3	16.7
24 August	0.0	0.0	0.0	0.0	0.0	0.0	31.3	0.0	0.0	29.1	39.6	0.0

. In addition, the temporal variations of ${\theta }_{\mathrm{xP}}$ and $\widehat{\alpha }$ are shown in Appendix A.2. One can observe that the temporal signature of these parameters differs for wheat and canola. Therefore, these descriptors can be efficiently utilized to distinguish wheat and canola fields using time series data. Furthermore, it might also be possible to distinguish these two crop fields from other land cover targets. The temporal variation of the clusters over the cultivation area is shown in Figure 11.

4. Conclusions

In this study, we propose a new polarimetric target characterization parameter, ${\theta }_{\mathrm{xP}}\in [-45{}^{\circ },45{}^{\circ }]$ for dual-pol Sentinel-1 Synthetic Aperture Radar (SAR) data. However, over natural surfaces, such as vegetation, ${\theta }_{\mathrm{xP}}\in [0{}^{\circ },45{}^{\circ }]$, where ${\theta }_{\mathrm{xP}}=0°$ denotes totally incoherent targets and ${\theta }_{\mathrm{xP}}=45°$ denotes totally coherent targets. Utilizing the scattering entropy parameter, $H_{\mathrm{xP}}$ along with ${\theta }_{\mathrm{xP}}$, we proposed a novel unsupervised classification scheme for dual-pol SAR data. We further utilized ${\theta }_{\mathrm{xP}}$ and the classification scheme to analyze the phenological development of canola and wheat over a Canadian test site. The results are promising and exhibit high sensitivity with the morphological changes of the crops at each phenological stage.

We observe a broad dynamic range of ${\theta }_{\mathrm{xP}}$ from the temporal analysis of both canola and wheat starting from their early vegetative stage to maturity and harvest stage. For both crops, the variation of ${\theta }_{\mathrm{xP}}$ within the desired phenology window is $\approx 40°$ to $\approx 15°$. High values of ${\theta }_{\mathrm{xP}}$ during the early stages are due to the low crop density within the fields. In contrast, the low ${\theta }_{\mathrm{xP}}$ values during the flowering stage are due to the complex canopy structure and additional appearance of flowers and grains. Moreover, we observe a trend reversal for both canola and wheat on 24 August. This trend reversal is due to the harvest of both crops during this period.

Within the scope of this study, we have characterized diverse crop phenological stages in terms of the physical scattering of the electromagnetic wave. The unsupervised classification scheme comprises twelve distinct zones, which represent these different physical scattering mechanisms. Thus, these clustering zones are beneficial for identifying the morphological status of the crop when a priori information is not available. Alongside this, the clustering scheme can also capture the difference in the growth stages simultaneously. Moreover, during the initial growth period, most of the points lay in the low entropy zones, while those points shifted towards high entropy zones at advanced phenological stages. This phenomenon is due to the changes in the canopy complexity with the advancement of the crop phenological stages. We already know that the dual-pol configuration does not have either the VV or HH polarization information. Hence, by utilizing ${\theta }_{\mathrm{xP}}$, we are unable to uniquely characterize scattering-type information, i.e., whether it is an even-bounce or odd-bounce mechanism. In addition, as stated earlier, if the SAR signal gets highly saturated due to the morphological characteristics of targets, ${\theta }_{\mathrm{xP}}$ might provide a similar saturation effect. This appearance is evident during the end of flowering to maturity stages of crops.

We can extend this study to several other crop types using different dual polarimetric SAR sensor configurations, such as HH-HV. The proposed parameter and the clustering framework will be valuable to analyze data from the upcoming dual-pol NASA-ISRO Synthetic Aperture Radar Mission (NISAR) and Copernicus Sentinel SAR missions for agricultural and environmental studies. The open-source code is available at: https://github.com/Subho07/dual_cross_pol_theta_clustering.