Spatio-Temporal Urban Change Mapping With Time-Series SAR Data

The strong urbanization impetus of developing countries leads to various urbanization phenomena such as building constructions, reconstructions, and demolitions. It is desirable to monitor and recognize these intraurban changes by utilizing temporal and spatial information in an automatic way. This may be useful, for example, to timely update urban information databases. The aim of this work is, therefore, to automatically extract first, and further recognize, change time series in sequences of SAR data with high-frequency acquisition. Specifically, SAR time-series segmentation and unsupervised classification are combined together to recognize areas with the same urban change pattern, by fully exploiting both the temporal and spatial dimensions. Experimental results in a fast-growing Chinese city show that the proposed approach is effective and able to characterize temporal patterns due to different kinds of intraurban changes.


I. INTRODUCTION
I N THE last years, there has been an explosion of interest in mining time series of remotely sensing databases thanks to the increasing accumulation of available historical data. Exploiting Landsat, MODIS, and other optical data series, a variety of methodologies have been proposed to obtain, for example, temporal characteristics of the vegetation and extract information about its phenological state along the timeline [1], [2], [3], [4]. Even if the revisit time has decreased from 30 days for Landsat to 10 days for Sentinel-2, it may take a few weeks or months to acquire cloud-free images. On the contrary, synthetic aperture radar (SAR) is a microwave active system able to work in all weather and time conditions. Spaceborne SAR sensors, thanks to their properties, provide uninterrupted data series, ideal to monitor change patterns and not only changes. Unfortunately, few works have been published in the field of quasi-continuous temporal observations using radar series, due to the lower availability of these data in the past. Luckily, and thanks to Sentinel-1 A and B, this is not a problem anymore, and SAR time series have now a higher acquisition frequency than Landsat or MODIS optical data series.
In parallel with what happens in SAR interferometry [5], [6] to monitor terrain or structure displacements, in this work, we would like to explore the use of SAR data series, both in amplitude and phase, to detect and recognize built-up temporal change patterns. Urban changes using SAR have been widely studied and several works were devoted in the past to spatially map changes due to urban expansion, urban land cover evolution and environmental degradation, but mostly using bitemporal images [7], [8], [9]. Very few studies have instead focused on continuously monitoring built-up changes through SAR time-series data with high frequency acquisition. Since high temporal-resolution SAR data are now available, a lot of approaches originally designed for change detection in optical data series were transferred with proper modifications to radar series, such as thresholding, differencing, segmentation, trajectory classification, regression, and decomposition [1], [10], [11], [12], [13], [14]. More recently, urban changes using SAR series have been considered in [15], where the temporal coefficient of variation introduced in [16] has been included in a powerful Markov random field (MRF) framework to characterize built-up land change from SAR stack. The MRF approach, paired with the homogeneous pixel selection algorithm in [17], allows the reduction of false positives in the change detection algorithm at the expense of a computationally intensive approach, tested on a limited geographical area. The use of similar information for wide area analysis is still an open issue, as well as the recognition of the change patterns in the temporal domain, fully exploiting the spatial information used to detect the change locations. Such rich information extraction would allow to timely update urban information databases and urban building geospatial layers, with change detection results, which is significant and necessary, especially in fast urbanizing countries [18].
Therefore, this article aims at designing a novel procedure to mine built-up changing patterns with the joint use of temporal and spatial properties, starting from high-frequency Sentinel-1 (S-1) SAR time series. Specifically, the contributions are as follows: 1) a coarse-to-fine approach applicable to wide geographical areas to estimate built-up change patterns, fully exploiting the contribution of temporal information to the recognition of spatial patterns of built-up changes; 2) a temporal change pattern recognition technique, able to characterize which type(s) of changes have happened, and when. As shown in Fig. 1, the developed methodology consists of five steps: the coherence and amplitude series are first extracted This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Fig. 1. Complete processing flow of the proposed approach consisting of five steps. Z refers to the Level-1 SLC series from which the coherence and amplitude series Γ and A are extracted. The temporal moments include the expected value E(·), the variance in time Var(·) and the ratio R(·) = Var(·)/E(·). The classes C c , C nc , C tv , and C uc refer to change, nonchange, time-varying natural land cover, and urban change; analogously, the time-series segments S nc , S tv , and S uc . VV and VH refer to the polarization. from the single look complex (SLC) series. Then, their temporal moments are computed. Steps 3-5 represent the core of the algorithm and are detailed in Section II-A-II-C. In particular, only the temporal moments of the coherence series are used for coarse detection of spatial changes (step 3) to extract urban changes and discard natural ones and unchanged areas by means of a decision tree. These moments are also employed in step 4 for the temporal segmentation and classification, i.e., to divide time series into subsequences and identify within them the actual built-up change segments. Both amplitude and coherence series are finally used in a segment clustering stage (step 5) for change-pattern discrimination. The final aim of the last phase is to achieve a finer mapping of urban changes as well as to detect different change patterns of urban activities through a hierarchical clustering strategy.
The rest of this article is organized as follows. Section II describes the developed methodology, detailing the steps of the algorithm, i.e., coarse spatial change detection, time-series segmentation, and clustering for change-pattern discrimination, while Section III presents the results obtained by the proposed method when using SAR scenes collected over Wuhan, in the central part of the China. Finally, Section IV concludes this article.

II. METHODOLOGY
We consider a sequence of S1 Level-1 SLC products with the resolution of 5 x 20 m, acquired in the interferometric wide (IW) swath mode. This time series is defined as where t is the number of the images within a series and Z k is a complex vector. We extract from this sequence the following two series: 1) an amplitude series 2) a locally estimated coherence series Γ = {γ 1 , γ 2 , γ 3 , . . ., γ t−1 }, whose generic element γ k is defined as (1) Note that γ k can be regarded as a similarity measurement between adjacent temporal images at time k and k + 1 and is computed in a M × N window of L samples. Moreover, γ k can be viewed as a local estimate of the module of the interferometric coherence γ in the following equation: It has been proved in [19] that the estimated coherence overestimates the actual coherence, especially when the value of the latter measure is very low. Although this bias can be attenuated by selecting a large window, this choice would cause a loss of resolution. Thus, we choose a relatively small window to estimate the coherence-in this work, a 3 × 10 pixels' window in the azimuth and range directions, respectively. The local estimation error will be compensated by using also amplitude series in the subsequent steps. In addition, the estimation bias from the flat earth phase component and the topographic phase has been removed using a digital elevation model (DEM) by means of the Sentinel Application Platform by European Space Agency. Finally, the variability of temporal decorrelation and seasonal volumetric decorrelation are used to differentiate between changes from natural distribution targets, urban buildings, etc. Since temporal decorrelation effect is proportional to the time span, the interval is kept at 12 days to ensure that each coherence estimate is consistent between adjacent SLC pair. Although the coherence Γ is a powerful parameter for change detection, following [20], we consider the two series A and Γ since they are both indicators for change. The coherence Γ is more sensitive to changes compared to the first-order difference of amplitude series, as experimentally tested in [20]. At the same time, A series data, with finer spatial resolution than locally estimated Γ, has relatively high values due to double-bounce backscattering from man-made structures and buildings in urban areas. Thus, the combination of coherence and amplitude helps to exclude time-varying natural landcovers, such as natural vegetation and crops. This operation is also beneficial to inspect temporal and spatial change patterns only for urbanized areas, and improves the efficiency of the proposed approach when exploring changes in wide areas. In particular, change areas of interest for urban area monitoring are estimated from the coherence series Γ, leveraging on the difference between low temporal coherence from time-varying distributed landcovers (e.g., natural vegetation, water bodies, and crops) and relatively high double-bounce backscattering of man-made structures and buildings.
The main idea of this research is to implement an unsupervised detection method to extract peculiar changes. Unfortunately, as described in [19], it is not possible to statistically disaggregate the contributions of the different classes from the distribution of the generic γ k . Analytic distribution models should instead be considered to approximate this statistical model. To this aim, we note that a Gaussian mixture model (GMM) is often applied to parameter estimation of multiple mixture statistical models and is used to discriminate among multiple classes in an unsupervised way. Accordingly, we employ this model to determine where the changes occurred. The distribution representing the entire coherence series in time, P (Γ), is indeed composed of nonchanged urban areas, time-varying natural landcovers, changes from urbanization and urban activities, and so on, each one with its own fraction.
When applying a GMM, the number of mixture members should be specified from the beginning, but the number of temporal patterns is difficult to be determined a priori. Therefore, as better detailed in Section II-A, the first step of the proposed procedure is the reduction of the complex n-component Gaussian mixture into a binary Gaussian mixture model by extracting the appropriate features. Since using directly the entire image sequence sample space would imply a large amount of computation to infer the model parameters of both the stable areas and the changing areas, especially when the image covers a large area, the proposed binarization, aimed at extracting only the possible changes, is applied in a patch-wise manner.
Once this task is completed, the recognition of the Γ sequence's time segments is performed in each patch related to the changes, which are then separated into patterns corresponding to different changes by means of an analysis in time domain. Finally, as last step of the proposed methodology, a tensor-wise procedure is employed to infer a more accurate statistical model. Specifically, combining the spatial coarse boundaries and the time segments determined previously, the spatio-temporal tensor containing the urban changes is extracted. Locally dependent statistical inference models built on these tensors, such as the aforementioned GMM, the minimum cross entropy thresholding, and variational autoencoder (VAE), are used to determine the spatial and temporal boundaries of fine urban changes. The overall procedure is represented in Fig. 1 and is detailed in the following sections.

A. Coarse Spatial Change Detection
The first step of the proposed approach allows to avoid the huge computational burden that would be required by time-series analysis for the whole image. It aims at a first discrimination between plausible changed and unchanged areas. Accordingly, a temporal averaging operation along time is performed. Note, however, that a word of caution is needed. Temporal averaging on SAR time series is often used for speckle denoising and is performed under the assumption that no change occurs within a time series. If applied to a sequence including changes, this operation tends to degrade the change intensities that occur at one or more instants of time along the series. Moreover, this makes it less easy to discriminate abrupt changes from naturally time-varying landcovers, like crops. It is, therefore, necessary to add some additional information about the differences between the series values and the averaged one.
Since the temporal average of the coherence series Γ is an estimate of its expected value E(Γ), its variance in time Var(Γ) is also computed. While the combination of E(Γ) and Var(Γ) may not be enough to discriminate between different types of changes, it is sufficient to recognize the changes from the invariant areas, which is a first step toward the recognition of more complex change patterns. To this aim, we consider the ratio R(Γ) = Var(Γ)/E(Γ). This quantity can be regarded as a normalized change intensity that contains a fraction of changing areas.
The statistical distribution of E(Γ), Var(Γ), and R(Γ) in the spatial domain is shown in Fig. 2 for a sample test area. The figure reveals that a bimodal Gaussian distribution is a sufficiently accurate model for the histograms of E(Γ), Var(Γ), and R(Γ). This is reasonable under the assumption that R(Γ) is the fusion of two coherence sets: one from time-varying change landcovers, such as water surfaces, forests and other natural covers, and another one from invariant landscape elements, e.g., unchanged built up areas.
Scatter plots in Fig. 4(b) refer to the divergence of different landcovers corresponding to Fig. 4(a). Among the considered temporal moments in the figure, R(Γ) provides a better discrimination between the mostly invariant built-up areas and time-varying natural landcover classes or urban areas subject to changes. This is exemplified by its scatter plot. Indeed, natural landcovers have a low temporal coherence because of weather short term events, such as wind and rain. Urban areas subject to  change activities present a similar pattern and a relatively high ratio due to a large variability in coherence, even if the average coherence value is larger than for natural targets. Building facilities with no change, on the contrary, present a large temporal coherence and very low variability, which translates into very small R(Γ) values.
Since the distribution of change and no-change areas is often very unbalanced in the spatial domain, it is not immediately feasible to recognize changed areas by a unique threshold for the whole scene. Thus, in this work, a segmentation procedure is applied to subdivide the scene into patches belonging to nochange built-up areas or to areas undergoing a change, which might include both natural and urban changes. The procedure works as follows.
1) The first step is a top-down quad-tree splitting algorithm as in [21], which generates homogeneous patches according to the value of R(Γ). To this aim, the aforementioned statistical measures are used. Specifically, to generate patches with an unimodal distribution, in addition to satisfying Ashman's D condition [22], the two distributions p 0 and p 1 that are extracted from a potentially bimodal, one must satisfy the additional condition that p 0 (μ 1 ) ≥ p 1 (μ 1 ) when σ 0 > σ 1 , where μ 1 and σ 1 are the mean and standard deviation of p 1 . Only in this case, the local p(R(Γ)) is considered to be bimodal. As it is shown in Fig. 3, if this condition is not satisfied, the mixed histogram is better represented by an unimodal distribution. The final patches, i.e., the leaf nodes of the quad-tree, are eventually used as input to the next step. Once the patches with unimodal distribution have been extracted as the leaves of the quad-tree, those leaves at level depth are merged to form a bimodal patch at the previous level depth − 1. If the new patch remains unimodal, the procedure is iterated at the upper level. The extracted bimodal patches, with different size in different parts of the scene, are collected to estimate the two statistical components, labeled as "invariant" and "mixed change" areas, using the GMM clustering algorithm. 2) After that the "invariant" built-up areas are excluded, this step focuses on a better characterization of the "mixed change" areas. Since mixed change can be visually separable by means of E(Γ) and Var(Γ), as shown in Fig. 4(b), another bimodal distribution fitting is applied to one of these quantities. The aim is to discriminate between urban and natural changes. The "urban change" component can be thus extracted, under the assumption that they have larger values of both E(Γ) and Var(Γ) than those of naturally time-varying landcovers. It is worth noting that the component-balancing strategy in the aforementioned step 1) still applies to both E(Γ) and Var(Γ). Finally, these coarsely mapped urban changes are objectified and used as input for the following mapping procedure, aimed at automatically determining the time period in which the change occurred within the entire time series.

B. Time-Series Segmentation
Once the areas with possible urban changes have been detected in the whole scene, the actual time frame of the urban change should be extracted. However, the aperiodicity is one of the most distinctive features of urban construction or inner urban activities with respect to other periodical phenomena such as farming or, more in general, vegetation phenology. This means that in many cases, the time interval over which the urban change occurs is only a small portion of the whole time sequence, and this portion should be singled out from the rest to better characterize and recognize different change patterns.
Although one pixel is spatially identified as potential urban change, it may include both natural and urban changes within the entire time series. To efficiently determine the time interval in which urban change occurs, before applying any clustering, as first step an averaged coherence series is generated by the spatial aggregation of the sequences that are all spatially identified as urban change and objectified in the previous step. Then, this sequence is temporally segmented by detecting its change and turning points, which correspond to relatively low and high coherence levels.
Specifically, segment sets S = {s 1 , s 2 , s 3 , . . .s m } are generated and the procedure in [23] is performed as summarized as follows. As a result, the extracted segments are assigned to three sets, namely invariant built-up S nc , time-varying natural landcovers S tv , and urban changes S uc . Only the last set is considered for further analysis.
1) The procedure in [23] starts with the shift-window algorithm [24] to generate "turning points" (points where the sign changes) of the difference series Φ, whose ith element is defined as φ i = |γ i −γ uc |, with i = 1, 2, . . ., whereγ uc is the expectation of the coherence for points labeled as possible urban changes. As an example, Fig. 5 reports the GMM distributions with three classes computed from E(Γ) in one patch. The series Φ is designed to detect valley points belonging to the urban change C uc , and peak points belonging to C nc and C tv classes. Then, these peak and valley points are used to detect the turning points at which C uc class transforms into other classes, or vice versa. The shift window size is set to 10 (only one peak is generated from ten adjacent points) and the deviation fromγ uc is considered significant if it is greater than 2σ, to ensure that no dense turning point sequence is produced. Finally, the regression line between peak and valley points in Φ (including the starting and ending points of the series) is used to find its turning points as those that are farther away from these regression lines than the others. 2) Once the sequence has been segmented, a classification is performed to select the segments of urban change, denoted as S uc . Since the conditional probabilities P (Γ|c i ) for the classes C = {c tv , c uc , c nc } are available, they are used to classify a generic segment s using its associate γ values, i.e., {γ m , γ m+1 , . . ., γ m+n }: 3) Finally, adjacent segments belonging to the same class are merged. The idea is to ensure the integrity of the different spans of urban activities, despite the fact that they can be fragmented by the previous steps. Moreover, to avoid the generation of very short segments, like those referring to the temporary shutdown during building construction that may be classified as no change, a window threshold is set up to limit the minimum segment length to five points. Segments shorter than the threshold are merged with the closest one.

C. Segment Clustering for Change-Pattern Discrimination
As a result of the procedure described previously, spatiotemporal segments are generated in the form of 3-D tensors. In this part, a fine-mapping procedure for urban changes is proposed on the one side to refine the previous coarse spatiotemporal mapping results, and on the other side, to detect different change patterns of urban activities.
Regarding the first point, the segment set including all intervals with no changes, S nc = {S i nc |i = 1, 2, 3, . . .}, can be utilized to estimate the boundaries of temporarily stable urbanchange targets by a two-component GMM clustering with recomputed values of E(Γ S i nc ) and R(Γ S i nc ). Accordingly, a set of boundaries B = {B i |i = 1, 2, 3, . . .} that correspond to temporary stable urban subareas before or after a change, or even in a stasis period within two changes, can be extracted and compared to identify a change of expansion/contraction or as a change inside an urban subarea.
Regarding the detection of change patterns, the final result cannot be limited to an entire quad-tree patch along the whole sequence, but must be spatially and temporally refined. Therefore, the previously extracted statistical model cannot be directly applied to each coherence image, because they were generated by means of a temporal and spatial average. Additionally, in this analysis, we do not know a priori whether we should consider two or three clusters. Finally, since urban changes involve many possible phenomena, it may be worth at this point to consider the intensity in addition to coherence, for a better characterization of punctual changes, e.g., changes to single buildings or small building blocks. To cope with the time-dependent and individual coherence modeling problem of the sequence, a Gaussian hidden Markov model (GHMM) is introduced to map gradual urban changes among segments [25].
When applying the GHMM method, the mixture of C uc , C nc , and C tv class is considered in this step. In case of two-component clustering for coherence modeling, the segments classified as S nc in time domain are spatially separated as C nc and C tv . When turned into three-component clustering for segments S uc , the amplitude difference sequence ΔA is introduced to enhance the new component of urban change class C uc . Indeed, the coherence shows high capability in discriminating between class C c and class C nc , but is weak in distinguishing C uc from C tv . In particular, the last coherence image in the previous S nc segment will be used to initialize the priori and Gaussian parameters of the hidden Markov model (HMM) with a zero prior probability for class C uc . Then, the log-ratio amplitude difference sequence ΔA is added.
At this point, to perform the clustering, and then, recognize a pattern in the sequences, two issues must be addressed. The first problem is the standardization of the change intensity along the time span, since the intensity of changes and the time span are not consistent across different locations of the subtensor. In order to compare these segments with each other, intensity characteristics of differently long sequences are roughly "equalized" by dividing the intensity of the change by the time span.
The second issue is the sequence alignment, since traditional unsupervised classifiers require to have features with the same dimension in input. Indeed, similarity between sequences is the base of clustering. However, unaligned time sequences with different duration, even in case of similar trends (e.g., in case of building construction), are likely to be considered as very different unless they are compared after resampling and scaling into those of same duration. Especially for nonperiodic urban activities, the presence of large shift variance, warp variance, and noise variance in the time series leads to large differences among similar patterns. To cope with these issues, the dynamic time warping (DTW) method is introduced to measure the similarity between two unaligned temporal sequences [26].
The sample space consists of spatially averaged sequences calculated during S uc segments in different locations. Specifically, the maximum temporally stable boundary B generated in the previous step is used to compute the average value for each tensor slice in each S uc segment in the form of a subtensor. Moreover, in order to include the starting and ending states of the urban change sequence, this sequence is extended by 3 time points on both sides. Then, the DTW method is applied to align the sequences. Finally, the clusters are separated using the K-means algorithm [27].

III. RESULTS
In order to validate the effectiveness of the proposed method, 93 SLC SAR scenes (ascending orbit) were collected over the city of Wuhan (see Fig. 6). The approximate viewing incidence angles of the series range from 40.6°to 40.8°, and the series covers from March 2018 to March 2021. Due to the fast urbanization speed of the city in recent years, even during the lockdown, it is feasible to obtain typical urban change patterns in these four years. Fig. 6 shows the mean complex coherence and intensity of the SAR sequence in the study area. It is clear that the former parameter is more sensitive to urban built-up than the latter one. Indeed, the bright areas in Fig. 6(b), which represent a high degree of coherence, highlight urban built-up areas with respect to natural landcovers.

A. Coarse Urban-Change Mapping
As mentioned in Section II-A, to cope with unbalanced mixedcomponent clusters, a quad-tree splitting strategy is exploited to generate bimodal patches. The patches for an area in Wuhan city center are shown in Fig. 7. The red and green rectangles, generated at the end of the top-down splitting phase are examples of homogeneous areas with unimodal distributions of R(Γ). Conversely, the blue patch, obtained after the bottom-up step, has a distinct bimodal distribution. The statistical distributions for R(Γ), E(Γ), and Var(Γ) in these three patches are shown in Fig. 7(d)-(f), respectively.
Subsequently, to unmix the clusters, bimodal tiles are decomposed using the GMM clustering algorithm. As shown in Fig. 8(a), the total mixed histogram is decomposed into "no change" and "mixed-change" components. Then, the mixedchange component is decomposed into urban changes and timevarying landcovers by applying again the GMM clustering algorithm, this time using E(Γ) as input [see Fig. 8(b)],which is more stable than Var(Γ) on mixed Gaussian modeling.
The coarsely mapped changes are given in Fig. 9. The areas in white in all the figures are those detected as no-change built-up areas C nc by decomposing R(Γ). After that a two-component GMM is applied to decompose E(Γ), urban and natural changes are extracted and shown in red and black, respectively.
To demonstrate the validity of the method, ground-truth samples were manually collected inside and outside the settlement areas as shown in Fig. 9, and results are shown in Table I. Considering the degradation of the change information due to  the temporal and spatial averaging, the validation is performed qualitatively at the object level. Specifically, a morphological filter is applied to the results in Fig. 9 to discard areas smaller than a single building. Some examples of the tested area are shown in Fig. 10. Changes like building construction activities or building demolition activities are detected, but the boundaries of these detected changes must be refined. At the same time, Table I shows that class C nc obtains higher detection accuracy  Fig. 6(b) (Fig. 7(a)-(c) refers to tiles 1, 2, and 3, respectively). (d)-(f) Histograms for R(Γ V V ), E(Γ V V ), and Var(Γ V V ) in three patches identified by a red, a green, and a blue square in the area depicted in Fig. 7(a). Fig. 8. On the left, an example of the decomposition of the distribution for one tile into "no change" in blue and "mixed changes" in yellow using R(Γ); on the right, the "mixed changes" in yellow histogram is decomposed into "urban change" in red and "time-varying natural change" in green using E(Γ).
because it is not affected by the temporal-averaging operation, while classes C tv and C uc have mutual confusion and relatively lower accuracy.

B. Segmenting the Time Series to Obtain Fine Urban-Change Maps
According to Section II-B, the next phase is the coherence time-series segmentation to characterize urban change areas. Some examples of segmentation and segment classification results are provided in Fig. 11, where each series plot corresponds to a subplot in Fig. 10. Unlike the cyclical characteristics of natural land type changes, asynchronous features in urban change patterns can be used to characterize periods of faster or slower activity.
For instance, shorter segments are more likely to be detected in time intervals when change activity increases. To this aim, the number of segments detected working on coherence sequences is significantly larger than those obtained from the backscattered amplitude sequences, which is a further evidence that coherent features are more sensitive to urban activities. Specifically, a backscattering sequence is able to show a clear trend of variation when reflecting urban high-rise construction activities, especially when a strong dihedral angular scattering variation happens. However, coherence holds a better potential for cases where weak amplitude changes are detected, e.g., in case of low-rise urban building activities.
It is also important to note that coherent sequences have good continuity and smoothness in the change trend, while backscattered sequences have some random fluctuations at adjacent time points. This advantage of coherent sequences is due to the use of spatial averaging in a window of fixed size during its computation.
The separation of the overall sequences into segments is based on break points, whereas peak and valley points indicate the increasing or decreasing trend of the change. For instance, in Fig. 11(b), corresponding to Fig. 10(b), the presence of two peak points in the sequence suggests a temporary pause in the construction activity.
After estimating the Gaussian parameters of the classes C = {C tv , C uc , C nc } by the GMM clustering algorithm, the segment classification results obtained by using the maximum likelihood method are shown in Fig. 11. Red, blue, and green colors identify segments assigned to different classes. For instance, segments with high coherence values are recognized as C nc and shown in blue. However, the short segments in Fig. 11(a), (b), (d), (f), and Fig. 9. Coarsely mapped "mixed changes" for three tiles through GMM clustering applied to E(Γ) (note that no-change areas have been detected using R(Γ) and are colored in white). Fig. 9(a)-(c) refers to the locations of tiles 1-3 in Fig. 6(b), respectively. The classes C uc , C nc , and C tv refers to urban change, nonchange, and time-varying natural landcover, respectively. The validation samples concern the ground-truth samples manually collected and used for the validation of the coarse mapping step of the proposed methodology (refer to Table I). (h), classified in this group, are peculiar. Indeed, the presence of coherence peaks in each of these segments indicates that the activity stops briefly in winter, in an interval of time around the Chinese New Year festival.

C. Segment Clustering and Change Pattern Mapping
Once the segmentation and classification are applied to all series, the sequences (in the form of tensors) are extracted and clustered. However, before clustering the sequences, the exact boundary of the change needs to be determined, as mentioned in Section II-C. After obtaining these boundaries, an objectoriented temporal clustering analysis is performed, as opposed to a pixel-level temporal clustering. This is expected to be advantageous, especially for noisy SAR images.
Since it is difficult to determine the number of clusters, a hierarchical approach is applied. In Fig. 12(d), the class-center sequences of Γ vv are obtained by calculating Fréchet mean of variable-length sequences clustered as one class. The length of the Fréchet mean sequence of one class is set to the averagedlength of all sequences belonging to the class. Leaf clusters, e.g., clusters 1 and 2, belonging to the same parent cluster, share similar sequence patterns. Analogously, clusters 1-4 have relatively similar temporal fluctuations, while clusters 5-8 have common trends in coherence and show ongoing building construction activity. In short, the differences between change patterns are mainly attributed to mean level of coherence, change trends, and cyclical fluctuations.
To effectively decipher these change patterns, we may refer to the patches with urban change shown in Section III-B. Patches Fig. 11. Graphical representation of coherence and amplitude time-series segmentation results. A maximum likelihood method is used to classify the segments and assign them to the classes C uc , C nc , and C tv , i.e., urban change, nonchange, and time-varying natural landcover identified with red, blue, and green color, respectively. (a) Series segmentation for Fig. 10(a). (b) Series segmentation for Fig. 10(b). (c) Series segmentation for Fig. 10 11 and 12(d). This pattern of variability is fluctuating, with an initial decrease, and then, an increase in coherence. There is a local peak of coherence in winter, after which another fluctuation is ushered in, and then, the final coherence reaches its highest level. Combined with Google Earth History Images, this pattern shows the seasonal variation in building construction activity, with overall construction activity taking several years to complete.
Patches (b) and (c) in Fig. 10 are assigned to cluster 7. Specifically, cluster 7 shows a large change in coherence throughout the sequence, with an initial coherence value much lower than the one of cluster 1, while at the end becomes much higher than the one of cluster 1. This pattern of variation suggests a short-term construction activity with a faster rate of construction.
The first subsequence of patches (f) and (h) are assigned instead to cluster 8. In Fig. 12(d), cluster 8 indicates an initial coherence level below 0.5 and a rapid increase in coherence after a short period of stabilization. This type of variation pattern shows only a nonseasonal construction activity, which can be considered as a fragment of a multiyear and cyclical construction activity.  Fig. 12(a)-(c). Fig. 12(a)-(c) refers to the locations of tiles 1, 2, and 3 in Fig. 6, respectively. Patch (e) belongs to cluster 3, which exhibits the lowest level of coherence in V V polarization. Considering that the spatial locations of cluster 3 in tile 1 [refer to Fig. 9(a)] are mainly concentrated in the natural landcovers, it potentially suggests that there are consistent changes involving natural covers and urban built-up structures. Conversely, cluster 4 is distributed within the built-up area and has a higher mean value of coherence sequence than cluster 3, indicating a multiyear, cyclical pattern of change in construction activity.
Finally, cluster 2, which corresponds to path (g), shows a decreasing trend of coherence, indicating a type of urban change of demolition activity. Demolition activities are likely also to be present in building reconstruction indicated by clusters 5 and 6.
In summary, analyzing the clustering results in Fig. 12(d), multiple urban change patterns can be effectively identified by the coherence's mean value, the fluctuation, the change trend and the range of sequences. Clusters 1-4 present multiple fluctuation, possibly due to a mix of built-up and natural landcovers, and a relatively low level of coherence, as well as a high coherence peak level in the short winter period. The coherence reduction of cluster 2 is very likely to correspond to building demolitions. On the contrary, clusters 5-8 present a significant coherence increase, which means that constructions are close to be completed in a relatively short period.

IV. CONCLUSION
In this article, a framework for temporal and spatial change pattern detection applied to SAR time series is proposed. The developed method has been tested in Wuhan City since this area experienced a fast urbanization speed in the last four years, with typical urban change patterns. Experimental results confirm that the proposed algorithm is able to discriminate among different urban change patterns in the time domain. The contributions of this research can be thus summarized as follows.
1) The proposed method employs quad-tree segmentation to identify local patches with bimodal properties, and then, applies to these patches a Gaussian mixture decomposition. The results reveal that the balanced bimodal split shows a high robustness and separability.
2) The developed solution effectively addresses the problem of jointly estimating the statistical parameters of change categories in the time and space domains. For the change information decay caused by time-domain averaging, the spatial and temporal tensor without decay is selected by a two-step method, assuming that temporal and spatial domain share the same statistical parameters of the model, and then, the parameters are estimated. The results show that the urban invariant class can be estimated without information loss using invariant segment, and the model parameter optimization of the urban change class can be effectively assisted by using the prior information provided by the unsupervised model to achieve the parameter optimization of the HMM model. 3) The proposed framework explores nonsupervised change patterns for nonperiodic urban activities. The developed algorithm successfully decomposes urban changes of different lengths, intensities, and randomness by splitting time series using turning points. The solution also accurately detects and explores changes in urban activities using a time-dependent model. An open issue that we aim to face is the computational time reduction of the algorithm. Indeed, three days are required to analyze one scene, of which, however, two days are needed as preprocessing to extract InSAR coherence and amplitude timeseries data. Currently, we exploit cloud computing platform to overcome this issue, anyhow future research directions include the optimization of the procedure to reduce its burden computational operations. Another future work direction involves a more accurate coherence estimation and the analysis of the temporal decorrelation effect. Indeed, it has been proved in [19] that estimated coherence overestimates the actual coherence, especially when the value of the latter measure is very low, which may increase the ambiguity between distributed timevarying targets and urban changes. Therefore, a more accurate coherence estimation is needed for distributed targets with improved estimators, e.g., the coherence estimator in DSIpro [17]. Finally, although from a preliminary study our results are not strongly affected by the building orientation, we plan to conduct a detailed analysis to examine this aspect.
The data and the methods involved in this article can be accessed at the dedicated github repository (https://github.com/ lastrye/Spatio_temporal_change_pattern).