1.

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,¹ 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 high-resolution Spotlight mode HS300 enables the acquisition of images with a geometric resolution of less than 1 m,⁴ 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–7 Focusing on land cover characteristics, the two items “conversion” and “modification” can be generally distinguished,⁸ 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.¹² 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.¹⁴ 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.¹³ In general, for this strategy, a classification of the input images has to be accomplished,¹⁸ which can be performed in a supervised or unsupervised manner.² 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.⁷ 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.⁷

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 pre-event 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 (${\sigma }^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. 2324.–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³⁴ 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 so-called 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.

2.

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.

Fig. 1

Concept of the proposed method for change analysis showing the two main steps: detection and categorization (cyan-colored) and classification (magenta-colored).

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 heterogeneous land use structure, consisting of urban and industrial areas [Figs. 2(b) and 2(c)], rural areas [Fig. 2(d)], and infrastructure [e.g., highways, Figs. 2(e) and 2(f)]. Hence, a relatively high rate of different change categories is expected. An additional dataset was considered for evaluation purposes. Results are discussed in Ref. 14.

Fig. 2

Test scene. (a) Time series mean SAR amplitude image; (b) mean amplitude image with WorldView-2 image. (A) Airport area, (B) urban structure, (C) industrial land use, (D) rural areas, (E) highway A8, and (F) federal highway B27.

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 (${\sigma }^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.³⁹

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)].

3.

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.

Fig. 3

Change detection flowchart. $A_{t1}$, $A_{t2}$: input amplitude image pair; LR: log-ratio calculation; $A_{mr}^{*}$: maximum ratio image; ${\mathrm{ASF}}_{\mathrm{CC}}$: connected alternating sequential filtering; $A_{mr,{\mathrm{ASF}}_{\mathrm{CC}}}^{*}$: denoised maximum ratio image; mean filtering: usage of a mean $3\times 3$ box filter; $A_{mr,{\mathrm{ASF}}_{\mathrm{CC}}}^{**}$: mean-filtered denoised maximum ratio image; Renyi’s entropy: application of Rényi’s unsupervised threshold calculation; CM: binary change map.

Fig. 4

Change detection workflow visualized by input and output examples. The subset shows a truck parking lot. (a) WV-2 image, (b) and (c) $A_{t1}$ and $A_{t2}$, (d) $A_{mr}^{*}$, (e) $A_{mr,{\mathrm{ASF}}_{\mathrm{CC}}}^{*}$, (f) $A_{mr,{\mathrm{ASF}}_{\mathrm{CC}}}^{**}$, (g) ChangeMap.

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

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 (${\mathrm{ASF}}_{\mathrm{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 ${\mathrm{ASF}}_{\mathrm{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.¹⁴ 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\mathrm{m}\times 1\mathrm{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⁴⁷ 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.⁴⁸ 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 (filter size $3\times 3$) version of the input image $A_{mr,{\mathrm{ASF}}_{\mathrm{cc}}}^{*}$ [Fig. 4(f)]. Afterward, the calculated threshold is applied on $A_{mr,{\mathrm{ASF}}_{\mathrm{cc}}}^{*}$ to produce a binary image (ChangeMap, CM) containing the changes at the foreground [Fig. 4(g)].

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.

3.1.

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.

Fig. 5

Time series change detection. (a) WV-2 image, (b) AM, and (c) high activity pixels. Color-coding: Low activity pixels (yellow), mean activity pixels (orange), and high activity pixels (red).

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.⁴⁹

A relatively new technique given by the DAPs⁵⁰ was adapted to our problem by integrating the SAR-specific attributes sigma nought (Sig0, ${\sigma }^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.⁵¹ 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.⁵⁰ DAPs provide information about the sensitivity of CCs against a certain attribute.⁵⁰ 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 ${\sigma }^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\times 10$-sized DAPs.

3.1.1.

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 ${\phi }_i^{Sig0}=[-22,-14,-6,210]$.

3.1.2.

Cov

Here, the Cov is defined as the ratio between the standard deviation and the mean value on the amplitude image after a $5\times 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 ${\phi }_i^{\mathrm{Cov}}=[0,0.5,1,1.5,2]$.

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.

Fig. 6

DAP signatures of two HAOs located in different high activity areas. (a) HAO at a construction site and (b) HAO at a car parking lot. The small cyan-colored object located in the center of the yellow-dotted boxes denotes the particular HAOs.

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\times 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.

3.3.

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.⁶⁰ Multiple algorithms were tested and the $k$-means clustering was found to be the optimal choice.⁵⁷ 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)

Eq. (2)

$$\mathrm{MQE}=N^{-1}·\sum _{i=1}^N\sum _{k=1}^M\xi (x_i\in C_k)·\delta (x_i-m_k),$$

where $N$ stands for the number of instances (here: HAOs) and $M$ for the number of cluster centroids $C$. With the function $\xi$, it is tested whether a specific HAO $x$ is part of the current cluster. The distance function $\delta$ 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;\dots ;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⁶¹).

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 ${\mathrm{FM}}_{M\times 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\times 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 $K{M}_{M\times N}$.

Fig. 7

Categorization flowchart to produce a categorized feature matrix $K{M}_{M\mathrm{x}N}$ by using an unsupervised configuration of $k$-means clustering. The ${\mathrm{MQE}}_i$ stands for the MQEs being used to extract an optimal number of clusters $k^{*}$.

4.

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 $K{M}_{M\times N}$ containing the cluster memberships for each HAO.

4.1.

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(b)].

Fig. 8

Exemplary selection of the initial class collecting point. (a) Selected polygon (green) in AM with high activity pixels (red), (b) matching WV-2 subset image, and (c) photograph from on-site inspection.

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.

4.2.

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.

Fig. 9

FP. (a) Initially selected class collecting point (green polygon) and (b) class-cluster histogram showing the percentages of HAO cluster memberships.

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-cluster-histogram 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

Eq. (3)

$${\mathrm{CL}}^{\tau }\le k^{*}/2$$

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\in \mathbb{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).

• • At most $k^{*}/2$ dominant clusters exist, when … (case b)

  • ○ … the sum of their HAO percentages is larger than the sum of the HAO percentages of all other clusters and … (case b.1)

  • ○ … the difference between these two sums is bigger than 20% (found by test runs); (case b.2).

In case of b., it is $1<\mathrm{CL}\le 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.

4.3.

BP

Besides the parameter $\mathrm{CL}$, the categorized feature matrix $K{M}_{M\times 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 $K{M}_{M\times 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.

Fig. 10

Schematically visualization of four clusters in 2-D feature space of the first two principal components PC1 and PC2, containing the BP history for an exemplary class (magenta-colored solid box). (a) Initial class covering parts of four clusters and corresponding HAOs (colored dots). (b) Additional visualized density buffers (dashed boxes in cluster-matching colors) around the cluster centroids (black dots). (c) Optimized final class.

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.

Fig. 11

Case analysis workflow. Bluish-colored circles denote the cases 0 (directly retaining) to 3b (directly deleting).

4.3.1.

Case 0 (class exists; $\mathsf{CL}=\mathsf{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.

4.3.2.

Cases 1, 1a, and 1b (class exists;$\mathsf{1}<\mathsf{CL}\le {\mathsf{K}}^{*}/\mathsf{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.

4.3.3.

Cases 2, 2a, and 2b (class exists;$\mathsf{1}<\mathsf{CL}\le {\mathsf{k}}^{*}/\mathsf{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\ne 100$, with $1\le i\le 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.

Fig. 12

Local shift of an initial class (magenta-colored solid box) in feature space. (a) Class before the optimization ($D_{\text{green}}=100$; $Z=1$). (b) Class after the optimization (density buffers of the involved red and green clusters are completely covered; $D_i=100$; $Z=0$). (c) Shifting area (hatched box around the class centroid) and positions of the shifted class (cyan-colored dots).

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.

4.3.4.

Cases 3, 3a, and 3b (class does not exist;$\mathsf{CL}>{\mathsf{k}}^{*}/\mathsf{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\le i\le 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.

4.4.

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.⁶⁴

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.

Fig. 13

Fuzzy-based assignment of HAOs. (a) Exemplary class $C_j$ surrounded by two Gaussian-shaped membership functions modeling the degree of membership in the directions of the first two principal components PC1 and PC2. (b) Example to show the principle for the HAO assignment, based on two HAOs [black- (marked by dotted arrow) and magenta- (marked by dashed arrow) colored dots] surrounded by the classes $C_1$, $C_2$, and $C_3$.

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). Black- and 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.

5.

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.

5.1.

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.

Fig. 14

Initial class catalog. The classes 3 and 5 being emphasized in this paper are marked by dashed red boxes.

Fig. 15

Initially selection of (a) class 3 (construction site) and (b) class 5 (collecting point). (a.1), (b.1) Class polygon in AM with reddish-colored high activity pixels; (a.2), (b.2) class area in WV-2 image subset; (a.3), (b.3) photographs of the class regions taken during on-site inspection (April 8, 2015).

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.

5.2.

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.

Fig. 16

(a) Plot of eigenvectors (principal components) against their percentage of variability. The vertical green dotted line represents 80% of total variability. (b) Loadings for the several features constructing the first principal component PC1. The horizontal black dotted line stands for 75% of the highest loading.

Fig. 17

(a) HAOs in 2-D feature space arranged to four clusters 1 to 4. (b) Initial classes visualized as colored boxes in 2-D feature space with categorized HAOs (colored dots).

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.

Fig. 18

Class-cluster histograms as results of the FP for the classes (a) construction site and (b) collecting point, showing one dominant cluster (yellow) for class 3 and two dominant clusters (blue and yellow) for class 5.

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: $\mathrm{CL}=1$ (class 3) and $\mathrm{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.

5.3.

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.

Fig. 19

BP history showing the state of key indicators before and after the BP. The right column contains the several applied cases.

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.

Fig. 20

Alternatively visualization of the BP history. (a) Initially and (b) optimized class catalog in 2-D feature space. The color-coding refers to Fig. 17.

5.4.

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).

Fig. 21

Final class catalog in relation to the initially selected classes. The arrows denote classes, which cover an identical area in feature space.

Fig. 22

Results of HAO assignment utilizing $z=50$. (a) WV-2 image subset, (b) assignment result for final class 3 (construction site) as overlay on time series mean amplitude image. (c) Assignment result for final class 5 (collecting point). (A) Container collecting point, (B) mixed area of collecting points and loading/unloading points, (C) car parking lot. The color-coding of the high activity pixels refers to Fig. 20 or Fig. 21.

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.

6.

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.

Fig. 23

Initial class catalogs containing (a) three coarsely defined classes and (b) 10 finely defined classes.

Fig. 24

(a) Coarsely and (b) finely selected initial class regions for class airport ramp. 1: Class polygons in AM with reddish-colored high activity pixels; 2: class areas in WV-2 image subsets.

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 $\mathrm{CL}\le 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).

Fig. 25

Final class catalog in relation to the initially selected classes (coarsely defined). The arrows denote classes, which cover an identical area in feature space.

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 selected class region [Fig. 24(A.1)] contains two dominant clusters [green and red, Fig. 27(a), marked by dotted arrows a1 and a2] that are clearly separable from each other (compare Sec. 4.3).

Fig. 27

Visualization of the BP history in 2-D feature space of both validation scenarios. Upper row: (a) coarsely and (b) finely defined initial class catalogs; lower row: final class catalogs.a) 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.

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 $\mathrm{CL}\le 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.

Fig. 26

Final class catalog in relation to the initially selected classes (finely defined). The arrows denote classes, which cover an identical area 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.

7.

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.

8.

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.