Practical approach for synthetic aperture radar change analysis in urban environments

Abstract. Change detection using remote sensing imagery is a broad and highly active field of research that has produced many different technical approaches for multiple applications. The majority of these approaches have in common that they do not deliver any detailed information concerning the type, category, or class of the detected changes. With respect to the extraction of such information, recent research often suggests that a land use classification is required. This classification can be accomplished in an unsupervised or supervised way, whereas the practicability of both strategies is more or less limited by the usage of reference or training data. Moreover, expert knowledge is needed to arrive at meaningful land use classes. An approach is presented that overcomes these drawbacks. A time series of synthetic aperture radar amplitude images is considered, enabling the detection of so-called high activity objects in urban environments. Such objects represent the basis of the investigations and denote the input for unsupervised categorization and classification procedures. The method supports even the unexperienced user in learning the actual information content leading to the capability to define a suitable scheme for change classification. Tests carried out on two different datasets suggest that the method is both practical and robust.


Introduction
In general, when speaking about change detection using remote sensing imagery, passive and active sensor systems can be mentioned for data acquisition. Especially for change detection and monitoring issues, where a regularly sampled data basis is necessary to produce meaningful analysis results, active sensors have some significant benefits. For example, current synthetic aperture radar (SAR) satellite sensors, such as the German TerraSAR-X (TSX; operating since 2007) and TanDEM-X (TDX; since 2010) missions, 1 illuminate the Earth using microwave radiation at a wavelength of about 3 cm. Due to their active transmission of microwave radiation, SAR systems can operate during day and night and independently of clouds, fog, and dust, providing 24/7 monitoring capabilities. 2,3 Focusing on the TSX and TDX mission, the highresolution Spotlight mode HS300 enables the acquisition of images with a geometric resolution of less than 1 m, 4 which allows detailed analysis of urban areas.
Change detection from remote sensing data is commonly described as a process that identifies changes of the Earth's surface or of objects and structures placed on it. [5][6][7] Focusing on land cover characteristics, the two items "conversion" and "modification" can be generally distinguished, 8 where conversion denotes the complete replacement of cover types, and modification stands for more subtle changes, affecting also the character of the specific land cover. In this paper, changes of land cover are not addressed. A method is presented that deals with the classification of land use, based on manmade change objects located at different structures of urban environment. Commonly utilized change detection techniques with their respective main applications can be additionally found in Refs. 9 and 10. For example, techniques are categorized into abrupt and progressive approaches. 8,11 Abrupt approaches identify changes in bitemporal data, whereas progressive methods allow for the extraction of change trajectories, delivering information about trends and impacts of change with respect to the land cover. 12 In Ref. 13, change detection approaches are sorted by their computational complexity, and in Ref. 7, pixel-based approaches are compared to object-oriented methods.
In this study, a method for change analysis is presented, which provides detailed information on the type or class of the detected changes. 14 Here, the two terms "change detection" and "change analysis," which are often used synonymously, are distinguished as follows: "change detection" represents the process of detecting changes leading to an explanation of change objects by their locality and date. Typically, any further information concerning the category, type, or class of the detected changes is not delivered. The term "change analysis" describes techniques that include both change detection and a more detailed explanation of type, category, or class of changed objects in the data processing flow.
An example of an SAR change analysis method was published in Ref. 15, where an initial change detection is combined with a classification strategy for the extraction of detailed change information. Focusing on SAR image exploitation, change detection as the first step of change analysis can be generally subdivided into coherent and incoherent techniques. 16,17 In our study, we aim at the analysis of urban environments, where change objects of interest are often embedded in a background with very low radar cross section (e.g., vehicles on a parking lot). Hence, a low coherence even for unchanged objects result, making a clear distinction between changes and their surroundings difficult. Therefore, we focus on an incoherent strategy.
With respect to change classification as second step of change analysis, postclassification comparison can be regarded as one typical representative. 13 In general, for this strategy, a classification of the input images has to be accomplished, 18 which can be performed in a supervised or unsupervised manner. 2 As central point of these approaches, the usage of reference and training data is crucial for the final analysis results. For example, for supervised classification, suitable training data have to be available or selected by the user, who consequently has to have a detailed knowledge about the method and the scenery. Moreover, the high amount of time-consuming work limits the practicability of this strategy. 7 In contrast, unsupervised classification methods automatically perform widely without user interaction. Nevertheless, they require more or less comprehensive training databases, which have to be as meaningful as possible with respect to the specific task. Moreover, if unsupervised clustering methods are considered, suitable reference data are needed to assign physical meaning to the automatically found data clusters. 2,19 Further strategies for classification are given by machine learning methods (e.g., support vector machines). Such approaches may hold the advantage of being independent from the statistical properties of the input data, but they are limited by the fact of requiring a more or less detailed training database, which directly affects the quality of the analysis result. 7 Referring to the task of change analysis using SAR imagery, different methods exist. Out of all published methods, we discuss those that derive sufficient information on the type or class of detected changes.
In Ref. 15, intensity image ratioing and coherence-based techniques are applied to two TSX Stripmap-type SAR images to identify land use and land cover changes. The type or class of the changes is subsequently derived by an object-based classification calculated on the preevent image.
Another SAR-based strategy for a postclassification of changes is described in Ref. 20, where short-term and long-term urban changes are investigated using ERS satellite imagery. As classifier, a neural network trained by radar image features is utilized. Although satisfactory results are achieved, it is noted that the neural network architecture requires care and knowledge in designing, training, etc.
As a feature for the analysis of crop changes in rice fields, the sigma nought (σ 0 ) backscatter coefficient and its variation over time is used in Ref. 21.
In Ref. 22, temporal variations of radar backscatter are used in combination with a long-term coherence analysis to arrive at an automatic identification of land cover changes. For urban change detection, a similar strategy is applied in Refs. 23-25, where backscatter and coherence variations are utilized for the calculation of a two-dimensional (2-D) histogram. This histogram is subsequently used as the input for a thresholding method. Another time series analysis of backscatter is given in Ref. 26.
A pixel-based extraction of change frequencies along a time series of Sentinel-1 amplitude images is performed in Ref. 27. In Ref. 28, a similar approach is used to categorize changes into the classes step, impulse, and cycle. Another similar concept to this approach can be found in Ref. 29, where four classes of change are distinguished: abrupt (urban), evolutional (vegetation), no change, and periodical (seasonal vegetation change).
Another concept for the identification of changes from SAR time series is available in the form of persistent scatterer interferometry (PSI). While PSI was originally developed for the detection and monitoring of surface deformations, 30,31 recent research has discovered its utility for change detection. For instance, in Ref. 32, PSI is combined with a conventional pixel-based method to detect changes at object level (buildings) and at surface level (deformations).
In Ref. 33, SAR change analysis is performed using a method consisting of three central steps: simple linear iterative clustering 34 superpixel segmentation, feature extraction, and classification. For the last step, a supervised perceptron learning algorithm is considered, which means that expert knowledge was required for the selection of training areas for each target class.
Some previous authors performed change analysis via a combined use of SAR and optical imagery. For instance, in Ref. 35, the classification results of multiple RADARSAT images and one Landsat image are compared.
In summary, one major aspect can be observed from the information provided above: a detailed analysis of detected changes is closely related to the classification task. From this, it follows that current and comprehensive reference and training data have to exist or have to be selected by the user. For both cases, an experienced user is necessary, having detailed knowledge of the dataset and the methods being considered. 2,9,13 Focusing on the change classification techniques, the postclassification comparison is often described as a suitable method if sufficient training data are available. Moreover, in Ref. 13, it is mentioned that a new method has to be easy to handle and as robust as possible. In addition, detailed information concerning the type or class of change has to be provided.
With the paper at hand, a method for SAR change analysis is proposed that overcomes the mentioned classification-based drawbacks. In contrast to existing approaches outlined above, the presented strategy is interpreted as a method, which supports an unexperienced user in learning the actual available information content with respect to a suitable classification of land use.
As input data, time series imagery is utilized, enabling the identification of areas, which are changing very frequently. Such high activity areas are characterized by manmade targets, the socalled high activity objects (HAOs) (Sec. 3), representing the basis of our investigations.
To ensure a high degree of portability, a sparse dataset only existing of SAR amplitude images is considered. The signal amplitude represents basic information that can be directly derived from all SAR sensors. Hence, in contrast to polarimetric-based approaches, 36,37 our strategy can be interpreted as more practical-oriented. This paper is organized as follows. The concept for change analysis and the considered dataset are described in Sec. 2. A more detailed description of the workflow is presented in Secs. 3 and 4. In Sec. 5, some experimental results are shown. A validation of the method is given in Sec. 6. The main characteristics of the presented method are critically discussed in Sec. 7. A conclusion is provided in Sec. 8.

Concept and Data
As already described above, the structure of a newly developed change analysis method should be as simple and portable as possible, which means that the workflow should contain a high degree of automation. This automation also implies a high level of usability and practicability.
The developed approach consists of an initial detection and categorization step followed by a classification of identified changes ( Fig. 1; cyan-and magenta-colored). Both steps can be completed without requiring extensive reference or training data.
Solely, time series of SAR amplitude images are used as input data, enabling the identification and analysis of regions that are changing frequently. Real-world examples for such "high activity areas" are parking lots, construction sites, or collecting points such as wasteyards. We will refer to change objects identified in high activity areas as so-called HAOs. Once identified, the HAOs are the input for a subsequent categorization step, which is accomplished on the basis of a previously performed feature extraction (Sec. 3.1).
The clustered or categorized HAOs represent the central input for the subsequent change classification step. Since no further training data is applied, the categorized HAOs are referred to as data reality. This data reality is compared to a class catalog that is initially selected by the user. Within the scope of a forward-projection (FP), all HAOs are assigned to a best fitting class and the membership composition of each class is assessed (which initial HAO categories are contained in each class). In doing so, it is tested whether an initial class contains a dominant HAO category or not. In that way, some significant prior information is extracted, which is considered in a subsequent iterative optimization step (backprojection, BP).
The core of the BP step is the analysis of predefined key indicators (Sec. 4.3) that are calculated for each class of the catalog. Based on this key indicator analysis, an optimized version of the class catalog can be calculated, resulting in class definitions that match best to the information content of the data reality. If requested by the user, this procedure may be passed through again as indicated by a dashed arrow on the right side of Fig. 1. Guidance for the user's interpretation of the BP results is provided via a semantic context analysis. For this step, a flexible and adaptive data basis can be included (e.g., OpenStreetMap or land use maps). In this paper, the semantic context is not addressed. For more information concerning this step, see Ref. 14. In summary and in contrast to strategies used so far, the categorization as well as the FP and BP achieve high levels of automation while relying on a very limited amount of input data.
The input data can be subdivided into SAR imagery and optical data. An area near Stuttgart (GER) airport [ Fig. 2(a)] has been chosen to demonstrate the developed technique due to its  The time series consists of 15 high resolution (HR) TSX amplitude images, acquired in the HR Spotlight mode HS300 in ascending orbit direction. The data are imaged in HH polarization and cover a time interval of one year (November 2013 to November 2014). The scene center is located at 48.68°latitude and 9.216°longitude (WGS84 ellipsoid). Due to the HS300 mode, a pixel spacing of 0.45 m (range) and 0.86 m (azimuth) exists. A digital terrain elevation data level 2 height model was used to orthorectify the TSX scenes. The benefit of using orthorectified data is given by the fact that further geodata such as GIS layers can be directly integrated in the processing scheme. Furthermore, geometric measurements concerning the shape and size of real-world objects can be easily accomplished (Sec. 3).
For the calculation of the so-called differential attribute profile (DAP) features (Sec. 3), all SAR images were sigma nought (σ 0 ) calibrated following an approach listed in Ref. 38. The mean amplitude image of the time series, determined for visualization purposes, is visualized in Fig. 2(a). Since multitemporal averaging was applied, a reduction of speckle was achieved without affecting the spatial resolution. This limitation occurs when using the conventional multilook strategy. 39 For evaluation and land use interpretation purposes, one WorldView-2 (WV-2) image was available. It was acquired in August, 2011 and does not cover the TSX scenery completely [ Fig. 2(b)].

Change Detection and Categorization
The flowchart of our change detection procedure including a detailed view on the input and output variables is shown in Fig. 3. A set of example input and output datasets of the procedure is shown in Fig. 4.
For incoherent SAR change detection, several methods exist, among which the so-called Log-ratio (LR) scaling can be seen as a common and frequently utilized algorithm. 11,40-44 LR does not require comprehensive parametrizations, which makes it well-transferable on other sensor products or scenes. Hence, LR is well-suited to serve as the first step of our automatic processing procedure. At least one SAR amplitude image pair containing two images A t1 and A t2 has to exist. A maximum LR image, which considers both appearing and disappearing changes between the two images, is then calculated on a pixel-by-pixel bases with E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 4 8 9 The resulting ratio image A Ã mr [ Fig. 4(d)] may contain a high level of noise, consisting of bright and dark structures of different size.
Without the application of an appropriate filtering procedure, false-detections would occur, limiting the accuracy of the detection result. At the same time, the filtering procedure should perform as structure-preserving as possible in order to support a detailed and meaningful change analysis.
These requirements are met by the commonly used morphological alternating sequential filters (ASFs) combined with the connected and attribute filter principle. Connected filters represent an advancement of classical morphological operations such as openings and closings since they operate on connected components (CCs) existing in the image and do not require a predefined structuring element (SE) as it is utilized in "conventional" ASFs. 45,46 The attribute filter principle implies testing the CCs against a specific criterion (e.g., area or size). When a CC fulfills the criterion, it is preserved-otherwise it is removed. Consequently, a connected ASF (ASF cc ) utilizing the area attribute provides both a structure-preserving noise reduction as well as the possibility to control the size of the resulting change objects. In other words, an ASF cc enables the detection of differently sized changes dependent on the used area attribute.
In this study, a parametrization is chosen that leads to the detection of small-scale changes as, for instance, related to vehicles or construction sites. 14 For example, as reference for the minimum size, eight pixels were chosen, representing an average size of a European compact car. Due to the fact that all input images were previously orthorectified to a quadratic raster with 1 m × 1 m pixels, the connection between real-world and imaged objects can be directly established.
Morphological operations often require imagery that is scaled to 256 gray values. Therefore, the ratio image A Ã mr is log-transformed and scaled to this 8-bit interval. After calculating and filtering the ratio image by applying the connected ASF [ Fig. 4(e)], the detection of changes is performed. For this step, the unsupervised entropy-based thresholding method of Rényi 47 is applied. Entropy-based thresholding methods have the benefit of calculating meaningful threshold values even if the gray values of the input image follow a unimodal distribution. 48 The method of Rényi, as described in Ref. 47, calculates entropies for the image fore-and background using the information of a 2-D histogram. Such entropies are summed up for each possible pair of image gray values and the optimal decision threshold is found by maximizing this sum of entropies. The background information is simulated by the mean-filtered  In the following, the strategy for the identification of HAOs is explained. Moreover, the main features that will be used for HAO categorization are introduced.

ActivityMap, HAOs, and Feature Extraction
The bitemporal change detection procedure as described above is adapted to be applicable to time-series data by applying it successively to each available amplitude image pair. The CMs are summed up and color-coded to construct a so-called ActivityMap (AM). At least five images are required to create the AM, resulting in a visualization to distinguish between low [ Fig. 5(b), yellow], mean [ Fig. 5(b), orange], and high-frequency [high activity, Fig. 5(b), red] change regions. Low activity changes correspond to image pixels that changed once within the observed time span. Mean activity changes represent pixels that changed a maximum of three-times and high activity changes were detected at least four-times along the time series.
For each high activity pixel, the matching HAOs are identified by evaluating the time series CM information at the location of the specific HAO. Since an HAO consists of multiple high activity pixels, a nonredundant extraction of the objects has to be accomplished to ensure that features are calculated for one HAO only once. For the dataset used in this study, a total of 34,363 HAOs were found. A set of features are calculated for each HAO to provide a detailed description of the object. We focus on commonly used features such as radiometry, local statistics, texture, and morphology. 49 A relatively new technique given by the DAPs 50 was adapted to our problem by integrating the SAR-specific attributes sigma nought (Sig0, σ 0 [dB]) and coefficient of variation (Cov). DAPs are the derivatives of the morphological attribute profiles (APs), which represent an advancement of the morphological profiles (MPs), originally introduced at the beginning of this millennium. 51 Both MPs and APs denote approaches for the hierarchical morphological filtering of image content according to a predefined value sequence. In contrast to MPs, APs utilize the attribute filter principle basing on CCs as already described in association with the connected ASF in Sec. 3. Hence, the selection of an SE of fixed size and shape is not necessary, and a more flexible image exploitation can be performed, modeling the image content in a more comprehensive way. 50 DAPs provide information about the sensitivity of CCs against a certain attribute. 50 In Ref. 45, the DAP principle was used to emphasize their advantage for the change categorization task. Furthermore, in Ref. 52, it was shown that promising results can be obtained using DAPs for SAR amplitude image segmentation.
Beyond the σ 0 and Cov attribute, which are below given in more detail, additional values were calculated for the HAO-based DAP information extraction [e.g., moment of inertia (Refs. 53 and 54); compare Ref. 14]. For each attribute sequence, five values are calculated leading to 1 × 10-sized DAPs.

Sig0
The SAR-specific attribute Sig0 stands for the mean calibrated backscatter intensity measured for an HAO. To derive a meaningful value sequence, the Sig0 backscatter information for different surface materials or structures published in Ref. 55 is utilized. From this, it follows that a maximum value of þ10 decibels can be considered as suitable mean value for surfaces with a high backscatter value. As a minimum, a value of −22 decibels was chosen by evaluating for the backscatter properties of asphalt surfaces. It was also observed that this value is in approximate correspondence with the noise equivalent sigma zero of most current X-and C-band sensor systems. 55,56 With these values, the sequence is defined as φ

Cov
Here, the Cov is defined as the ratio between the standard deviation and the mean value on the amplitude image after a 5 × 5-sized moving average kernel was applied. Consequently, it can be interpreted as a measure of local heterogeneity, leading to relatively high values in manmade environments (e.g., urban structures). In Ref. 55, suitable minimum and maximum Cov values were already described. The resulting sequence is φ Cov For change categorization, the benefit of applying the DAP technique using the value sequences outlined above is demonstrated by the comparison of two HAOs located in different high activity areas within our SAR scene (Fig. 6). Here, all considered DAP attributes are visualized.
It can be seen that the signatures (bottom part of Fig. 6) allow a clear distinction between the two different HAOs. The signatures summarize the DAP values calculated for the HAOs with respect to the defined value sequences. A peak in a signature indicates the sensitivity of an HAO to a specific attribute. In other words, it shows when an HAO is deleted during the attribute filtering. Using all features with their individual value sequences, a feature matrix of size 34; 363 × 72 finally results for the Stuttgart airport dataset. Due to the high dimensionality of this feature matrix, a dimension reduction is performed in our approach to ease subsequent data analysis steps.

Dimension Reduction
Different algorithms for dimension reduction were tested. From this analysis, the well-known principal component analysis (PCA) was identified as optimal method. 57 PCA searches for a multidimensional linear transformation, in which the input data points (here: HAOs) are optimally separable from each other. This transformation is defined by the eigenvectors (principal components) and eigenvalues of the dataset, which are calculated from its normalized covariance matrix. The principal components correspond to the axes of the new coordinate system. The eigenvalues describe the amount of data variance represented by each of the new data dimensions. 58 Each of the resulting principal components is a linear combination of the original input variables (here: features). Consequently, the influence of the original features on the principal components can be directly evaluated by calculating the so-called loadings.
One aspect of the PCA-transformed dataset is the count of the principal components that are sufficient to describe the dataset without significant loss of information. One strategy for identifying required principal components is the evaluation of the cumulative variance of the participating eigenvectors. For this, the percentages of the total variance are summed up for each principal component until a defined threshold is reached. Common thresholds hold at least 60% of the total variance. 59 In this work, a threshold of 80% is used to ensure that no important information is lost. Moreover, with respect to our dataset, it was found that only the first two principal components are needed to represent the desired amount of information.

Categorization
The categorization concept is closely related to information clustering in the sense that it attempts to aggregate HAOs with respect to their feature characteristics. Hence, a categorization can be seen contrary to a classification, where a predefined system of classes has to exist. 60 Multiple algorithms were tested and the k-means clustering was found to be the optimal choice. 57 Furthermore, an unsupervised version of k-means as described in Ref. 52 was implemented to support a high level of automation. This version of k-means is based on the minimization of the mean quantization error (MQE) E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 3 0 2 where N stands for the number of instances (here: HAOs) and M for the number of cluster centroids C. With the function ξ, it is tested whether a specific HAO x is part of the current cluster. The distance function δ represents the calculation of the Euclidian distance between the HAO and a particular cluster centroid. The clustering is performed for different number of clusters k (e.g., k ¼ 1; 2; : : : ; 10) in ascending order. Afterward, the extracted MQE for each iteration is plotted against the number of clusters. A gradient analysis is performed on the plot and an optimal k is found where the difference between two adjacent gradients is at a maximum (elbow method 61 ). Due to the fact that k-means is based on a random initialization of cluster centroids, one would get different clustering results for each iteration. Hence, the clustering is performed for one specific number of clusters multiple times (e.g., 10 iterations). The overall categorization strategy used in this study is visualized in Fig. 7. At the beginning, the reduced feature matrix FM M×N (M: Count of HAOs; N: count of principal pomponents) is used as input for the k-means clustering. A MATLAB implementation of k-means is used with the settings replicates = 10 and emptyaction = drop. The first setting indicates that 10 iterations are conducted. The second setting avoids the generation of clusters with no HAOs. The clustering is executed 20 times while consecutively increasing the number of clusters leading to a 1 × 20-sized mean MQE vector. From this vector, the gradients are calculated and a k Ã is extracted, which represents the optimal number of clusters. With k Ã , 10 repeated k-means runs are initiated to ensure representative results. The output is the reduced, categorized feature matrix KM M×N .

Classification
As shown in Fig. 1, the classification step mainly consists of an FP and a backprojection (BP) step. As input for the FP, an initial predefined class catalog is jointly used with the categorized feature matrix KM M×N containing the cluster memberships for each HAO.

Initial Class Catalog
The selection of an initial class catalog consisting of high activity areas of interest represents the only interaction of the method with the user. Since the AM containing the high activity pixels is available to the user as guidance for catalog creation, even an unexperienced user should be able to define a class catalog without difficulty. Each initial class is marked as a polygon around a region of interest (e.g., parking lot, etc.). Figure 8 shows an example of the class catalog creation process. A user is presented with identified AM information [ Fig. 8(a)] that can be compared to underlying image data [e.g., Fig. 8 In this example, HAOs cluster at a location belonging to the land use class collecting point. Once labeled by the user, this class would include open spaces in urban environments used for the collection of waste, garden rubbish, building materials, or similar items.

FP
The FP delivers information about the existence or nonexistence of a specific class of the initial class catalog. Within the FP, the categorized HAOs are visualized in the image space, combined with the user-defined initial classes [ Fig. 9(a)]. For example, Fig. 9(a) shows categorized HAOs projected onto the SAR image information. The color assignment of the HAOs corresponds to their cluster membership as determined during the categorization step (e.g., blue HAOs belong to cluster #3). Also shown are the outlines of the class sites that were initially defined by the user (here, site for class collecting point). Based on this information, the relationship between the class catalog and the HAO categories can be analyzed.
Class-cluster-histograms are created that indicate the contributions of each individual HAO cluster to a class. For each class, one class-cluster-histogram results, containing the percentages of the HAO cluster memberships. To demonstrate this process, Fig. 9(b) shows the class-clusterhistogram for the class collecting point.
In this case, two clusters (3 and 4) dominate the class definition, leading to the assumption that this class actually exists in the data reality. In this study, a class is considered as existing when the threshold criterion E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 3 4 5 is fulfilled. This means that for an example of k Ã ¼ 4, a class is considered to be existing if it contains at most two dominating clusters. With i, ii, j, jj ∈ N, the dominance of clusters is defined as follows: • Exactly one dominant cluster c i exists for class C j , when it holds more than twice as much HAOs as the second-most relevant c ii . In that case, the number of dominant clusters CL equals one (case a). . . . the difference between these two sums is bigger than 20% (found by test runs); (case b.2).
In case of b., it is 1 < CL ≤ k Ã ∕2. The parameter CL, determined within the FP step, is an input variable for further processing. It represents the basis for the case analysis conducted as part of the BP.

BP
Besides the parameter CL, the categorized feature matrix KM M×N is used as input for the BP. As mentioned earlier, it was sufficient only to apply the first two principal components, resulting in a feature matrix of size KM M×2 . Hence, the extraction of empirical key indicators for each class can be performed in a 2-D feature space. Such indicators serve as the basis for an optimization of the initially selected class catalog. The aim is to arrive at a final catalog that matches optimally to the actual existing data reality. For the calculation of the key indicators, each initial class is first visualized in the 2-D feature space of clustered HAOs (FP step). Figure 10 exemplifies this process. Figure 10(a) shows the initial class boundary (magenta-colored solid box) displayed in the 2-D feature space. It can be seen that this particular class overlaps mostly with two of the four HAO categories.
Density buffers D i [symbolized by rectangular-shaped regions around the HAO cluster centers; Fig. 10(b)] are constructed by a region growing procedure, starting around their centroids. The optimal case exists when the involved clusters represented by their density buffers are completely lying within the boundaries of an analyzed class [D i ¼ 100; e.g., red cluster in Fig. 10(b)]. The partitioning degree Z indicates whether at least one cluster is cut by the class or not. An optimal configuration of Z exists when no involved cluster represented by its density buffer is cut by the class boundaries (Z i ¼ 0). In Fig. 11, the case analysis workflow is visualized. Subsequently, its individual steps are outlined.

Case 0 (class exists; CL ¼ 1)
This case is indicates that a class has already been optimally defined. The class covers the density buffer of exactly one cluster completely (D ¼ 100). From this, it follows that the partitioning degree Z is Z ¼ 0. In consequence, not further optimization is necessary.

Cases 1, 1a, and 1b (class exists; 1 < CL ≤ K Ã ∕2)
In this case, the class covers at maximum k Ã ∕2 clusters completely (D i ¼ 100). The differentiation of the cases 1a and 1b is realized by proving that a split of an initial class into two subclasses is meaningful. Such a split is warranted if a class has been defined too coarsely, containing other, actually independent classes. The test for a possible class split is accomplished by analyzing the mutual overlap of the density buffers of involved clusters in the 2-D feature space.

Cases 2, 2a, and 2b (class exists; 1 < CL ≤ k Ã ∕2)
If an initial class contains a maximum of k Ã ∕2 dominant clusters, but does not cover at least one of the involved density buffers completely (D i ≠ 100, with 1 ≤ i ≤ k Ã ∕2), an optimization is initiated. This optimization is performed by locally shifting the class in the feature space.
With this, it is tested whether a better configuration of the key indicators is achieved while retaining the original rectangular size and shape of the class. This process is schematically shown in Fig. 12. In this example, the local shift leads to a better configuration of the originally extracted key indicators.
After the optimization, it is further tested whether the class should be split (case 2a) or not (case 2b). The procedure is the same as for cases 1a and 1b.

Cases 3, 3a, and 3b (class does not exist; CL > k Ã ∕2)
Case 3 represents situations where an initially selected class contains more than a maximum of k Ã ∕2 dominant clusters. Hence, the class is considered to not exist in the data reality. Nevertheless, it is tested whether the class can be automatically redefined (case 3a). If a redefinition is not possible, the class is deleted (case 3b).
For a redefinition, the class has to fulfill the criterion of covering at least one of the involved density buffers completely (D i ¼ 100 with 1 ≤ i ≤ k Ã ∕2). If the class fulfills the criterion, a  region growing process is initiated, starting with a small box around the centroid of the cluster being completely covered. The region growing stops, when at least one of the key indicators becomes suboptimal. With this redefinition, the user is supported in learning the actual available information content of the dataset. The algorithm proposes an alternative class definition that better matches to the data reality. In general, it should be kept in mind that the redefined class is no longer connected to the real-world semantics of the original class. Current research is focused on the development of a procedure that supports the user in interpreting redefined classes.

Final Class Catalog and Class Assignment
As output of the case analysis, a final catalog is created, containing classes that optimally match to the data reality. To assign each HAO to these optimized classes, the class boundaries are projected into the 2-D feature space. Since final classes can overlap with each other in feature space, a fuzzy-based strategy is chosen to handle possible ambiguities of class membership.
The core of the fuzzy logic 62,63 approach is the calculation of degrees of membership of each HAO with respect to the optimized classes. In contrast to other classification techniques, this approach provides the capability to deal with datasets that have not been precisely formulated. 64 The concept for the fuzzy-based assignment of HAOs to final classes is shown in Fig. 13. In Fig. 13(a), the feature space extent of an exemplary class C j is presented with its centroid being indicated by the circle in the class center. The class area is divided into five so-called fuzzy buffers, whose widths are determined based on Gaussian-shaped membership functions. Choosing five buffer regions was found to provide best performance in test runs.
The membership functions are defined along both class axes and are centered on the class centroid. Hence, high degrees of membership are assigned to HAOs that lie near the class centroid while low degrees of membership result for HAOs located near the edges of the class. The assignment principle in the case of multiple overlapping classes is shown in Fig. 13(b). Blackand magenta-colored dots (marked by dotted and dashed arrow, respectively) denote the locations of two different HAOs in the 2-D feature space. Three partially overlapping classes C 1 , C 2 , and C 3 are visualized with their buffer regions. Focusing on the black-colored HAO, an assignment to class C 1 would occur, due to the higher degree of membership (yellow buffer: 11% to 25%) in relation to class C 2 (red buffer: 0% to 10%.). If an HAO is located in the same buffer zone of multiple classes (here: magenta-colored HAO in cyan buffer of classes C 1 and C 3 ), the Euclidian distances between the HAO and the involved class centroids are evaluated, and the HAO is assigned to the class with the smallest distance (here: class C 3 ).
The assignment process is additionally guided by a minimum degree of membership z, meaning that an HAO is assigned to a specific class only if its degree of membership lies above this threshold. If a membership level of z is not reached for any class, the HAO remains unclassified. Consequently, the higher the threshold z, the higher the potential rate of unclassified HAOs. The selection of a single z for all datasets is not recommended since it depends on the structure of the 2-D feature space and the neighborhood of clusters arranged within.

Experimental Results
In this section, we apply our approach to the SAR image time series described in Sec. 2. Concerning the change classification step (see Fig. 1), we focus on two classes that are described in more detail in the following paragraphs.

Initial Class Catalog
As basis for the classification task, an initially class catalog consisting of six classes is created (Fig. 14). In Fig. 15, the selection of the classes construction site and collecting point (Sec. 4.1) are visualized to exemplify the class selection process.
At the construction site [ Fig. 15(a)], a construction project for the transformation of a former parking garage [ Fig. 15(A.2)] into a new business complex [ Fig. 15(A.3)] is shown. The construction project was in progress during the acquisition dates of all images of the time series. Hence, the movement of construction machines, trucks, cranes, and debris create a high activity change environment in this area.

Reduction of Feature Matrix, Categorization, and Forward-Projection
After calculating a feature vector for each HAO leading to a high-dimensional feature matrix, the reduction strategy described in Sec. 3.2 was applied. Eleven principal components were selected for further processing, with the first two representing more than 40% of the whole original data variability [ Fig. 16(a)]. Through an evaluation of the loadings (Sec. 3.2), it was found that the DAP features (Sec. 3.1) provide important information for the distinction between change categories. An example of the importance of DAP features is shown in Fig. 16(b), where the loadings on the first principal component are presented. For the dataset presented here, the k-means clustering approach as described in Sec. 3.3 leads to four clusters, which are visualized in the 2-D feature space represented by the first two principal components [ Fig. 17(a)]. The boundaries of the six initial classes are displayed as colored rectangles [ Fig. 17(b)]. The two classes 3 (dashed green box) and 5 (dotted dark yellow box) are highlighted.
The FP step evaluates the percentages of HAO cluster memberships in the initial class regions (Sec. 4.2). Based on this information, the HAOs being part of the initial classes are identified and their cluster memberships are extracted. Focusing on the two classes construction site and collecting point, this leads to the class-cluster-histograms as shown in Fig. 18.
Considering the criterion given in Sec. 4.2, one can observe that both classes have a maximum of k Ã ∕2 dominant clusters, meaning that both classes exist. It follows that: CL ¼ 1 (class 3)  and CL ¼ 2 (class 5). Consequently, the cases 0 (class 3), 1, or 2 (class 5) of the workflow in Fig. 11 come into consideration for further processing.

Backprojection
For all initial classes, the key indicators as described in Sec. 4.3 are calculated and iteratively evaluated. The aim is to produce an optimized class description with respect to the indicators. In Fig. 19, the BP history for all classes is listed, providing information about the initial and final configuration of the indicators. The applied case is also listed, providing information on how each initial class was treated.
Focusing on the two highlighted classes construction site and collecting point, it can be observed that the cases 0 and 2b were applied. This means that the class construction site exists in the data as initially selected and does not require any further optimization (case 0). For the fifth class (collecting point), an optimization by locally shifting the class center in the 2-D feature space was necessary (case 2b).
It can be seen that this shifting led to a better class configuration. Moreover, the values of the parameter D indicate that the two dominant clusters 3 and 4 are better covered by the modified class. The partitioning degree Z has been also reduced to one. In Fig. 20, the optimization process is visualized in the 2-D feature space. It is observable that the optimized classes are better matched to the data reality described by the four change clusters. Since the final classes 1b, 2, and 6b cover the identical area in feature space, only class 1b is visualized.

Final Class Catalog and Class Assignment
The optimized version of the class catalog consists of seven classes (1a, 1b, 2, 3, 5, 6a, and 6b). In Fig. 21, these optimized classes are listed including their relationships to the initially selected classes. For HAO assignment, a minimum degree of membership of z ¼ 50% was chosen. In several test runs, it was found that this value produces best results for the analyzed dataset. In Fig. 22, the assignment results for the two highlighted classes are shown. Here, an industrial region is displayed consisting of three areas of interest: a container collecting point (A), a mixed area of collecting points and loading/unloading points (B), and a car parking lot (C).
From Fig. 22, two main conclusions can be drawn: First, the spatial consistency of HAO class memberships suggests that the achieved classification results are meaningful and match the local data semantics. For example, the class collecting point (Fig. 22, lower row) dominates the small container collecting point (A) indicating the robustness of the HAO assignment step. Second, the classes construction site and collecting point dominate area B (mix of loading/unloading points; Fig. 22, upper row (right) and lower row). Also, this result meets the expectations and can be regarded as exact and correct. The performance of the method is also visible in the third area of interest, a car parking lot for employees. In this area, neither the construction site nor the collecting point class is prevalent.

Validation
In the following, a validation of the presented method is accomplished by analyzing its robustness against different initial class catalogs. As the selection of the initial class catalog is done by the user, it is important to assess the robustness of the algorithm against variability in user input.
To analyze the impact of different initial class catalogs on the final classification result, we consider two constellations: in a first run, we simulate a small catalog consisting of only three classes, whose exemplary sites were chosen to have considerable spatial extent (left-hand side of Figs. 23 and 24). In a comparison run, we simulate a user that created a more diverse initial catalog by defining most spatially distinct sites (right-hand side of Figs. 23 and 24). Based on these data, we assess whether the developed algorithm can cope with this diversity of inputs by producing consistent final classification results.
In Fig. 23, the initial class catalogs for both scenarios are given. The difference between coarsely and finely selected class sites is shown in Fig. 24, where the sites for the classes 3 and 10 are visualized (class airport ramp in both catalogs). It can be seen that this class was  defined radically differently in the two processing runs, representing a much more extended area in the coarse catalog containing different HAOs. These are, for example, airplanes, baggage cars, loading/unloading points, mobile stairways, etc. In contrast, the more detailed catalog chose a more strict definition of the class airport ramp (Fig. 24, right).
Referring to the constraint given by the criterion CL ≤ k Ã ∕2 (Sec. 4.2), a coarser selection of initial classes has to lead to deletions or redefinition of initial classes. This assumed behavior is confirmed with the BP results for the coarse class catalog (Fig. 25).
It can be seen that case 3a was applied to the first two initial classes, resulting in a redefinition of the initial class. These class redefinitions were possible as at least one involved cluster was covered completely by a specific class. In contrast to these redefinitions, class 3 airport ramp has been retained and split into two new classes 3a (Fig. 27, A, lower row, dark blue dashed box) and 3b (Fig. 27, A, lower row, dark brown dotted box). This is caused by the fact that the initially  Much more and finer selected classes are processed with the second catalog [ Fig. 23(b)]. Contrary to the coarser defined class catalog, a high rate of retained classes is assumed. The selection of a smaller class region representing a very local image area causes a much lower mixture of involved clusters. Consequently, a dominance of CL ≤ k Ã ∕2 clusters is much more likely to exist, leading to a relative small rate of class deletions or redefinitions.  With Figs. 26 and 27(b), it can be seen that the algorithm provides results, which match to these assumptions. Only two initial classes (6 and 9) have been completely deleted. Focusing on the case analysis, case 2b dominates the BP processing results, meaning that these classes were optimized by local shift in feature space.
For one initial class (truck parking lot), case 3a is applied, which means that it is considered to not exist in the data reality. Since at least one involved density buffer is covered completely (Sec. 4.3.4), the class is redefined by the algorithm. Focusing on Fig. 27(b), lower row, it is observable that the redefined version of this class (dotted dark green box) covers the core of the cluster green (marked by dotted arrow a1).
In summary, the results show that the proposed BP procedure works robustly for different initial class catalogs. Dominant BP cases were deletions or redefinitions (coarser selected class catalog) and optimization by local shifts in feature space (finer selected class catalog). In both scenarios, the final classes are in good accordance with real-world semantics. For future work, an automatic procedure, which provides support in interpreting such final classes, will be considered.
Additional evaluation results are given in Ref.
14. There, for example, the method is applied on a complementary dataset, leading to similar and meaningful conclusions. Therefore, we are convinced that our method provides robust and precise results. The dotted arrows mark the two clusters green (a1) and red (a2), which are dominant for the initial class 3 (airport ramp). (b) lower row: The dotted arrow a1 marks the cluster green, from which the core is covered completely by the redefined class 3.

Discussion
In this section, it is discussed whether the proposed method meets the requirements for a new change analysis method, as outlined in Secs. 1 and 2. There, it was mentioned that such a method has to be as portable as possible. Moreover, in an ideal case, a classification-based analysis of the changes has to be performed without any reference or training data. A high degree of automation and usability has to be pursued.

Automation and Portability
To ensure sufficient portability, a high degree of automation was targeted for the development of our processing flow. The smaller the number of parameters to be set by the user, the better the compatibility to other datasets.
Starting with the first steps of the workflow (Fig. 1), the LR method combined with the connected ASF (ASF cc ) represents a fast, simple, and robust change detection method for SAR time series data. The usage of the ASF cc filter results in the detection of specific-sized changes. At the current state, an appropriate configuration of the ASF cc has to be done by the user. Nevertheless, it is conceivable to automatize this step by, for example, considering metadata information concerning the spatial resolution of the sensor system.
Focusing on the categorization of the detected changes, an unsupervised version of the k-means clustering algorithm was applied. With the help of the MQE calculation (Sec. 3.3), the optimal number of clusters for a specific dataset is extracted. The categorization itself is performed on the basis of a previously calculated feature matrix, being determined during the feature extraction step (Sec. 3.1). Among all extracted features, the morphological feature, given by the DAPs was found to be the most relevant in change categorization. For the DAP calculation, several value sequences have to be set according to the current available data. Consequently, the DAP feature has to be adapted when applying the method to different datasets. An automaton of this adaption might be achieved by including metadata information.
Concerning the change classification procedure, the two main steps FP and BP work fully automatic. The only input parameter to be set by the user is the initial class catalog (Sec. 4.1). Key indicators were defined to produce a final version of the initial catalog, optimally matching to the actual data reality (Sec. 4). Since the key indicators serve to evaluate the structure of the feature space with respect to the existing clusters, they are directly portable to other datasets. The predefined parameter CL in combination with its threshold criterion CL ≤ k Ã ∕2 to analyze the existence of a specific class, has undergone robustness tests that were published in Ref. 14. For the last step of the change classification, the user has to choose a value z, representing the minimum degree of membership of an HAO with respect to its final class (Sec. 4.4). Since it is not possible to determine a value for z that is suitable for all datasets, this parameter has to be chosen depending on the considered remote sensing scene.

Usability and Practicability
Due to the small number of input parameters, a high degree of automation has been achieved leading to a high usability of the algorithm, especially for inexperienced users. A sparse input dataset is used to demonstrate that relevant change classification results can be achieved without requiring additional reference or training data. Hence, it can be stated that the method is an easyto-use approach, fast, and easily applicable.
From this, it also follows that the method holds a high degree of practicability. In Sec. 1, it was determined that no method previously existed that was capable to provide detailed change analysis without the use of comprehensive reference or training data. It was found that many of the conventionally used methods required time-consuming preclassification of the scene based on a certain degree of expert knowledge. Our approach copes with these limitations. Moreover, it supports the user in learning the actual information content of the scene and in formulating a class structure that optimally matches to the data. In addition, we are convinced that the basic principle of our method can be applied to other land use classification issues not strictly connected to the change detection problem. This might be of interest for situations, where a classification should be performed while sufficient reference or training data are not available.

Limitations
Concerning the limitations, two main aspects should be mentioned. The first one relates to the interpretation of the final change classes. Since, in some cases, the final classes may differ from the initially selected scene semantics, it might be hard to draw conclusions about the land use they represent. For this, the user is recommended to utilize complementary data (e.g., optical space-or airborne imagery). Also the context classes that are calculated by the algorithm for each HAO may be of use here (see Ref. 14).
The second limitation of the approach is the requirement for a selection of the minimum degree of membership z. At the current state, this parameter has to be manually set by the user and may require a bit of trial-and-error testing. Although a generally suitable value for z does not exist, recommended ranges for z may be proposed by the method with the help of reference data. Such data might be land use maps representative for the acquisition dates of the SAR data. Unfortunately, suitable land use maps may not be available in all cases.

Conclusion
In this paper, an innovative approach for SAR change analysis is presented. In general, change analysis has been more or less strongly connected with the issue of land use classification. However, previous change analysis methods required a significant amount of reference information to result in meaningful classification results. Hence, it was the goal of this work to develop a method that can provide classification of detected changes without requiring comprehensive reference or training data, resulting in an approach that is easy to use and has high practical relevance.
The method can be divided into two major steps, the first of which is the detection and categorization of changes. In our study, the focus is put on manmade change objects located at different structures of urban environment. Since time series data are used, areas can be detected that change very frequently along time. Such high activity areas (e.g., parking lots) consist of so-called HAOs (e.g., vehicles), which represent the basis of our investigations. Nevertheless, the method is not strictly fixed to the analysis of HAOs. Other change objects might also be considered but are not addressed in this paper.
HAOs are detected by applying the LR method to pairs of SAR images. For ratio image denoising, ASFs are applied. The filtering with ASFs has the advantage that change objects of specific size ranges are retained, whereas changes beyond set size limits are rejected. Doing so, a high flexibility with focus on the change detection results is achieved.
After their detection, HAOs are utilized as input for a feature extraction step, a preprocessing step applied before change categorization. Commonly used image features are extracted, among which the so-called DAPs were found to be most useful in change categorization.
The classification of HAOs represents the second major aspect of the presented workflow. Since no additional reference or training database is considered, the categorized HAOs are regarded as data-driven reality. On the basis of this data reality, an initial class catalog is selected by the user. Within the scope of the so-called FP procedure, each selected class is evaluated by identifying the HAOs being part of it. The count of dominant clusters is used as input parameter for the BP step, which represents an iteratively optimization procedure. Several key indicators are evaluated for each class to produce a final class catalog, which matches best to the actual available data.
The final step of the classification procedure is given by the assignment of the HAOs to the final class catalog. For this, a fuzzy-based strategy is applied to handle possible overlapping final classes in feature space. Doing so, a degree of membership is calculated for each HAO, which is further proved against a threshold criterion.
To produce experimental results, the method was applied on a time series of 15 TerraSAR-X amplitude images. The scenery contains the Stuttgart (GER) airport area, surrounded by heterogeneous land use structures, leading to multiple change categories and to a high number of HAOs. With a comparison to real-world semantics, it was shown that the proposed method creates robust and precise classification results. The robustness of the approach was also shown by applying different initial class catalogs to simulate variability in user input. It was observed that the algorithm produced meaningful results for different initial class catalog configurations.
Regarding our method in general, it is also conceivable to adapt the basic principle on all land use classification tasks where no change information is available. The principle is not strictly fixed to the use of SAR images and might be applied on, e.g., optical data after an appropriate adjustment of the processing steps, for example, the configuration of the DAP attributes. It can be concluded that the innovative method presented for SAR change analysis provides a high degree of practical relevance and usability and produces robust and meaningful results.