A Comparison of Alternative Approaches to Temporal Matching of Spectral Clusters in Post-Catastrophic Landscape Dynamics

1. Introduction

Assessing the spatiotemporal dynamics of landscape and vegetation cover is a necessary prerequisite for understanding the directions and rates of ecosystem transformation, identifying centres of degradation, detecting disruptions in successional processes, and predicting further changes in territorial structure [1,2]. The analysis of spatiotemporal dynamics in landscape and vegetation cover becomes especially important under conditions of ecological catastrophes [3], when natural and anthropogenically disturbed complexes undergo rapid, large-scale, and often nonlinear transformations affecting both species composition and vegetation productivity, as well as the spatial organisation of landscape cover [4]. In territories affected by ecological catastrophes, it is essential both to document the state of the area at a given moment in time and to track trajectories of change, their intensity, spatial heterogeneity, and temporal persistence [5]. Remote sensing is an effective tool for analysing the spatiotemporal dynamics of landscape and vegetation cover, because it enables the assessment of temporal changes in spatial patterns over large areas, ensures repeatability of observations, and provides a basis for comparing the state of landscape and vegetation cover before, during, and after the impact of destabilising factors [6].

The spatiotemporal dynamics of landscape and vegetation cover can be investigated using both continuous indicators, particularly spectral indices that reflect gradients of surface condition, productivity, or moisture, and classification methods, which enable the distinction of discrete cover types and the analysis of their spatial relationships [7,8]. While continuous indicators are well-suited to capturing gradual changes in the intensity of individual properties [7,9], classification-based approaches provide an interpretation of structural transformations within a territory through changes in land-cover categories [10,11]. Obtaining classifications for individual dates and ensuring their accurate comparability over time is critically important. This challenge is addressed through temporal matching, which harmonises classification results obtained at different observation dates based on phenological, seasonal, and image-acquisition conditions [12,13]. The necessity of temporal matching of land-surface classification results stems from the fact that remote sensing data reflect the class of an object and its state at a specific moment in the seasonal cycle [14]. The spectral responses exhibit significant variation across the distinct phenological phases within vegetation cover [15]. The same vegetation type may appear as different classes on different dates, while different vegetation types may, during certain periods of the year, become spectrally similar [8,16]. These effects are further compounded by variations in moisture conditions, inundation, weather, cloud cover, acquisition geometry, and sensor characteristics, leading to the accuracy of single-date classifications depending strongly on the acquisition date [17,18]. Maps produced for different seasons, years, or satellite platforms may differ not because of actual landscape change but because of differences in observation times, creating a risk of pseudo-change, reducing model transferability, and complicating the accurate comparison of spatiotemporal series [13,19,20]. The temporal harmonisation of land-surface classification results is therefore a necessary condition for reliable monitoring, the detection of genuine cover changes, and the construction of comparable multi-year classifications [12,21].

During the temporal matching of land-surface classification results, several challenges are typically encountered. The phenological non-stationarity of vegetation types means that the same vegetation type may exhibit different spectral responses at different phases of the seasonal cycle. Classes that are well separated on one date may overlap on another. In many studies, accounting for phenological dynamics is regarded as an essential component of matching multi-date classifications [21,22,23]. Even in the absence of real change, spectral differences between scenes may arise from atmospheric conditions, illumination, acquisition geometry, and sensor differences. Consequently, inter-date differences often reflect variations in observation conditions rather than actual cover transformation, making reflectance normalisation a critical step before matching [24,25]. Even slight misregistration between scenes alters the configuration of object boundaries and may generate false changes, especially along the contours of fine-grained or mosaic landscapes. Although object-based approaches are less sensitive to misregistration, precise co-registration remains essential [26,27]. Reference samples are often derived from field descriptions, thematic maps, or aerial photographs acquired on different dates, under different phenological conditions, or at different spatial resolutions. This introduces label noise into the training data, associated with misregistration, delayed updating, and the intrinsic complexity of the cover itself [28]. Such errors affect both classifier training and validation procedures [29]. Errors also arise from mismatches in resolution and semantic inconsistencies when older, coarser-resolution maps are used for matching or as training labels for newer, more detailed products. This issue is particularly important when comparing classifications produced by different sensors or by different generations of mapping products [30]. A further challenge is the methodological incomparability of classifications produced in different years: changes in interpretation rules, spatial detail, feature sets, processing algorithms, and legend structure may generate discrepancies between maps that are unrelated to real cover change. Consequently, part of the detected dynamics may represent an artefact of differences in mapping procedures rather than actual transformation of landscapes or vegetation [31].

Temporal Cluster Matching (TCM) has been proposed for temporal matching and change detection of structures based on changes in the spectral relationship between an object footprint and its surrounding neighbourhood [32]. In object-based pasture mapping, temporal correspondence can be established by matching pasture objects according to the similarity of their multivariate Sentinel-2 time-series trajectories, particularly using Dynamic Time Warping to account for shifts in phenological timing [33]. Matching approaches have also been used to analyse spatial correspondence between pixel-based remote-sensing classifications and polygon-based forest inventory maps [34]. Existing matching approaches in remote sensing include temporal matching of spectral trajectories, object-based correspondence analysis, and hierarchical reconciliation between heterogeneous land-cover representations [35]. In forest monitoring, temporal matching can be formulated as the correspondence between abrupt shifts in NDVI time series and specific management events, such as clear-cuts or partial forest cuts [36]. For dynamic landscape monitoring, matching may be based on repeated image-pattern correspondence, allowing displacement or trajectory information to be extracted from temporally separated satellite scenes [37]. The use of multi-temporal Sentinel-2 imagery improved vegetation classification accuracy because seasonal phenological changes provided additional discriminatory information between spectrally similar vegetation types [38]. Thus, despite the broad use of matching in remote sensing, there remains a methodological gap in temporally linking independently derived vegetation or habitat clusters across multitemporal image series under phenological and interannual spectral drift. This is the problem addressed in the present study.

Phenological and successional dynamics operate on different temporal scales under normal conditions [39]. Phenological changes reflect intra-annual variability in the state of vegetation cover, whereas successional changes characterise the longer-term reorganisation of its structure and composition [40]. However, under catastrophic impacts, the issue of temporal alignment becomes considerably more complex, as these two types of dynamics begin to overlap [38] partially. Following environmental disturbance, phenological variability may manifest as quantitative changes in the spatial proportions of qualitatively similar cover types, where their density, productivity, or degree of development changes without a distinct transition to another type [41]. Catastrophic processes may sharply accelerate successional transformations, so that a qualitative reorganisation of the cover accompanies quantitative changes–that is, a transition to other vegetation types or forms of landscape organisation [42]. Under conditions of ecological catastrophe, interpreting remote-sensing time series requires distinguishing short-term phenological variability from accelerated successional change, because spectrally or structurally similar shifts may have fundamentally different ecological origins.

Our hypothesis is that temporal matching between spectral clusters representing the same landscape-cover type can be established from the internal geometry of spectral point clouds, particularly from their elongation and principal-axis orientation in spectral space. This hypothesis is based on the fact that, even within a single observation date, individual patches of the same landscape-cover type may occupy different phenological or successional states. Such within-type heterogeneity produces elongated spectral point clouds whose dominant orientation reflects the principal direction of internal spectral variability. Through time, successive states of the same landscape-cover type form a continuous spectral trajectory that, because of the discrete character of remote-sensing observations, appears as a sequence of partially separated elongated clouds. Under these conditions, the local geometry of a spectral point cloud constrains the probable direction of its subsequent displacement in spectral space and therefore provides a basis for temporal matching between successive states of the same landscape-cover type. This geometrically oriented approach is then compared with alternative matching methods, including centroid distance, Mahalanobis distance, Linear Discriminant Analysis, and Random Forest, in order to evaluate whether the geometry of spectral point clouds provides additional information for preserving temporal correspondence among landscape-cover states.

2. Materials and Methods

2.1. General Study Design and Strategy for Evaluating Temporal Matching Methods

The study aimed to evaluate the extent to which different temporal-matching methods can preserve the identity of habitat types across a multi-date series of satellite images. The overall analytical framework combined two complementary evaluation modes: a quasi-supervised mode, in which habitat types were defined based on reference geobotanical surveys, and an unsupervised temporal matching mode, in which the ability of the methods to maintain a stable cluster structure without direct use of external labels for each date was assessed. In the quasi-supervised mode, reference habitat types were identified for a fixed period of field surveys, after which the quality of their recognition was evaluated when transferring these correspondences to other dates within the satellite time series. This approach enabled assessment of the “ageing” effect of reference information, i.e. the decline in the reliability of correspondences with increasing phenological or interannual divergence from the survey date. In the unsupervised mode, methods were compared based on how consistently they reproduced the temporal continuity of clusters across consecutive dates. This evaluation included the emergence of new codes, the stability of chains of correspondences, the structure of cluster-transition graphs, and the consistency of the resulting superclusters among alternative approaches.

2.2. Study Area as a Model Landscape System

The floodplain ecosystems of Khortytsia Island on the Dnipro River were selected as the model study area, as they represent a dynamic landscape system that has undergone pronounced post-catastrophic restructuring following the destruction of the Kakhovka Reservoir in 2023 [43]. Khortytsia is the largest island on the Dnipro River, located within the city of Zaporizhzhia, Ukraine, downstream of the Dnipro Hydroelectric Power Plant [44]. Before the catastrophe, the island was connected with the northern part of the former Kakhovka Reservoir water area and its floodplain–riparian zone [45]. The study covered the period from 2021 to the most recent available satellite observations, which made it possible to compare the state of the landscape and vegetation cover before the catastrophe, during the phase of abrupt hydrological restructuring, and throughout subsequent post-catastrophic development [43]. Following the destruction of the Kakhovka Hydroelectric Power Plant dam on 6 June 2023, the rapid decline in water level caused substantial changes in the hydrological regime, shoreline morphology, moisture conditions of floodplain soils, and the spatial structure of habitats [43,44,46,47]. Within the floodplain system of Khortytsia, these processes were manifested in the disconnection of some water bodies from the main channel, shallowing, exposure of bottom substrates, formation of new sandy shores and shoals, and rapid colonisation of newly exposed surfaces by pioneer and riparian-aquatic vegetation. At the same time, changes in relatively elevated parts of the floodplain may proceed more slowly, but remain ecologically significant because of changes in the moisture regime of the root zone and increased summer stress on woody stands [48,49]. Khortytsia is therefore a suitable model site for analysing the temporal dynamics of landscape cover, as it combines sharply transformed, transitional, and relatively stable areas. This mosaic structure provides conditions for testing remote-sensing methods capable of tracking both gradual changes in known habitat types and the emergence of new, previously unrepresented cover classes [50].

2.3. Acquisition and Primary Processing of Satellite Imagery

A multi-date series of Sentinel-2 MSI (Multispectral Instrument) Level-2A satellite images acquired during 2021–2026 was used for the analysis [51]. Only scenes with low cloud contamination were included in the processing, with total cloud cover not exceeding 5%. Each scene was processed separately within the study-area polygon. Spectral bands were extracted from SAFE archives, clipped to the study area, masked by the polygon boundary, and aligned to a common spatial grid with an identical coordinate reference system, spatial extent, and resolution. The radiometric harmonisation of the imagery was performed, including BRDF (Bidirectional Reflectance Distribution Function) correction to reduce inter-date differences associated with varying acquisition conditions, aimed at reducing the influence of sun–sensor geometry on surface reflectance [52]. Cloud- and shadow-affected pixels were subsequently masked using Sentinel-2 quality layers. To preserve the continuity of the multi-date series, missing values were later reconstructed along the temporal trajectory of each pixel using regression-interpolation gap filling based on neighbouring observation dates [53,54]. Based on the spectral bands, 36 spectral indices were calculated for each date. These indices characterised different aspects of surface condition, including greenness, moisture, vegetation structure, hydrological regime, soil exposure, and integrated spectral brightness components. The formulas of the indices and their ecological interpretation are provided in our dataset [55]. All index layers were geometrically harmonised across dates and aligned to a common raster structure.

2.4. Identification of EUNIS (European Nature Information System) Habitats Based on Plant Communities

Geobotanical surveys, vegetation physiognomy, and the landscape position of patches were used for the interpretation of spectral classes in terms of EUNIS (European Nature Information System) habitats. Habitat identification was primarily based on plant associations established from field relevés, whereas vegetation physiognomy, moisture conditions, relief, and spatial context were used as complementary criteria. This approach follows the general logic of the EUNIS classification system, in which habitat types are associated with vegetation but are not direct equivalents of phytosociological syntaxa, as their delimitation also depends on ecological and geographical context [56]. The assignment of spectral classes to habitat types was performed by analysing the dominant and accompanying plant associations recorded within each class. This allowed avoiding a rigid “one class—one habitat” correspondence and retaining several spectrally distinct classes within the same EUNIS habitat type when they represented different structural or phenological states of the same habitat. The relationship between spectral classes and plant associations was evaluated using contingency-table analysis [57]. A contingency matrix was constructed in which rows corresponded to plant associations and columns to spectral classes. The overall strength of association was assessed using Pearson’s χ² test with a Monte Carlo permutation procedure because of the presence of low expected frequencies in part of the matrix [58,59]. The strength of association was quantified using Cramér’s V. Individual classes were additionally evaluated using one-versus-rest comparisons with Holm correction for multiple testing [60,61]. Standardised residuals were used to identify associations that were statistically overrepresented or underrepresented within particular spectral classes [62,63]. The detailed procedure of vegetation classification, including the diagnostic species composition of the identified floristic groups and their assignment to EUNIS habitat types, is presented in Supplementary Material S1.

2.5. Feature space for Temporal Matching

Calculations were performed in a common ordination space constructed for the full set of observations from all dates, using spectral indices shared by all temporal slices. For each spectral index $x_{ij}$, the seasonal component was described using a cyclic GAM:

$$x_{ij}=s_j\left({\mathrm{D}\mathrm{O}\mathrm{Y}}_i\right)+{\epsilon }_{ij},$$

where ${\mathrm{D}\mathrm{O}\mathrm{Y}}_i$is the day of year for the $i$-th observation, $s_j(\cdot )$is a cyclic smooth function, and ${\epsilon }_{ij}$is the residual component. The residuals ${\epsilon }_{ij}$were used as seasonally adjusted values. For each date $d$, the robust centre of feature $j$was then calculated as either the median or the trimmed mean:

$$c_{jd}=\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}{x}_{ij}:i\in d,c_{jd}={\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}_{trim}{x}_{ij}:i\in d.$$

The values were subsequently centred as:

$$x_{ij}^{\mathrm{\text{'}}}=x_{ij}-c_{jd}.$$

Within each “class × date” combination, quantile trimming was then applied. Observations were retained if:

$$Q_{q,j}^{\left(g,d\right)}\le x_{ij}^{\mathrm{\text{'}}}\le Q_{1-q,j}^{\left(g,d\right)},$$

where $Q_{q,j}^{\left(g,d\right)}$and $Q_{1-q,j}^{\left(g,d\right)}$ are the lower and upper quantiles for feature $j$in group $g$on date $d$, respectively. In the analysis, $q=0.01$was used. This sequence of transformations produced a comparable set of features for constructing the common matching space.

The prepared observations from all dates were analysed in a unified ordination space using a two-step principal component approach. First, principal component analysis of the total spectral variability (hereafter, total PCA) was applied to the complete multi-date dataset. This step extracted the dominant gradients shared across all dates, which for vegetation cover mainly reflected broad variation in greenness, moisture conditions, and surface brightness. In the second step, principal component analysis of residual variation (hereafter, residual PCA, or rPCA) was applied after removing the variation explained by the dominant total PCA gradients. This step was used because the strongest common gradients may mask weaker, but ecologically meaningful, differences among habitat states. The residual components, therefore, represented supplementary spectral structure that was less expressed in the main axes but could be important for distinguishing temporal trajectories and spectrally similar habitat states. Each observation was then represented by its coordinates in this combined component space. All component coordinates were used in the matching procedures, while the first two total PCA coordinates were also used to visualise spectral clouds and temporal trajectories.

2.6. Geometric Properties of Spectral Clusters

The geometric approach was based on the assumption that the shape of a spectral cluster reflects the temporal trajectory of change in the corresponding cover type. Accordingly, the geometry of the cluster can be interpreted as an indicator of the direction of subsequent spectral change and used to establish correspondences between clusters on consecutive dates. Temporal matching was therefore based on the shapes of spectral clusters, the orientations of their principal axes, and their relative positions in the common ordination space (Figure 1).

For each “class $k$× date $t$” combination, the spectral cluster in the common ordination space was described by a robust geometric descriptor:

$$D_{k,t}=\{{\mu }_{k,t},{\mathrm{u}}_{k,t},{\lambda }_{1,k,t},{\lambda }_{2,k,t},{\kappa }_{k,t}\},$$

where ${\mu }_{k,t}$ is the robust centroid of the cluster, $u_{k,t}$ is the unit vector of the first principal axis, ${\lambda }_{1,k,t}$ and ${\lambda }_{2,k,t}$ are the first and second eigenvalues of the robust covariance matrix, and:

$${\kappa }_{k,t}={\displaystyle \frac{{\lambda }_{1,k,t}}{{\lambda }_{2,k,t}}}$$

is the elongation, or anisotropy, index of the cluster.

The centre and covariance structure were estimated robustly using the minimum covariance determinant method, followed by eigen decomposition of the robust covariance matrix. The first eigenvalue was interpreted as longitudinal variability along the dominant direction of the cluster, whereas the second eigenvalue represented transverse variability. Thus, $u_{k,t}$ defined the orientation of the main gradient of within-cluster variability, while ${\kappa }_{k,t}$ described the degree of cluster elongation.

Correspondence between clusters on two consecutive dates was established asymmetrically. For each cluster on the current date $t$, its descriptor $D_{k,t}$ was compared with the descriptors of all clusters $D_{j,t-1}$ present on the previous date $t-1$. For a pair consisting of reference cluster $j$ at $t-1$ and current cluster $k$ at $t$, the displacement vector between centroids was first calculated as:

$${\delta }_{j\to k}={\mu }_{k,t}-{\mu }_{j,t-1}.$$

Its projection onto the principal axis of the reference cluster was defined as:

$$s_{j\to k}={\delta }_{j\to k}^{\mathsf{T}}{\mathrm{u}}_{j,t-1},$$

whereas the perpendicular component was calculated as:

$${\delta }_{j\to k}^{\perp }={\delta }_{j\to k}-s_{j\to k}{\mathrm{u}}_{j,t-1}.$$

The transverse separation between clusters was then described by the Euclidean length of this perpendicular component:

$$d_{j\to k}^{\perp }=\parallel {\delta }_{j\to k}^{\perp }\parallel .$$

Thus, $s_{j\to k}$ represents the longitudinal displacement along the principal axis of the reference cluster, whereas $d_{j\to k}^{\perp }$ represents the transverse deviation of the centroid of the current cluster from this axis.

For geometric interpretation, an axial “tube” region was defined around the principal axis of each cluster. For the reference cluster, its transverse radius was defined as:

$$r_{j\to k}^{\mathrm{r}\mathrm{e}\mathrm{f}}=\alpha \sqrt{{\lambda }_{2,j,t-1}}+\beta \mid s_{j\to k}\mid ,$$

where $\alpha$is the radius scaling coefficient and $\beta$is the axial expansion parameter. For the current cluster, the transverse radius was defined as:

$$r_{k,t}^{\mathrm{c}\mathrm{u}\mathrm{r}}=\alpha \sqrt{{\lambda }_{2,k,t}}.$$

In the calculations, $\alpha =1.5$ and $\beta =0$ were used. Under this setting, the axial region was effectively cylindrical, i.e. its transverse radius did not depend on longitudinal displacement.

The unnormalised transverse gap between two clusters was defined as:

$$g_{j\to k}=d_{j\to k}^{\perp }-\left(r_{j\to k}^{\mathrm{r}\mathrm{e}\mathrm{f}}+r_{k,t}^{\mathrm{c}\mathrm{u}\mathrm{r}}\right),$$

and the normalised transverse gap as:

$$g_{j\to k}^{*}={\displaystyle \frac{g_{j\to k}}{\mathrm{m}\mathrm{a}\mathrm{x}\left(r_{j\to k}^{\mathrm{r}\mathrm{e}\mathrm{f}}+r_{k,t}^{\mathrm{c}\mathrm{u}\mathrm{r}},\epsilon \right)}},$$

where $\epsilon$is a small positive constant introduced to avoid division by zero. If

$$g_{j\to k}\le 0,$$

The clusters were considered to overlap or touch geometrically in the transverse direction. Accordingly, negative values of $g_{j\to k}^{*}$indicate overlap of the axial regions, whereas positive values indicate a spatial gap between them. This value was subsequently used as the normalised gap (gap norm).

The angular difference between the orientations of two clusters was evaluated using their first principal axes:

$${\theta }_{j\to k}={\displaystyle \frac{2}{\pi }}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mid {\mathrm{u}}_{j,t-1}^{\mathsf{T}}{\mathrm{u}}_{k,t}\mid \right).$$

Because the scalar product is taken in absolute value, this measure is independent of the sign of the eigenvector. The value of ${\theta }_{j\to k}$ ranges from 0 to 1, where 0 indicates complete alignment of orientations and larger values indicate stronger angular mismatch. This measure was subsequently used as the angle between clusters.

The integrated geometric similarity score between clusters was defined as:

$$S_{j\to k}=w_{\theta }{\theta }_{j\to k}+w_g{g}_{j\to k}^{*},$$

where $w_{\theta }$and $w_g$are the weights of the angular and spatial components, respectively.

In the analysis, equal weights were used:

$$w_{\theta }=1,w_g=1.$$

Smaller values of $S_{j\to k}$ corresponded to greater temporal geometric similarity between clusters. For each cluster $k$ on date $t$, the best match on the previous date was defined as the cluster $j$ for which $S_{j\to k}$was minimal:

$$\widehat{j}(k,t)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{j\in {\mathcal{C}}_{t-1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{S}_{j\to k},$$

where ${\mathcal{C}}_{t-1}$ is the set of all clusters present on the previous date. This integrated score was subsequently used as the combined score.

A separate test was performed to determine whether the geometric relationship between temporally adjacent states of the same class was stronger than that between the current cluster and its nearest between-class alternatives. For each within-class temporal pair $\left(k,t-1)\to (k,t\right)$, the values ${\theta }_{k\to k}$, $g_{k\to k}^{*}$, and $S_{k\to k}$ were calculated and compared with the best between-class alternatives on the previous date:

$${\theta }_{\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k,t\right)=\underset{j\ne k}{\mathrm{m}\mathrm{i}\mathrm{n}}{\theta }_{j\to k},$$

$$g_{\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}}^{*,\mathrm{m}\mathrm{i}\mathrm{n}}\left(k,t\right)=\underset{j\ne k}{\mathrm{m}\mathrm{i}\mathrm{n}}{g}_{j\to k}^{*},$$

$$S_{\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k,t\right)=\underset{j\ne k}{\mathrm{m}\mathrm{i}\mathrm{n}}{S}_{j\to k}.$$

For each within-class pair, the following logical criteria were then recorded:

$${\mathrm{p}\mathrm{a}\mathrm{s}\mathrm{s}}_{\theta }=I\left({\theta }_{k\to k}<{\theta }_{\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}}^{\mathrm{m}\mathrm{i}\mathrm{n}}(k,t)\right),$$

$${\mathrm{p}\mathrm{a}\mathrm{s}\mathrm{s}}_g=I\left(g_{k\to k}^{*}<g_{\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}}^{*,\mathrm{m}\mathrm{i}\mathrm{n}}(k,t)\right),$$

$${\mathrm{p}\mathrm{a}\mathrm{s}\mathrm{s}}_S=I\left(S_{k\to k}<S_{\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{e}\mathrm{e}\mathrm{n}}^{\mathrm{m}\mathrm{i}\mathrm{n}}(k,t)\right),$$

where $I(\cdot )$is the indicator function. The proportions of cases in which the within-class pair had a smaller angle, a smaller normalised gap, or a smaller total geometric score than the nearest between-class alternative were then calculated. In the results, these values were interpreted as the proportions of successful matches according to the corresponding criterion.

The separation margin between the true within-class correspondence and the nearest external class was described as:

$$margi{n}_{score\hspace{0.17em}k,t}=S_{between}^{min}\left(k,t\right)-S_{k\to k}.$$

Positive values of $\mathrm{m}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{i}\mathrm{n}\_\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{i}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{e}$ that the true within-class correspondence has a smaller total geometric score than any between-class alternative and is therefore geometrically preferable. Negative values indicate cases in which at least one external class is geometrically closer to the current cluster than its own previous state. This value was subsequently reported as the margin score.

The statistical evaluation of the geometric criterion of temporal correspondence was performed nonparametrically. For $\theta$, $g^{*}$, and $S$, within-class values were compared with the corresponding best between-class alternatives using the paired Wilcoxon test with a one-sided alternative, according to which within-class values were expected to be smaller. For the $margi{n}_{score}$, a one-sample Wilcoxon test against zero was applied with the alternative hypothesis of a positive shift. In addition, the proportion of within-class pairs for which the axial regions overlapped, i.e. $g_{k\to k}\le 0$, was calculated and used as an indicator of geometric connectivity between consecutive temporal states of the same class.

2.7. Axis-Distance and Geometric Methods

The axis-distance and geometric matching methods were based on the assumption that a cluster in spectral space should be considered not only as a set of points or as a centroid, but also as an oriented geometric structure described by the centre of the cluster, the direction of its principal axis, and its transverse width. The following cluster descriptor was used as the basic reference for comparison:

$$D_{j,t-1}=\{{\mu }_{j,t-1},{\mathrm{u}}_{j,t-1},{\lambda }_{1,j,t-1},{\lambda }_{2,j,t-1}\},$$

where ${\mu }_{j,t-1}$is the centroid of the reference cluster, $u_{j,t-1}$is the unit vector of its first principal axis, and ${\lambda }_{1,j,t-1}$and ${\lambda }_{2,j,t-1}$are the first and second eigenvalues of the robust covariance matrix. This description made it possible to move from point-based to axis-based comparison of clusters between consecutive dates.

In the simplest centroid-based variant, correspondence between clusters was established only by the Euclidean distance between their centroids:

$$d_{j\to k}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}}=\parallel {\mu }_{k,t}-{\mu }_{j,t-1}\parallel .$$

The best match for cluster $k$ on date $t$was then defined as:

$${\widehat{j}}_{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}}(k,t)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{j\in {\mathcal{C}}_{t-1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{d}_{j\to k}^{\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}},$$

where ${\mathcal{C}}_{t-1}$ is the set of all clusters present on the previous date. In this approach, the reference for comparison is point-based and reflects only the current position of the cluster in feature space.

The axis-distance method represents a geometric extension of this approach. Here, the reference for comparison is not a single point but the principal axis of the reference cluster, defined by the line:

$$L_{j,t-1}\left(\tau \right)={\mu }_{j,t-1}+\tau {\mathrm{u}}_{j,t-1},\tau \in \mathbb{R}.$$

Unlike a centroid, this axis may pass both through and beyond the point cloud; it therefore represents not only the current state of the system but also the dominant direction of its internal variability, which can be interpreted as a hypothetical trajectory of subsequent spectral change. For a pair consisting of reference cluster $j$at $t-1$and current cluster $k$at $t$, the centroid displacement vector was calculated as:

$${\delta }_{j\to k}={\mu }_{k,t}-{\mu }_{j,t-1},$$

with its longitudinal projection onto the principal axis of the reference cluster given by:

$$s_{j\to k}={\delta }_{j\to k}^{\mathsf{T}}{\mathrm{u}}_{j,t-1},$$

and its perpendicular component by:

$${\delta }_{j\to k}^{\perp }={\delta }_{j\to k}-s_{j\to k}{\mathrm{u}}_{j,t-1}.$$

The axis-distance score was defined as the transverse distance from the centroid of the current cluster to the principal axis of the reference cluster:

$$d_{j\to k}^{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}}=\parallel {\delta }_{j\to k}^{\perp }\parallel .$$

Accordingly, the best match was defined as:

$${\widehat{j}}_{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}}\left(k,t\right)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{j\in {\mathcal{C}}_{t-1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{d}_{j\to k}^{\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{s}}.$$

Thus, whereas the centroid method minimises the full distance between cluster centres, the axis-distance method minimises the transverse deviation from the axial direction of the reference cluster. This means that clusters are considered close primarily when their new state is positioned along the expected axis of change.

The geometric method used the same axial principle, but complemented it by assessing the angular alignment of the principal axes and the transverse gap between the axial regions of the two clusters. For this purpose, a transverse radius was defined around the principal axis of the reference cluster as:

$$r_{j\to k}^{\mathrm{r}\mathrm{e}\mathrm{f}}=\alpha \sqrt{{\lambda }_{2,j,t-1}}+\beta \mid s_{j\to k}\mid ,$$

and, for the current cluster, as:

$$r_{k,t}^{\mathrm{c}\mathrm{u}\mathrm{r}}=\alpha \sqrt{{\lambda }_{2,k,t}},$$

where $\alpha$is the scale of the transverse radius and $\beta$ is the parameter controlling expansion along the axis. In the calculations, $\alpha =1.5$and $\beta =0$ were used; therefore, the axial region was effectively cylindrical.

The transverse gap between the two clusters was defined as:

$$g_{j\to k}=\parallel {\delta }_{j\to k}^{\perp }\parallel -\left(r_{j\to k}^{\mathrm{r}\mathrm{e}\mathrm{f}}+r_{k,t}^{\mathrm{c}\mathrm{u}\mathrm{r}}\right),$$

and its normalised form as:

$$g_{j\to k}^{*}={\displaystyle \frac{g_{j\to k}}{\mathrm{m}\mathrm{a}\mathrm{x}\left(r_{j\to k}^{\mathrm{r}\mathrm{e}\mathrm{f}}+r_{k,t}^{\mathrm{c}\mathrm{u}\mathrm{r}},\epsilon \right)}},$$

where $\epsilon$is a small positive constant. If $g_{j\to k}\le 0$, the axial regions of the two clusters were considered to overlap or touch in the transverse direction.

The angular mismatch between the principal axes of the two clusters was calculated as:

$${\theta }_{j\to k}={\displaystyle \frac{2}{\pi }}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mid {\mathrm{u}}_{j,t-1}^{\mathsf{T}}{\mathrm{u}}_{k,t}\mid \right),$$

a quantity ranging from 0 to 1 and independent of the sign of the eigenvector. The integrated geometric score was then defined as:

$$S_{j\to k}=w_{\theta }{\theta }_{j\to k}+w_g{g}_{j\to k}^{*},$$

where $w_{\theta }$and $w_g$are the weighting coefficients of the angular and spatial components, respectively. In the calculations, $w_{\theta }=1$ and $w_g=1$ were used. The best match for cluster $k$ on date $t$was defined as:

$${\widehat{j}}_{\mathrm{g}\mathrm{e}\mathrm{o}\mathrm{m}}\left(k,t\right)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{j\in {\mathcal{C}}_{t-1}}{\mathrm{m}\mathrm{i}\mathrm{n}}{S}_{j\to k}.$$

Thus, unlike the axis-distance method, which considers only the transverse deviation from the principal axis of the reference cluster, the geometric method uses a combined criterion that simultaneously evaluates axial proximity, angular alignment, and the degree of transverse overlap between axial regions. The axis-distance and geometric methods were therefore treated as two distinct but hierarchically related approaches: the former as a minimal axial alternative to the centroid method, and the latter as its more complete geometric generalisation.

2.8. Alternative Matching Approaches

For comparison with the geometric method, four alternative approaches to cluster matching between consecutive dates were used: centroid distance, Mahalanobis distance, Linear Discriminant Analysis (LDA), and Random Forest (RF) [32]. All methods evaluated how a cluster on date $t_2$ corresponded to the classes defined on the preceding date $t_1$, after which the class with the best score was selected. In the centroid method, the Euclidean distance was calculated between the centroid $\mu$of the cluster on date $t_2$ and the centroids of all reference clusters on date $t_1$. The class with the minimum distance was considered the best match. The Mahalanobis method used both the centroid and the cluster’s covariance structure. For each class on the reference date, a distributional model was constructed, and for the query cluster, the median Mahalanobis distance of its points to this model was calculated. The class with the minimum distance was then selected as the best match. Linear Discriminant Analysis (LDA) and Random Forest (RF) were used as supervised classification-based alternatives. These models were trained on points from the ref. For the points belonging to a cluster on date(lower index 2) $t_2$, posterior probabilities of membership in each reference class were then obtained. These probabilities were averaged across all points within the query cluster, and the class with the highest mean probability was taken as the corresponding class. The comparison involved geometric axis-tube matching against centroid distance, Mahalanobis distance, LDA, and RF.

2.9. Strategy for Comparing Temporal Matching Methods for Spectral Clusters

The comparison of temporal matching methods for spectral clusters was performed in three consecutive stages (Table 1). In the first stage, the geometric assumption underlying the axis-distance and geometric approaches was tested. Whether point clouds of the same class on neighbouring dates preserve geometric connectivity, that is, whether their orientational and spatial similarity is greater than that of the nearest between-class alternatives. To this end, the relative position, direction of the principal axis, and degree of overlap of spectral clusters were evaluated in the common ordination space. In the second stage, the practical quality of matching between neighbouring dates was assessed. The axis-distance and geometric approaches were compared with alternative methods based on centroid distance, Mahalanobis distance, Linear Discriminant Analysis, and Random Forest. The comparison was performed in a quasi-supervised mode, in which habitat types were defined from reference geobotanical surveys conducted at a fixed date, after which the transferability of these correspondences to other temporal slices was analysed. In the third stage, the temporal stability of the methods was analysed, that is, changes in matching quality with increasing phenological and interannual drift relative to the reference period. The analysis evaluated how rapidly identification accuracy degraded with increasing temporal distance and whether geometrically oriented approaches were better able to preserve correspondence with the true class under conditions of post-catastrophic landscape cover restructuring. The evaluation scheme comprised three interconnected levels: testing the geometric hypothesis, benchmark comparison of practical matching quality, and analysis of the temporal stability of alternative methods.

2.10. Criteria for Comparing Alternative Matching Methods

The comparison of cluster-matching approaches between neighbouring dates was performed in a benchmark mode, in which, for each cluster on the current date, $t$, its correspondence to all candidate clusters presents on the previous date $t-1$was evaluated. The true correspondence was defined as the cluster belonging to the same reference habitat class on the previous date. For each matching case, the predicted class, the rank of the true class among all candidates ($rank\_true$), the $correct$indicator for a correct top-ranked correspondence, and the $margin$representing the separation between the true class and the nearest alternative candidate were determined. The overall performance of the methods was evaluated using top-1 and top-3 accuracy, the mean and median ranks of the true class, and the mean and median values of m a. r g i. n. These metrics were calculated for the full set of transitions between neighbouring dates, as well as separately for each date pair and for each class. For methods based on distance or geometric scores, the best correspondence was defined as the candidate with the minimum score, whereas for Linear Discriminant Analysis (LDA) and Random Forest, the best correspondence was defined as the class with the highest mean posterior probability.

2.11. Evaluation of the Sensitivity of Matching Methods to Phenological and Interannual Temporal Drift

Cluster matching across consecutive dates was considered a problem of establishing temporal invariance of vegetation types relative to the reference surveys conducted in August 2024 and 2025. The rationale for this approach was that a reference survey captures the state of vegetation cover at a specific point in time, whereas before and after this date the spectral structure of the same vegetation type changes systematically with seasonal phenology and interannual variability. Accordingly, the comparison of matching methods was aimed at evaluating the extent to which each approach preserves the correspondence between clusters of the same type over time under conditions of such drift. The sensitivity of matching methods to temporal drift was analysed separately in two aspects. The first aspect, represented within-year phenological variability, was defined as the distance between the satellite observation date and the reference survey date of the corresponding year. For each observation associated with survey year i., the reference centre was defined as:

$$d_i^{\mathrm{r}\mathrm{e}\mathrm{f}}=\mathrm{D}\mathrm{a}\mathrm{t}\mathrm{e}(y_i,\mathrm{08,15}),$$

that is, 15 August 2024 or 15 August 2025, depending on the value of $Year$. Geobotanical surveys were conducted in August 2024 and 2025, and the reference date was conventionally set to 15 August. For an observation acquired on date i, the signed and absolute temporal offsets were then calculated as:

$${\mathsf{\Delta }}_i=d_i-d_i^{\mathrm{r}\mathrm{e}\mathrm{f}},$$

$$\mid {\mathsf{\Delta }}_i\mid =\mid d_i-d_i^{\mathrm{r}\mathrm{e}\mathrm{f}}\mid ,$$

where $d_i$ is the date of the $i$-th satellite observation, $d_i^{\mathrm{r}\mathrm{e}\mathrm{f}}$is the reference date corresponding to 15 August of the calendar year to which the observation belongs, ${\mathsf{\Delta }}_i$is the signed temporal offset in days, and $\mid {\mathsf{\Delta }}_i\mid$ is its absolute magnitude. At the level of each “class × date” combination, these values were summarised using the mean and median.

The second aspect represented the interannual component of variability. It was defined as the absolute difference between the year of the satellite observation and the modal year of the reference geobotanical surveys forming the corresponding “class × date” combination. In addition, the number of represented years, as well as the minimum, maximum, and modal survey year, were determined to characterise the composition of the reference dataset. This formulation enabled separating within-year phenological drift from the interannual component of variability. In this way, interannual drift was not an abstract temporal distance but specifically the discrepancy between the image year and the dominant year of the reference material underlying the corresponding cluster. The phenological and interannual drift metrics were subsequently joined to the benchmark results for each date pair. Matching quality was evaluated using three indicators: correctness of the top-ranked correspondence, rank of the true class among all candidates, and the margin between the true class and the nearest alternative candidate.

For the correctness of the top-ranked correspondence, a generalised linear model (GLM) with a binomial distribution was used. In contrast, linear models (LM) were applied to the log-transformed rank of the true class and to the margin metric. All models included fixed effects of method, the corresponding temporal-drift component, and the interaction between method and drift (method × drift). A statistically significant interaction was interpreted as evidence that sensitivity to temporal drift differs among methods, i.e. that the rate of deterioration in matching quality depends on the selected approach. Thus, the analysis tested not only the presence of an overall effect of temporal distance, but also whether its strength varies across matching methods. This formulation directly corresponds to the testing strategy defined in the manuscript, in which the key criterion is the interaction $method\times temporaldistance$.

For phenological drift, separate models were evaluated in which the predictor was the absolute deviation of the image date from 15 August of the same year. For interannual drift, analogous models were fitted using the absolute difference between the image year and the modal year of the reference surveys as the predictor. Within both drift components, comparisons focused on how rapidly increasing temporal distance reduced the probability of correct correspondence, increased the rank of the true class, and altered the margin metric. Slower degradation of these indicators was interpreted as evidence of greater temporal stability of the corresponding method.

2.12. Unsupervised Mode: Temporal Matching Without Reference Labels

In the unsupervised mode, clustering was performed independently for each date; thus, the resulting clusters initially lacked predefined correspondence across temporal slices. To establish temporal continuity, a global cluster library was constructed sequentially through the time series. Each local cluster identified on the current date was compared with the existing elements of the global library and either absorbed into an existing global code or assigned a new code if no sufficiently similar correspondence was found. This approach avoided forcing newly emerged post-catastrophic surface states into previously existing classes. Temporal matching was based on a combined criterion integrating spectral similarity and spatial continuity. Spectral similarity evaluated how closely the current cluster corresponded to previously observed clusters in the common spectral space. In contrast, spatial continuity quantified the extent to which the current cluster occupied areas previously associated with a given global code. If the combined similarity score exceeded a predefined threshold, the cluster inherited the corresponding global code; otherwise, a new global code was created. Consequently, temporal matching represented not only changes in spectral properties but also the spatial memory of landscape structure. After construction of the temporal cluster library, the resulting global codes were transferred to the full raster domain. In the present study, this transfer was performed using a Random Forest classifier trained on the cluster library. In contrast, an earlier implementation of the procedure used Linear Discriminant Analysis (LDA) for the same purpose. The classifier was then applied to all pixels for each date, after which the temporal stability of codes, the emergence of new codes, and the structure of superclusters were analysed.

2.13. Identification of Superclusters Based on the Graph of Temporal Correspondences

A graph of temporal correspondences was constructed to aggregate the primary temporal codes into more stable higher-level units. The vertices of the graph represented the global cluster codes produced during the temporal matching procedure, whereas the edges represented recurrent correspondences between codes on consecutive dates. Edge weights characterised the strength of these temporal correspondences. To remove random or weak connections, a minimum edge-weight threshold i.i. ($min\_edge$) was selected in advance by scanning a range of threshold values and identifying the breakpoint beyond which further increases substantially altered the graph structure. After threshold filtering, the graph was partitioned into communities using the Walktrap algorithm [64]. The resulting communities were interpreted as superclusters, i.e. groups of temporal codes that were repeatedly linked throughout the temporal sequence and therefore represented more stable trajectory-level units of landscape cover. Based on the community membership of graph vertices, a correspondence table between the original temporal codes and superclusters was constructed and subsequently used to map and compare methods. To evaluate the structural stability of the superclustering procedure, several integral characteristics of the temporal correspondence graph were additionally analysed, including the number of vertices and edges, graph density, weighted density, graph modularity, the number of connected components, the proportion of the largest connected component, and the entropy of the supercluster distribution. Modularity characterised the strength of graph partitioning into communities, i.e. the degree to which superclusters were separated from one another. The parametric stability of the graph structure was additionally assessed by varying the minimum edge-weight threshold (min_edge) and quantifying its effects on the number of superclusters, graph modularity, and graph connectivity. This enabled the determination of whether the identified superclusters reflected a stable temporal organisation of cluster trajectories rather than artefacts of a particular threshold value.

2.14. Wall-to-Wall Prediction of Cluster Classes and Transfer Quality Assessment

After construction of the cluster library and the correspondence table linking the primary codes to superclusters, classification decisions were transferred from the reference-pixel sample to the full raster for each date. For this purpose, the saved cluster-library object was used, containing the predictor set, standardisation parameters, cluster codes, and the trained classifier. In the present study, the transfer was performed using Random Forest, which was applied to the complete set of pixels for each scene within the same feature space in which the cluster decisions had originally been generated. For each date, a raster map of the primary cluster codes was produced, followed by a supercluster map after application of the correspondence table. In this way, local cluster decisions obtained from the pixel sample were propagated to the entire study area while preserving the consistent temporal numbering of clusters and superclusters. Transfer quality was evaluated through internal validation using the same temporal slices and cluster codes used to construct the cluster library. For each date, the cluster codes generated during the library-construction stage were compared with the codes predicted by the classifier for the corresponding pixels. The overall reliability of transfer was characterised using the mean values and standard deviations of accuracy, balanced accuracy, and macro-F1 [65]. Accuracy represented the overall proportion of correctly reproduced codes, balanced accuracy accounted for unequal class representation, and macro-F1 characterised the balanced quality of recognition across all clusters irrespective of their spatial extent.

2.15. Integration of Temporal Maps Into a Spatio-Temporal Model

The series of supercluster maps was integrated into a unified spatio-temporal model to obtain a unified spatial model of landscape dynamics through time and to enable comparison of these dynamics among alternative matching procedures. After normalising these values, the relatively dominant supercluster was identified, and surfaces of relative membership were generated, reflecting the proximity of each pixel to the characteristic structure of its corresponding supercluster. At the next stage, contextual superclustering was performed. For each pixel, the local spatial context was evaluated by measuring distances to the nearest patches of different superclusters, and a profile of relative contextual proximity was constructed. Based on these profiles, second-order superclusters were identified, characterised not only by the temporal membership of pixels themselves but also by the similarity of their surrounding spatial context. Additional measures of classification uncertainty were also calculated, including the margin between the nearest alternatives and the entropy of fuzzy membership. To analyse the dynamics of second-order superclusters, temporal changes in their compositional structure were evaluated. For each supercluster, the proportions of the primary superclusters across different dates were determined, and the contributions of the long-term trend and seasonal components were analysed using partial RDA. Based on the relative importance of the trend and seasonal components, superclusters were classified into seasonal, directional-trend, or mixed-dynamic types.

2.16. Comparison of Supercluster Maps Derived from Different Matching Methods

Correspondence among supercluster solutions derived from different matching methods was evaluated through the spatial overlap of identically labelled or functionally equivalent superclusters. Because different methods could produce different numbers of superclusters and assign incompatible numerical codes to them, labels were first harmonised relative to the Random Forest reference map. For each supercluster produced by an alternative method, the closest supercluster in the reference map was identified using the maximum Jaccard index. The maps were then compared within this common label system. Based on the harmonised maps, overall agreement among methods was assessed using Fleiss’ kappa [66], pairwise agreement between methods, and the hierarchical structure of method similarity based on the distance $1-agreement$. In addition, consensus strength was calculated for each pixel as the proportion of methods assigning that pixel to the dominant harmonised supercluster. This metric was used as a spatial measure of interpretative reliability: high values indicated stable areas reproduced consistently by different methods, whereas low values indicated zones of methodological uncertainty or transitional landscape structure.

3. Results

3.1. Geometric Structure of Spectral Trajectories

At the level of individual habitat types, temporal changes exhibited an ordered trajectory-like structure within the common ordination space. However, the shape and direction of these trajectories differed substantially among classes (Figure 2). For most habitat types, successive states formed elongated point clouds with a consistent orientation of the principal axis, indicating a pronounced dominant direction of spectral variability.

3.2. Geometric Properties of Spectral Clusters

The geometric criterion of temporal correspondence showed that same-class transitions between neighbouring dates had statistically significantly smaller angles between the principal axes of spectral clusters than the nearest between-class alternatives(Wilcoxon test, $V=67496$, $p<2.2\times {10}^{-16}$). Similarly, the normalised transverse gap between spectral clusters was statistically significantly smaller for within-class pairs than for the nearest between-class pairs ($V=65210$, $p<2.2\times {10}^{-16}$). The combined geometric score integrating angular and spatial components was likewise statistically significantly smaller for within-class transitions than for the nearest between-class alternatives ($V=66782$, $p<2.2\times {10}^{-16}$). In addition, margin score values were statistically significantly greater than zero ($V=153334$, $p<2.2\times {10}^{-16}$), indicating the presence of a positive separation margin between the true within-class correspondence and the nearest alternative class. The distribution of points relative to the equality line showed that, in most cases, within-class temporal pairs had lower combined scores than the nearest between-class alternatives (Figure 3). This was reflected by the predominance of points above the diagonal. However, the distribution was not homogeneous, with a noticeable proportion of cases below the equality line. The proportion of successful matches according to this criterion was 72.6%, indicating a predominant, although not universal, ability of the geometric approach to correctly distinguish temporally adjacent states of the same class. Consistently, the median $margin\_score$value was positive (0.097), and the one-sample Wilcoxon test showed that $margin\_score$was statistically significantly greater than zero. The true within-class correspondence tended to have a lower combined score than the nearest between-class alternative.

At the level of the normalised transverse gap, within-class temporal pairs in most cases had lower gap norm values than the nearest between-class alternatives. Such cases accounted for 69.9% of all comparisons, indicating a predominant, although not universal, advantage of within-class correspondences. The median normalised transverse gap for within-class pairs was lower than that of the nearest between-class alternatives ($-0.824$ versus $-0.747$), and the paired Wilcoxon test confirmed a statistically significant shift towards smaller within-class values ($V=65210$, $p<2.2\times {10}^{-16}$).

Because negative gap norm values in the applied metric correspond to spatial overlap or contact between the axial “tubes”, these results indicate that temporally adjacent states of the same class were, on average, characterised by stronger geometric proximity than the nearest between-class alternatives. At the same time, the distribution of points was not homogeneous, with approximately 28.3% of cases failing to meet this criterion, indicating a substantial proportion of ambiguous transitions. Among within-class pairs, overlap of the axial regions was observed in 95.3% of cases, further supporting a high, although not absolute, degree of temporal connectivity of spectral clusters in the gap norm feature space.

3.3. Comparison of the Axis-Distance and Geometric Matching Methods with Alternative Approaches

The Random Forest method achieved the highest top-1 accuracy (0.72) in a benchmark evaluation of six approaches to cluster matching between adjacent dates. In contrast, both proposed approaches also demonstrated strong performance: the Geometric method reached a top-1 accuracy of 0.70, and the Axis-distance method achieved 0.68 (Table 2). At the same time, the highest top-3 accuracy was obtained by the Axis-distance method (0.93). In contrast, the Geometric method and Random Forest both showed identical values of 0.92, indicating stable inclusion of the true class among the nearest candidates. Both proposed methods substantially outperformed the Centroid and LDA methods, and especially the Mahalanobis approach. According to ranking-based metrics, the Axis-distance method also occupied one of the leading positions: the median rank of the true class was equal to 1, and the mean rank was 1.65, which was better than those of the Geometric method (1.68), Centroid (1.79), LDA (1.84), and Mahalanobis (2.44), and was surpassed only by Random Forest (1.59). Thus, Random Forest remained the best-performing approach under the strict top-1 accuracy criterion; however, both proposed geometrically oriented approaches demonstrated high competitiveness. The Geometric method achieved nearly maximal top-1 accuracy, whereas the Axis-distance method achieved the best top-3 accuracy and one of the strongest overall ranking performances.

Matching accuracy generally increased during the year for all tested methods, from the winter–spring period towards the second half of summer and autumn. However, the pattern of this increase differed substantially among approaches (Figure 4). The largest amplitude of seasonal change was observed for the axis-distance and geometric methods. For axis-distance, predicted top-1 accuracy ranged from 0.315 to 0.840, with an amplitude of 0.525. For the geometric method, the corresponding values were 0.348, 0.872, and 0.524. Random Forest also showed a pronounced seasonal increase, although slightly weaker than that of the two proposed geometrically oriented approaches: predicted accuracy ranged from 0.470 to 0.861, with an amplitude of 0.391. For LDA and centroid distance, the seasonal trajectory was smoother and closer to monotonic: predicted accuracy ranged from 0.346 to 0.727 for LDA and from 0.425 to 0.737 for centroid distance, with amplitudes of 0.380 and 0.311, respectively. The Mahalanobis method showed the lowest accuracy throughout the year, with a predicted range of only 0.265–0.524 and the smallest seasonal amplitude among all methods, 0.259.

GAM analysis confirmed a statistically significant seasonal trend for all methods. At the same time, the effective degrees of freedom ($edf$) indicated that the trend shape differed across approaches. For LDA ($edf=1.00$) and centroid distance ($edf=1.00$), the seasonal relationship was effectively close to linear. For the geometric method ($edf=2.23$) and Random Forest ($edf=2.41$), moderate nonlinearity was observed. In contrast, the seasonal trajectory was more complex for axis-distance ($edf=4.16$) and especially for Mahalanobis distance ($edf=2.89$). Thus, the proposed approaches, particularly the axis-distance and geometric methods, showed the most pronounced seasonal restructuring of accuracy, indicating high sensitivity to within-year changes in the spectral structure of vegetation cover. At the same time, Random Forest also achieved high accuracy in the second half of the season. In contrast, centroid distance, LDA, and especially Mahalanobis distance were inferior in both the level and amplitude of the seasonal maximum.

3.4. Matching Error Under Increasing Temporal Divergence from the Reference Surveys

For phenological drift, the global effect of temporal distance from the reference date was not statistically significant for the correctness of the top-ranked match. Still, it was statistically significant for the rank of the true class ($p<0.001$) and for the separation margin between the true class and the nearest alternative ($p<0.001$) (Table 3). The $method\times drift$interaction was statistically significant for top-ranked correctness ($p=0.019$) and separation margin ($p<0.001$), but was not significant for the rank of the true class ($p=0.72$). This indicates that sensitivity to phenological drift differed among methods primarily in the probability of correctly matching the top-ranked class and the degree of separation between the true class and the nearest alternative. In contrast, the rank of the true class showed greater deterioration across methods. For top-ranked correctness, centroid distance showed the highest stability ($-0.0052$), followed by Mahalanobis distance ($-0.0059$) and LDA ($-0.0067$), whereas Random Forest ($-0.0094$), axis-distance ($-0.0103$), and especially the geometric method ($-0.0117$) exhibited a faster decline in correctness with increasing phenological drift. For the rank of the true class, the slowest rate of deterioration was observed for LDA ($0.0014$), followed by centroid distance ($0.0016$) and Random Forest ($0.0016$), whereas axis-distance ($0.0018$) and Mahalanobis distance ($0.0018$) were slightly less stable. The geometric method ($0.0020$) showed the fastest increase in the true class’s rank. For the separation margin, the highest stability was observed for axis-distance ($-0.0005$); very similar results were obtained for LDA ($-0.0008$), centroid distance ($-0.0009$), and Random Forest ($-0.0009$), whereas the geometric method ($-0.0020$) lost separation margin more rapidly. Mahalanobis distance ($-0.0132$) differed strongly in the direction of the lowest stability. Thus, under phenological drift, axis-distance did not show an advantage in terms of top-ranked correctness or rank of the true class, but proved to be the most stable according to the separation-margin criterion, i.e. it best preserved the separation between the true class and the nearest alternative. In contrast, the geometric method showed lower phenological stability across all three metrics.

For interannual drift, the global effect of temporal distance was not statistically significant for the correctness of the top-ranked match. Still, it was clearly expressed for the rank of the true class ($p<0.001$) and for the separation margin between the true class and the nearest alternative ($p<0.001$) (Table 4). The $method\times drift$interaction was likewise not significant for the correctness of the top-ranked match ($p=0.63$). Still, it was statistically significant for the rank of the true class ($p=0.04$) and especially for the separation margin ($p<0.001$). This indicates that increasing interannual temporal distance was not associated with a consistent decline in the probability of correct top-ranked matching across all methods, but was associated with a deterioration in the rank position of the true class and a reduction in its separation from the nearest alternative. At the same time, the significant interaction effects for the rank of the true class and the separation margin indicate that sensitivity to interannual drift differed among methods, i.e. that the rate of deterioration in matching quality depended on the selected approach. For the rank of the true class, the slowest rate of deterioration was observed for centroid distance (0.01), followed by LDA (0.03), Random Forest (0.04), and axis-distance (0.04), whereas the rank of the true class increased more rapidly for the geometric method (0.05) and especially for Mahalanobis distance (0.07).

Note: Slopes represent the rate of change in matching quality with increasing interannual drift. More negative slopes for top-ranked correctness and separation margin indicate faster deterioration, whereas larger positive slopes for the rank of the true class indicate lower temporal stability. Stability ranks are ordered from most stable (1) to least stable (6) for each metric.

For the separation margin, the axis-distance method showed the highest stability (minus 0.02). The very similar results were obtained with Random Forest (minus 0.03) and LDA (minus$-0.03$). In contrast, the loss of separation margin was noticeably greater for centroid distance ($-0.04$), even more pronounced for the geometric method ($-0.10$), and extremely strong for Mahalanobis distance ($-1.18$). The axis-distance method combined a relatively moderate deterioration in top-ranked correctness with the highest stability in separation margin, indicating better preservation of the separation between the true class and the nearest alternative even under conditions of deterioration in the true class’s rank position.

3.5. Wall-to-Wall Prediction Accuracy for Different Matching Approaches

Model transfer over the full raster shows similar values of accuracy and balanced accuracy across most approaches, while differences between methods are more clearly expressed in macro-F1 (Table 5). The axis-distance method achieves the highest accuracy and macro-F1 among the tested approaches. In contrast, the Mahalanobis method exhibits lower overall accuracy but higher macro-F1, indicating improved class balance at the expense of overall correctness. The centroid, LDA, and Random Forest approaches show comparable performance, with no pronounced advantages.

3.6. Internal Stability of Scene-Matching Methods

The internal stability analysis revealed substantial differences among the scene-matching approaches (Table 6). Random Forest and LDA produced the fewest new codes, indicating the most stable temporal matching structure. Axis-distance also showed relatively high stability, with only 32 new codes across 93 dates. In contrast, the Mahalanobis approach generated a markedly higher number of new codes and a much higher mean number per date, indicating strong temporal fragmentation and lower internal stability. The geometric and centroid approaches showed intermediate and very similar behaviour. The results suggest that Mahalanobis distance is the least stable matching strategy, whereas Random Forest, LDA, and Axis-distance provide more temporally consistent code propagation.

3.7. Graph Structure and Robustness of Supercluster Aggregation Among Matching Approaches

The supercluster graphs were constructed from recurrent transitions between codes across consecutive dates, reflecting their mutual temporal correspondence (Figure 5). Using a common thresholding rule, these transitions were integrated into a unified matching structure, from which the final supercluster grouping was derived. The resulting graphs represent the aggregated matching solutions, characterised by specific numbers and configurations of superclusters for each method.

The comparative analysis reveals clear differentiation among the evaluated methods in terms of graph structure and parametric robustness, with consistent patterns observed across metrics and their integrated representation (Table 7). The constructed graphs differ markedly in their size, connectivity, and degree of structural organisation, reflecting differences in how temporal correspondences between clusters are aggregated. At the structural level, methods such as Geometric and Centroid produce relatively dense and well-connected graphs, characterised by a small number of components and a high share of nodes concentrated within the largest component. In contrast, Axis-distance, LDA, and Random Forest yield more fragmented structures, with multiple disconnected components and lower largest-component dominance, indicating a more conservative matching behaviour.

Differences are also evident in modular organisation. LDA exhibits the highest modularity, suggesting clearer separation between superclusters, whereas Centroid shows lower modularity, reflecting less distinct community structure. The number and entropy of superclusters further indicate variation in the distribution of cluster relationships across the graph. The most pronounced contrast is observed for the Mahalanobis approach, which demonstrates extreme behaviour across nearly all metrics. It produces a substantially larger graph with a very high number of nodes and edges, resulting in a fully connected structure (single component) with maximal largest-component share. Despite relatively high modularity, this configuration is associated with excessive connectivity and a lack of structural sparsity, resulting in markedly elevated fragility indices and substantial parameter-dependent variability. These results indicate that methods differ in matching accuracy and in the structural properties of the resulting temporal correspondence networks. Approaches such as Axis-distance and Random Forest maintain more fragmented and stable graph structures, whereas Mahalanobis introduces excessive connectivity and reduced robustness under parameter variation.

The ordination provides a more coherent interpretation of these relationships by projecting the multidimensional metric space into a continuous structure, where distances among methods reflect their overall dissimilarity (Figure 6).

The configuration shows a distinct separation of Mahalanobis along the gradient defined by graph fragility and range-based instability, confirming its deviation from the other methods. In contrast, LDA and Random Forest are aligned with higher modularity and component dominance, indicating a tendency toward more consolidated but less fragmented structures. Axis-distance occupies an intermediate position, without strong alignment with any single metric vector, indicating a balanced trade-off between structural and stability-related properties. This is consistent with its minimal distance to the ideal solution, suggesting that it does not maximise individual metrics but avoids extreme values across the entire set. Centroid and Geometric methods are located near the density-related vectors, reflecting their association with increased edge density and weighted connectivity, while remaining separated from the instability gradients. The orientation of vectors further clarifies the underlying structure of the metric space. Variables related to graph size and connectivity (nodes, edges, component range) form a coherent gradient opposing modularity-related measures, while fragility and variability indicators define an independent axis of instability. The relative positions of the methods along these gradients demonstrate that differences among them arise not from a single dominant property, but from distinct combinations of structural compactness, modular organisation, and robustness.

3.8. Ordination of Graph Structure and Robustness Profiles of Matching Methods

The overall agreement among the six alternative supercluster maps was moderate and statistically significant (Fleiss’ κ = 0.363, z = 1332, p < 0.001). This indicates that the methods reproduced a common spatial structure, although certain discrepancies remained in the assignment of pixels to superclusters. Pairwise comparisons between methods were further visualised using hierarchical clustering based on the distance measure 1—agreement (Figure 7a). This representation allowed the grouping structure among the alternative classification approaches to be assessed. The resulting dendrogram showed that the two geometric methods, spectral ellipse distance and axis distance, formed the closest pair, whereas Random Forest was the most distinct from the remaining methods.

Different scene-matching methods produced alternative supercluster solutions, and the relationships among these solutions were quantified using consensus strength. This index describes the proportion of methods that assigned the dominant supercluster label to each pixel. Higher values, therefore, indicate stronger agreement among methods, whereas lower values indicate spatially localised uncertainty in supercluster assignment. The mean consensus strength was 0.612 ± 0.225, with a median of 0.500. The interquartile range was 0.500-0.833, indicating spatial heterogeneity in agreement among the alternative matching methods. Spatial variation in consensus strength was clearly structured by the main geomorphological gradients of the study area (Figure 7b). The highest values were associated with the highest-relief positions and open water surfaces, indicating that these spatial units were consistently identified across different matching methods. In contrast, the dynamically changing transitional part of the relief showed substantially lower consensus strength, reflecting bigger method-dependent differences in supercluster assignment. This pattern suggests that disagreement among methods was concentrated mainly in ecotonal and geomorphologically transitional zones, where landscape-cover conditions are more heterogeneous and temporally unstable.

3.9. Spatio-Temporal Consistency of Aligned Supercluster Solutions Across Matching Methods

The aligned supercluster solutions showed a generally consistent spatial structure across scene-matching methods and reproduced similar temporal patterns. Although the relative areas and model fits varied across methods for individual superclusters, the dominant dynamic types were broadly comparable. Directional trends were repeatedly identified for several superclusters, whereas seasonal or mixed patterns occurred mainly in more transitional units (Table 1). This indicates that the alternative matching methods differed in the precise spatial allocation of some pixels but converged in identifying the major spatiotemporal organisation of the landscape.

The aligned supercluster maps represent an integrated spatial synthesis of the alternative scene-matching approaches (Figure 8).

Table 8. Dynamic characteristics of aligned superclusters across scene-matching methods (all reported R² values are statistically significant).

Statistic*	Method**	Supercluster
		1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
Area (%)	AD	11.60	15.00	5.30	8.50	14.30	6.90	5.90	6.40	–	9.30	3.70	8.40	4.70	–	–
	GE	9.50	11.00		11.00	11.50	11.00	4.60	4.90	–	8.40	8.20	–	7.10	7.50	5.20
	LDA	6.60	6.20	12.70	6.40	7.00	8.10	7.20	–	–	9.50	5.40	6.60	12.10	12.30	–
	CEN	11.86	13.08	3.16	10.78	8.90	7.40	4.67	7.10	–	8.88	–	7.48	6.50	10.17	–
	MAH	11.30	10.90	4.30	–	3.60	13.30	9.80	9.70	8.20	8.60	–	–	7.10	4.90	8.30
	RF	10.90	15.50	7.30	9.00	13.90	6.50	5.80	5.40	8.00	7.40	4.70	5.60	–	–	–
Trend R²	AD	0.49	0.27	0.40	0.26	0.27	0.15	0.14	0.43	–	0.41	0.32	0.32	0.21	–	–
	GE	0.31	0.34		0.31	0.33	0.20	0.03	0.26	–	0.37	0.28	–	0.41	0.30	0.44
	LDA	0.14	0.55	0.26	0.09	0.32	0.15	0.07	–	–	0.27	0.50	0.39	0.24	0.32	–
	CEN	0.30	0.19	0.43	0.19	0.34	0.14	0.07	0.43	–	0.25	–	0.32	0.19	0.27	–
	MAH	0.26	0.20	0.30		0.13	0.26	0.06	0.18	0.16	0.29	–	–	0.24	0.14	0.20
	RF	0.38	0.21	0.37	0.29	0.34	0.14	0.01	0.29	0.36	0.40	0.30	0.34	–	–	–
Season R²	AD	0.05	0.06	0.12	0.21	0.26	0.25	0.13	0.12	–	0.13	0.17	0.19	0.19	–	–
	GE	0.14	0.02	–	0.10	0.05	0.11	0.14	0.16	–	0.15	0.14	–	0.12	0.07	0.05
	LDA	0.21	0.09	0.19	0.14	0.09	0.19	0.23	–	–	0.21	0.11	0.17	0.21	0.00	–
	CEN	0.27	0.23	0.23	0.37	0.26	0.31	0.27	0.13	–	0.33	–	0.25	0.42	0.32	–
	MAH	0.24	0.13	0.21	–	0.34	0.21	0.10	0.31	0.22	0.22	–	–	0.14	0.35	0.27
	RF	0.11	0.04	0.20	0.19	0.14	0.25	0.22	0.17	0.17	0.17	0.20	0.15	–	–	–
Type***	AD	DT	DT	DT	MIX	MIX	SEA	MIX	DT	–	DT	DT	DT	MIX	–	–
	GE	DT	DT	–	DT	DT	DT	SEA	DT	–	DT	DT	–	DT	DT	DT
	LDA	SEA	DT	DT	SEA	DT	SEA	SEA	–	–	DT	DT	DT	MIX	DT	–
	CEN	MIX	MIX	DT	SEA	DT	SEA	SEA	DT	–	SEA	–	DT	SEA	MIX	–
	MAH	MIX	DT	DT		SEA	MIX	SEA	SEA	SEA	DT	–	–	DT	SEA	SEA
	RF	DT	DT	DT	DT	DT	SEA	SEA	DT	DT	DT	DT	DT	–	–	–

Note: * – statistic denotes the type of supercluster-level characteristic reported in each row: Area (%) is the relative spatial extent of the supercluster; Trend R² is the adjusted R² of the directional temporal trend model; Season R² is the adjusted R² of the seasonal model; ** – AD – Axis-distance; GE – Geometric (spectral ellipse distance); LDA – Linear Discriminant Analysis; CEN – Centroid (Euclidean); MAH – Mahalanobis; RF – Random Forest.; *** – Type indicates the dominant dynamic pattern: DT – directional trend; SEA – seasonal; MIX – mixed.

Table 8. Dynamic characteristics of aligned superclusters across scene-matching methods (all reported R² values are statistically significant).

Statistic*	Method**	Supercluster
		1	2	3	4	5	6	7	8	9	10	11	12	13	14	15
Area (%)	AD	11.60	15.00	5.30	8.50	14.30	6.90	5.90	6.40	–	9.30	3.70	8.40	4.70	–	–
	GE	9.50	11.00		11.00	11.50	11.00	4.60	4.90	–	8.40	8.20	–	7.10	7.50	5.20
	LDA	6.60	6.20	12.70	6.40	7.00	8.10	7.20	–	–	9.50	5.40	6.60	12.10	12.30	–
	CEN	11.86	13.08	3.16	10.78	8.90	7.40	4.67	7.10	–	8.88	–	7.48	6.50	10.17	–
	MAH	11.30	10.90	4.30	–	3.60	13.30	9.80	9.70	8.20	8.60	–	–	7.10	4.90	8.30
	RF	10.90	15.50	7.30	9.00	13.90	6.50	5.80	5.40	8.00	7.40	4.70	5.60	–	–	–
Trend R²	AD	0.49	0.27	0.40	0.26	0.27	0.15	0.14	0.43	–	0.41	0.32	0.32	0.21	–	–
	GE	0.31	0.34		0.31	0.33	0.20	0.03	0.26	–	0.37	0.28	–	0.41	0.30	0.44
	LDA	0.14	0.55	0.26	0.09	0.32	0.15	0.07	–	–	0.27	0.50	0.39	0.24	0.32	–
	CEN	0.30	0.19	0.43	0.19	0.34	0.14	0.07	0.43	–	0.25	–	0.32	0.19	0.27	–
	MAH	0.26	0.20	0.30		0.13	0.26	0.06	0.18	0.16	0.29	–	–	0.24	0.14	0.20
	RF	0.38	0.21	0.37	0.29	0.34	0.14	0.01	0.29	0.36	0.40	0.30	0.34	–	–	–
Season R²	AD	0.05	0.06	0.12	0.21	0.26	0.25	0.13	0.12	–	0.13	0.17	0.19	0.19	–	–
	GE	0.14	0.02	–	0.10	0.05	0.11	0.14	0.16	–	0.15	0.14	–	0.12	0.07	0.05
	LDA	0.21	0.09	0.19	0.14	0.09	0.19	0.23	–	–	0.21	0.11	0.17	0.21	0.00	–
	CEN	0.27	0.23	0.23	0.37	0.26	0.31	0.27	0.13	–	0.33	–	0.25	0.42	0.32	–
	MAH	0.24	0.13	0.21	–	0.34	0.21	0.10	0.31	0.22	0.22	–	–	0.14	0.35	0.27
	RF	0.11	0.04	0.20	0.19	0.14	0.25	0.22	0.17	0.17	0.17	0.20	0.15	–	–	–
Type***	AD	DT	DT	DT	MIX	MIX	SEA	MIX	DT	–	DT	DT	DT	MIX	–	–
	GE	DT	DT	–	DT	DT	DT	SEA	DT	–	DT	DT	–	DT	DT	DT
	LDA	SEA	DT	DT	SEA	DT	SEA	SEA	–	–	DT	DT	DT	MIX	DT	–
	CEN	MIX	MIX	DT	SEA	DT	SEA	SEA	DT	–	SEA	–	DT	SEA	MIX	–
	MAH	MIX	DT	DT		SEA	MIX	SEA	SEA	SEA	DT	–	–	DT	SEA	SEA
	RF	DT	DT	DT	DT	DT	SEA	SEA	DT	DT	DT	DT	DT	–	–	–

Note: * – statistic denotes the type of supercluster-level characteristic reported in each row: Area (%) is the relative spatial extent of the supercluster; Trend R² is the adjusted R² of the directional temporal trend model; Season R² is the adjusted R² of the seasonal model; ** – AD – Axis-distance; GE – Geometric (spectral ellipse distance); LDA – Linear Discriminant Analysis; CEN – Centroid (Euclidean); MAH – Mahalanobis; RF – Random Forest.; *** – Type indicates the dominant dynamic pattern: DT – directional trend; SEA – seasonal; MIX – mixed.

Despite local differences among methods, the resulting spatial patterns remained broadly consistent. They highlighted the same major geomorphological and temporal structures of the landscape. These aligned solutions therefore provide a common analytical framework that can serve as a basis for subsequent quantitative assessments using landscape-ecological approaches, including analyses of spatial heterogeneity, connectivity, ecotonal structure, and temporal landscape dynamics.

4. Discussion

4.1. Temporal Matching as a Problem of Continuity in Dynamic Landscapes

Temporal matching in satellite image analysis is a complex methodological problem because spectral correspondence between images depends not only on acquisition conditions but also on the temporal variability of the observed landscape system [16,41,67]. The problem of spectral matching of satellite imagery encompasses two interrelated dimensions: observational and conceptual. The observational dimension relates to the temporal gap between acquisition dates, during which spectral continuity between images progressively deteriorates [68]. As the interval between observations increases, differences caused by phenological development, moisture conditions, illumination, and other acquisition-related factors may accumulate, altering not only individual spectral values but also the overall configuration of spectral relationships among classes [69,70,71]. Hypothetically, if the temporal resolution of the image series were sufficiently high to render this degradation negligible, different matching approaches would be expected to yield similarly high matching accuracy. In terms of its observational design, the Sentinel-2 constellation is intended to provide a revisit frequency sufficient for accurate and reliable temporal matching, because substantial phenological changes in vegetation cover are generally unlikely to occur within such short intervals. However, meteorological conditions significantly constrain the availability of high-quality, cloud-free imagery, thereby introducing considerable irregularity into the dataset’s effective temporal structure [72,73]. Our image series for 2021–early 2026 confirms this limitation: despite the nominal revisit capacity of Sentinel-2, the effective multi-temporal dataset comprised 93 usable observation dates, with a mean interval of 20.1 days between consecutive images and a median interval of 10 days. The minimum interval was 2 days, whereas the maximum interval reached 100 days. The shortest intervals between high-quality images were most typical in the second half of summer and early autumn, whereas the longest gaps occurred during winter and early spring. This pattern is consistent with the general seasonal availability of usable optical satellite observations in temperate Europe, where cloud-, shadow-, and snow-free observations are unevenly distributed through the year and often require compositing or interpolation to obtain uninterrupted time series [74,75]. Importantly, the most pronounced phenological transitions occur during spring, precisely when the availability of suitable imagery is often at its lowest. This makes temporal matching particularly challenging and methodologically relevant under conditions where the continuity of satellite image time series is substantially interrupted.

Gaps between observation dates can be partly compensated for by using observations from other years acquired at similar days of the year, thereby reconstructing an average or expected annual phenological trajectory from multi-year satellite time series [40,74,76]. The calendar year is therefore ignored as an independent temporal dimension. In contrast, the day of the year is used as the primary reference for aligning observations along a common seasonal trajectory [77]. However, under highly dynamic or post-disturbance conditions, interannual variability becomes ecologically significant in its own right. In such systems, combining observations from different years into a single phenological trajectory may obscure abrupt transitions, accelerated succession, and regime shifts. Therefore, unsupervised classification is particularly important for catastrophe-induced dynamic landscapes, as it allows newly emerging surface types to be incorporated without forcing them into predefined classes. Nevertheless, this flexibility does not eliminate the problem of temporal gaps: when the continuity of satellite image series is interrupted, spectral states must still be linked, separated, or registered as new across unevenly spaced observation dates.

A key conceptual issue concerns the landscape unit’s category. Supervised classification implicitly treats landscape units as a priori-defined entities that objectively exist as discrete, identifiable classes [78,79,80]. Although this assumption simplifies reality, it is methodologically convenient and practically useful. However, it also imposes an important limitation: the set of possible system states is effectively treated as fixed. Even if this set can hypothetically be expanded, the spectrum of admissible states remains constrained by prior knowledge and predefined categories. As a result, such an approach leaves little room for genuinely new system states to emerge. For modelling processes in dynamically changing systems, it is essential to retain the ability to identify novel states and incorporate them into the system’s temporal trajectory. In this context, the concept of landscape continuity becomes particularly important [81]. How can relationships between different system states be established through time when the states themselves are changing? One possibility is system stationarity, in which the system largely preserves its properties over time. Another is directional transformation, where change proceeds along one dominant trajectory or several alternative trajectories. A further possibility is the disappearance of one system and the emergence of another in its place. Under such conditions, a fundamental problem arises: how can one determine whether the original system has merely changed its properties or whether an entirely new system has emerged? This conceptual issue gives rise to the natural spatial heterogeneity of the entity that we identify as a landscape-cover type. In our view, the concept of landscape continuity explains why natural spatial heterogeneity contains the key to resolving the problem of temporal matching.

The two complementary criteria were therefore employed for the spatio-temporal matching of evolving landscape-cover states: spatial position, which reflects spatial continuity and the inheritance of landscape structure, and spectral similarity, which indicates correspondence in surface properties. Natural processes result from long-term successional dynamics of landscape cover, in which some system states are successively replaced by others, forming ordered sequences of transformations [82]. Catastrophic events of natural or anthropogenic origin may abruptly interrupt such continuity [83]. However, they themselves initiate new chains of interconnected processes, thereby spatially inheriting their effects. For this reason, the spatial criterion is critically important for matching procedures. Spectral similarity reflects the qualitative identity of a particular cover type and its distinction from other types. At the same time, spectral properties may change while the system still preserves its qualitative identity. A territory characterised by a specific set of spectral properties, therefore, represents a particular type of landscape cover. Temporal matching becomes complicated because both the spatial extent of such cover types and their spectral characteristics are subject to spatial and temporal variability.

4.2. Geometric Interpretation of Spectral Heterogeneity

The results support the interpretation that the spectral representation of a given cover type forms a structured manifold of admissible states rather than a compact and internally homogeneous entity. From the perspective of landscape continuity, asynchronous phenological and successional dynamics within the same cover type produce ordered internal variability, which is expressed geometrically as cluster elongation. Differences in moisture conditions, microrelief, species composition, and disturbance history produce continuous variation among patches. Therefore, the spectral representation of a landscape-cover type forms an elongated cloud rather than a compact spheroidal one. Under this interpretation, the internal geometry of the spectral point cloud may contain information about the temporal trend of transformation within a landscape-cover type. Our findings support this assumption: the same landscape-cover types form inherited trajectories in spectral space. Even when their typical spectral values change substantially through time, temporally adjacent states tend to remain connected within the same structured manifold rather than shifting to unrelated regions of spectral space.

The principal axis of the spectral point cloud can be regarded as a local trace of the broader trajectory followed by a landscape-cover type in spectral space. It is important to distinguish between the global structure of the spectral space and the local geometry of individual spectral point clouds. The ordination space derived from PCA/rPCA characterises the dominant structure of spectral variability across the entire dataset, whereas the principal axes used in the matching procedure are estimated separately for each landscape-cover type and describe the internal geometry of individual spectral clouds within that shared coordinate system. Thus, the directional properties used for matching are not imposed directly by the ordination itself, but emerge from the organisation of independently formed clusters in the common spectral space. It therefore provides a basis for reconstructing the main direction of landscape-cover transformation across successive temporal states. Similar to ecological trajectories in ordination analysis, it represents the main gradient of within-cluster variability and therefore approximates the most probable direction of subsequent temporal displacement. This becomes particularly important during phenological and post-disturbance changes. The results of this study indicate that the identity of a landscape-cover type is not always preserved by centroid position alone: even when the centroid shifts substantially, the orientation of the spectral point cloud may retain information about the inherited trajectory of the same type in spectral space. These findings show that temporally adjacent states of the same landscape-cover type may shift substantially in spectral space while their spectral point clouds remain aligned along the same directional trajectory. The preservation of within-class directional similarity indicates that the observed spectral trajectories were not arbitrary but were constrained by the internal structure of the corresponding landscape-cover type.

For clusters belonging to the same landscape-cover type, overlap of axial regions between consecutive dates provides the geometric basis for reliable axis-based matching, because temporal correspondence is preserved through inherited alignment along the dominant axis of spectral variability. Ecological systems may undergo substantial spectral shifts without losing their qualitative identity, particularly during seasonal or post-disturbance dynamics [41,84]. Under such conditions, centroid-based approaches may interpret directional displacement within the same landscape-cover type as class replacement, whereas geometric overlap preserves information about trajectory inheritance within that cover type. These findings indicate that the predominance of overlap among within-class transitions supports the hypothesis that spectral trajectories retain geometric inheritance across successive states. Conversely, incomplete overlap in some transitions may mark episodes of accelerated restructuring or the emergence of transitional landscape states under post-catastrophic conditions.

The findings of this study provide empirical support for the hypothesis of geometric connectivity among spectral clusters, although this connectivity was not universal across all transitions. Dedicated evaluation of the geometric criterion demonstrated that within-class transitions between neighbouring dates were characterised by smaller angular and spatial differences than the nearest between-class alternatives, together with a positive separation margin between the true correspondence and the nearest competing class. This confirms that temporally adjacent states of the same class indeed preserve geometric similarity within the common ordination space, and therefore that the shape of the spectral cluster and the orientation of its principal axis can be used as an informative basis for temporal cluster matching. In a benchmark comparison of six approaches to cluster matching across neighbouring dates, Random Forest achieved the highest top-1 accuracy, although both proposed geometrically oriented approaches also performed well. The geometric method reached near-maximum top-1 accuracy and generally confirmed high practical suitability for matching neighbouring temporal states. The axis-distance method was not the best according to the strict top-1 accuracy criterion. Still, it demonstrated the highest top-3 accuracy and one of the best rank-based performances, indicating a stable ability to retain the true class among the nearest candidates. Thus, in practical terms, the proposed methods proved competitive with standard alternatives, although their strengths differed: the geometric method performed better under the strict first-choice criterion, whereas the axis-distance method showed greater rank stability and a broader margin of correct separation. Drift analysis showed that a single dimension cannot describe the temporal stability of the proposed approaches. For phenological drift, the global effect of temporal distance was significant for the rank of the true class and the separation margin. In contrast, the $method\times drift$ interaction was significant for top-ranked correctness and separation margin. Under these conditions, the axis-distance method did not show an advantage in terms of top-ranked correctness or rank of the true class. Still, it proved to be the most stable according to the separation-margin criterion, i.e. it best preserved separation between the true class and the nearest alternative. In contrast, the geometric method showed lower phenological stability and was inferior to more stable alternatives across all three metrics. For interannual drift, the global effect of temporal distance was likewise pronounced for the rank of the true class and the separation margin. At the same time, the $method\times drift$ interaction was significant specifically for these two metrics. In this case, the axis-distance method again demonstrated the highest stability according to the separation-margin criterion. However, for the rank of the true class, it was inferior to centroid distance, LDA, and Random Forest. The geometric method again showed no advantage and, under interannual variability, deteriorated more rapidly than axis-distance, centroid distance, LDA, and Random Forest. Thus, the axis-distance method should be interpreted as a structurally oriented and interpretable approach whose main strength lies not in maximising top-ranked correctness under all conditions, but in preserving the geometric identity of the true class and maintaining its separation from competing alternatives under temporal drift. The geometric method, in turn, is best regarded as an effective benchmark-level approach for short-term matching between neighbouring dates but offers little advantage in terms of temporal stability. Overall, both proposed approaches confirmed the value of geometric descriptions of spectral clusters for temporal matching tasks. However, among interpretable approaches, the axis-distance method proved the most promising when preservation of separation between the true class and the nearest alternatives is of primary importance.

4.3. Why Do the Geometric Method and Axis-Distance Behave Differently?

The comparison of the proposed approaches indicates that the geometric and axis-distance methods should not be interpreted simply as competing alternatives in which one universally outperforms the other. Instead, they represent two different matching regimes based on different assumptions about temporal continuity in spectral space. The geometric method demonstrated high matching accuracy for neighbouring dates, indicating that detailed cluster geometry contains substantial information about short-term temporal correspondence. Because the method simultaneously considers cluster orientation, transverse overlap, and local geometric configuration, it is highly sensitive to the fine structure of spectral clouds. This allows precise identification of temporally adjacent states when phenological displacement remains relatively limited. It should be emphasized that the dependence on detailed local geometry makes the method more sensitive to temporal drift. As phenological or interannual divergence increases, the spectral cloud’s shape changes, reducing the stability of local geometric correspondence. The geometric approach, therefore, behaves as a high-resolution local matching method that performs best when temporal continuity is preserved over short intervals. In contrast, the axis-distance method relies less on the exact shape of the spectral cloud and instead emphasises the dominant direction of temporal variability. Although its strict top-ranked accuracy was slightly lower, the method more consistently preserved separation between the true class and competing alternatives under both phenological and interannual drift. This suggests that the principal axis captures the general direction of admissible temporal displacement, even as the cloud’s detailed geometry changes over time. The method, therefore, acts as a low-dimensional temporal regulariser that constrains matching along coherent spectral trajectories rather than local cluster configurations. In this sense, axis-distance operates primarily as a trajectory continuity constraint, preserving the structural identity of landscape-cover types despite progressive temporal transformation.

4.4. Temporal Drift and Ecological Meaning of Matching Instability

The results demonstrate that temporal matching is strongly influenced by seasonal restructuring of spectral space. During the vegetation season, the relative positions and separability of habitat states change continuously due to phenological development, moisture fluctuations, and variation in vegetation productivity. This effect is especially pronounced during spring, when rapid transitions between phenological phases reduce the stability of spectral relationships and increase overlap between vegetation states. Consequently, the separability of landscape-cover types should be regarded as seasonally dynamic rather than constant. Certain landscape-cover types become more spectrally distinct during periods of maximal structural differentiation, whereas at other stages their trajectories partially converge in spectral space. The pronounced seasonal variation in matching accuracy observed for the geometric and axis-distance approaches indicates that these methods are sensitive to seasonal changes in the geometry of spectral trajectories formed by landscape-cover types. Because these methods explicitly use the geometric shape and orientation of spectral point clouds, they capture seasonal reorganisation in the spectral representation of landscape-cover types rather than relying solely on static boundaries between predefined cover categories.

In contrast to within-year phenological drift, interannual drift was used to evaluate how different matching methods respond to longer-term divergence between observation years and reference years. In the present dataset, this divergence occurred against the background of the 2023 catastrophe and subsequent landscape restructuring, but the main analytical focus was the differential stability of the matching approaches. The interannual-drift analysis therefore provided a diagnostic basis for distinguishing matching approaches according to their robustness: some methods were more prone to losing separation between the true correspondence and competing alternatives, whereas others preserved this separation more effectively. In this respect, axis-distance showed the highest stability according to the separation-margin criterion, indicating that it better maintained the distinction between the corresponding landscape-cover type and neighbouring alternatives under interannual divergence.

4.5. Structural Properties of Temporal Correspondence Graphs

The obtained graph structures show that each temporal matching method produces a distinct mode of temporal organisation of landscape dynamics, expressed jointly through classification accuracy, graph connectivity, fragmentation, and modularity. The graphs integrate recurrent correspondences between clusters through time and their topology reflects how each method reconstructs continuity, fragmentation, and transition within the evolving landscape system. More fragmented graph structures indicate a conservative matching strategy that accepts temporal continuity only for strongly supported correspondences. Such approaches tend to preserve sharper separation between trajectories and reduce the probability of merging ecologically distinct states. In contrast, denser and more connected graphs reflect more permissive matching behaviour, where broader ranges of spectral transitions are interpreted as belonging to the same temporal structure. Excessive graph connectivity may indicate loss of discriminative organisation within the temporal correspondence network. When too many transitions become connected, the graph increasingly reflects generalised spectral similarity rather than structured ecological trajectories. Under such conditions, temporal matching becomes less selective, and distinctions between alternative developmental pathways begin to collapse. The strongly connected structure produced by the Mahalanobis approach appears to illustrate this effect, where high connectivity is associated with reduced structural sparsity and lower robustness of the resulting temporal organisation. These results suggest that graph topology itself contains ecologically meaningful information about the reconstruction of landscape continuity through time. Fragmentation, modularity, and connectivity should therefore be interpreted not merely as technical properties of the matching procedure, but also as indicators of how different methods conceptualise temporal relationships between evolving landscape states.

4.6. Why Does Mahalanobis Become Unstable?

The instability of the Mahalanobis approach stems from its strong dependence on the covariance structure. This method assumes that the dispersion pattern of a reference class provides a reliable description of the expected variability of the same class at subsequent dates. However, in post-catastrophic landscapes, this assumption is weak. Hydrological restructuring, rapid colonisation of exposed substrates, and accelerated successional shifts alter not only the centroid position of spectral clusters but also their internal dispersion, elongation, and orientation. As a result, covariance-sensitive matching may begin to overconnect states that are statistically compatible in a broad distributional sense but ecologically belong to different trajectories. When the covariance structure becomes unstable, the Mahalanobis distance loses selectivity: it may treat divergent or transitional states as acceptable variants of the same class. This explains the excessive graph connectivity and high fragmentation of temporal codes observed for this method. In such systems, instability of covariance is not only a statistical problem but also an ecological signal that the landscape is undergoing structural reorganisation.

4.7. Axis-Distance as a Balanced Graph Solution

In contrast to both highly fragmented and excessively connected matching structures, the axis-distance approach produced a comparatively balanced graph organisation. Its topology combined moderate fragmentation with a relatively stable modular structure, indicating that the method preserved separation among temporal trajectories without collapsing them into a single, densely connected network. It avoided the excessive proliferation of disconnected components characteristic of overly conservative matching. This balance is reflected in the relatively low graph fragility and the proximity of the method to the multidimensional ideal solution. Ecologically, such behaviour suggests that the axis-distance approach preserves continuity of landscape trajectories while remaining sufficiently selective to maintain distinctions between alternative developmental pathways. Because the method primarily constrains matching along dominant temporal variability rather than detailed local geometry, it reduces the risk of both artificial fragmentation and excessive graph consolidation. The resulting graph structure, therefore, appears to represent a compromise between flexibility and structural stability. The method avoids forcing spectrally divergent states into a common trajectory while simultaneously preserving continuity among states connected by coherent directional change.

4.8. Consensus Structure and Uncertainty of Landscape Interpretation

The consensus analysis was used to reduce the complexity of multiple matching-based supercluster solutions and to represent them within a common spatial framework. This step was necessary because each matching method produced its own organisation of temporal codes and superclusters, making direct comparison difficult. By aligning the supercluster solutions and evaluating agreement among methods, it became possible to identify spatial structures that were reproduced consistently and to distinguish them from areas whose interpretation depended more strongly on the matching procedure. At the same time, the consensus map should not be interpreted as a single definitive spatio-temporal model of the landscape. Rather, it represents an integrative approximation that summarises common patterns across alternative models, while each individual matching solution retains information about a particular aspect of landscape dynamics. The aligned superclusters also provided a basis for quantifying the temporal behaviour of landscape units, including the relative contribution of seasonal and directional-trend components. In this sense, the procedure simplified the highly detailed set of cluster trajectories into a smaller number of interpretable spatial units, while preserving information about their dominant temporal dynamics. However, the differences among methods remain ecologically meaningful. Lower agreement was concentrated mainly in transitional zones, where relief, moisture conditions, vegetation structure, and disturbance effects interact most strongly. These areas are heterogeneous and dynamically unstable, so their ambiguous assignment reflects the complexity of landscape transformation rather than only classification error. Thus, consensus strength can be interpreted as a measure of interpretative stability. High-consensus areas correspond to robust landscape entities whose spectral, spatial, and temporal identities are reproduced across alternative matching approaches. Low-consensus areas indicate ecotonal or actively transforming zones where several developmental trajectories may overlap or where landscape-cover states have not yet stabilised. Therefore, disagreement among matching methods becomes a spatial indicator of transformation fronts and helps identify areas where post-catastrophic reorganisation remains most active.

4.9. Implications for Remote Sensing of Post-Catastrophic Landscapes

The proposed approach appears particularly suitable for the remote sensing of rapidly changing landscape systems, where temporal continuity is only partially preserved. Such conditions are typical for post-disturbance succession, floodplain restructuring, wetland transformation, fire scars, shoreline dynamics, and reservoir drawdown zones, where landscape-cover states may change faster than conventional classification frameworks can accommodate. A major advantage of the approach is that it does not require a fixed, predefined legend of classes. Instead, temporal matching is based on continuity of spectral trajectories and geometric relationships between states, allowing the system to accommodate the emergence of previously absent configurations. This is especially important under post-catastrophic conditions, where entirely new habitat structures may arise during ecological reorganisation. The results further suggest that temporal matching in dynamic systems should not be interpreted as strict preservation of invariant classes. In many cases, continuity is partial and trajectory-dependent rather than categorical. The proposed framework, therefore, provides a more flexible representation of landscape dynamics, in which transformation, divergence, and emergence of new states can be incorporated into the temporal structure of the analysis rather than treated solely as classification error.

4.10. Limitations and Future Directions

Several limitations of the proposed framework should be considered. The analysis depends on the structure of the common ordination space. Although PCA provides a convenient representation of dominant spectral gradients, nonlinear ecological trajectories may become distorted after projection into a linear component space. Consequently, the geometry of local manifolds may differ from the simplified structure represented in the ordination.

The approach remains sensitive to temporal gaps in the satellite series. Long cloud-induced interruptions may exclude critical phenological phases and artificially increase the apparent distance between neighbouring temporal states. This problem is particularly important during periods of rapid spring development, when spectral trajectories may change substantially within short intervals.

The resulting geometry depends strongly on the clustering procedure itself. Because spectral trajectories are reconstructed from independently derived clusters, different clustering strategies may alter the shape, orientation, and continuity of the resulting temporal structures. The interpretation of geometric correspondence, therefore, partly reflects assumptions embedded in the initial partitioning of spectral space. Future development of the framework may involve nonlinear manifold representations that can better preserve complex temporal trajectories. Other promising directions include trajectory-learning approaches, graph-based or graph neural methods for temporal correspondence analysis, and probabilistic matching frameworks that explicitly model uncertainty in state transitions. Integration of optical imagery with SAR data may further improve temporal continuity under cloud-contaminated conditions and increase sensitivity to vegetation structure and surface moisture.

5. Conclusions

Unsupervised classification is particularly important for monitoring dynamic landscapes, especially post-catastrophic systems, as it enables the detection of newly emerging surface states without forcing them into predefined categories. Within this approach, temporal matching becomes a critical methodological component, as it determines whether independently derived clusters can be linked into coherent trajectories of landscape-cover change. In addition to conventional matching methods, this study evaluated geometrically oriented techniques based on the shape, orientation, and axial structure of spectral point clouds. The combined set of tests demonstrated that these approaches provide a reliable basis for temporal matching by preserving correspondence among landscape-cover states, maintaining temporal stability, and achieving competitive matching accuracy under conditions of phenological and interannual variability.