Region-based automatic building and forest change detection on Cartosat-1 stereo imagery

ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

Contents lists available at SciVerse ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

Region-based automatic building and forest change detection on Cartosat-1 stereo imagery J. Tian a,⇑, P. Reinartz a,b, P. d’Angelo a, M. Ehlers b a b

German Aerospace Center (DLR), Remote Sensing Technology Institute (IMF), 82234 Wessling, Germany University of Osnabrueck, Institute for Geoinformatics and Remote Sensing (IGF), 49076 Osnabrück, Germany

a r t i c l e

i n f o

Article history: Received 4 July 2012 Received in revised form 15 February 2013 Accepted 24 February 2013 Available online 1 April 2013 Keywords: Stereo imagery Digital Surface Model (DSM) Change detection Forest change Industrial area change

a b s t r a c t In this paper a novel region-based method is proposed for change detection using space borne panchromatic Cartosat-1 stereo imagery. In the first step, Digital Surface Models (DSMs) from two dates are generated by semi-global matching. The geometric lateral resolution of the DSMs is 5 m  5 m and the height accuracy is in the range of approximately 3 m (RMSE). In the second step, mean-shift segmentation is applied on the orthorectified images of two dates to obtain initial regions. A region intersection following a merging strategy is proposed to get minimum change regions and multi-level change vectors are extracted for these regions. Finally change detection is achieved by combining these features with weighted change vector analysis. The result evaluations demonstrate that the applied DSM generation method is well suited for Cartosat-1 imagery, and the extracted height values can largely improve the change detection accuracy, moreover it is shown that the proposed change detection method can be used robustly for both forest and industrial areas. Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved.

1. Introduction Cartosat-1 (also called IRS-P5) was launched by the Indian Space Research Organisation (ISRO) in May 2005. The camera system on board the satellite acquires at the same time two 2.5 m spatial resolution stereo panchromatic images using a forward view camera with + 26° viewing angle and a backward view with 5° viewing angle (Kumar et al., 2006). The main purpose of this system is to acquire along-track stereo data for generating DSMs in large areas. It acquires images from the entire earth during a 126 day revisit cycle. The swath width is approximately 30 km (http://www.nrsc.gov.in). These characteristics of Cartosat-1 allow the monitoring of land cover changes for large areas with a relatively high repetition rate. However, Cartosat-1 has not yet been widely used for change detection until now due to the lack of multi-spectral channels. Using the stereo imagery together with Rational Polynomial Coefficients (RPCs) from the provider, Digital Surface Models (DSMs) can be generated using photogrammetric techniques. The generated DSMs provide height information for land cover analysis. DSMs from Cartosat-1 with 25 m  25 m lateral resolution were generated by Rao et al. (2007) for updating 1:25,000 or 1:10,000 scale topographic maps, in this research the RPCs were refined by using manually collected Ground Control Points (GCPs). ⇑ Corresponding author. Tel.: +49 8153 283588; fax: +49 8153 281444. E-mail address: [email protected] (J. Tian).

Martha et al. (2010) analysed landslide volumetric changes based on 10 m  10 m resolution DSMs, which were generated with SAT-PP photogrammetric software from ETH Zürich. With a newly developed dense matching method (Hirschmüller, 2008), DSMs with 5 m  5 m resolution can be obtained (d’Angelo et al., 2008). Therefore change detection for small forest patches and industrial buildings becomes possible. Numerous change detection methods using remote sensing with different kinds of images have been proposed by many researchers (Coppin and Bauer, 1996; Mas, 1999; Lu et al., 2004; Berberoglu and Akin, 2009; Blaschke, 2010; Klonus et al., 2011; Chen et al., 2012). The investigations can mainly be divided into pixel-based and region-based methods. In the first case, the change features from two images are compared for each pixel independently. In the second case, the images are segmented into disjoint and homogeneous regions, and then change features are extracted and compared for these regions (objects, segments). Literature research shows that the region-based remote sensing image analysis is attracting more interest approximately since the year 2000 (Blaschke, 2010). Many researches have been performed in comparing region and pixel-based change detection methods for remote sensing data. Walter (2004) used region-based classification results from two dates to generate a land cover change map. The original regions were obtained from existing GIS data, different land cover classes were better separated with region-based features. Desclée et al. (2006) used the region-based approach for forest change detection, the region-based change detection

0924-2716/$ - see front matter Ó 2013 International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS) Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.isprsjprs.2013.02.017

227

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

method exhibited a much higher Kappa Index of Agreement (KIA) (KIA = 60%) than the pixel-based method (KIA = 49%) using the same features from multi-spectral satellite images. Im et al. (2008) compared region and pixel-based change detection methods, the results showed that the region-based change detection method could also reach a higher KIA (about 90%) than pixel-based methods (KIA = 80–85%). Duveiller et al. (2008) applied the forest/ non forest classification method from Desclée et al. (2006) in region-based deforestation detection from 571 image pairs, where an initial image selection was performed to exclude all bad quality images and image pairs with no forest change. Aguirre-Gutiérrez et al. (2012) compared pixel-based, region-based methods and their combination in land cover classification. A higher accuracy was achieved based on a combined classification method along with a region-based change detection analysis. A brief look at the above comparisons between the pixel-based and region-based change detection reveals that the region-based method performs generally better. However, all of these papers were focusing on using multi-spectral information. Due to the lack of multispectral information not many tests have been performed on using single channel Cartosat-1 stereo data and no literature has been found for combined evaluation for automatic change detection including the derived DSMs. The research using Cartosat-1 images for automatic and semi-automatic change detection is still limited to large scale topographic monitoring or manual interpretation (Martha et al., 2010; Prabaharan et al., 2010; Kamini et al., 2006). Kamini et al. (2006) chose onscreen digitization to obtain land cover maps from Cartosat-1 data to monitor the land use in frequently flooded areas. Kumar et al. (2007) merged 2.5 m resolution Cartosat-1 data to 5.8 m multi-spectral Resourcesat-1 images in a visual interpretation procedure. Prabaharan et al. (2010) also used visual image interpretation to get land use and land cover (LULC) classification maps in the year 2008 from Cartosat-1 images, and compared them with the land cover maps from images of other sensors. In this paper, we develop a change detection method using a combination of height changes from DSMs, optimized region boundaries and spectral change from panchromatic images. Firstly, DSMs with 5 m  5 m resolution are generated with semi-global matching (SGM) from Cartosat-1 stereo imagery as described in Section 2. After orthorectification of the near-nadir images the initial regions are obtained from the segmentation results of these orthorectified images of two dates. The change analysis regions are generated through a combination of these segmentation results. A rule based region merge approach is provided to reach a better consistency between the homogeneous regions of the two dates. Several height and spectral features are then extracted at region level. Additionally to the above mentioned methods, DSMs of the two dates are used in a newly developed combined change detection procedure as described in Section 3. This approach is tested for two different tasks. The first task is forest change detection, performed in a test area in Bavaria, Germany, and the other task is building change detection in an industrial area, where we choose Istanbul, Turkey as our test area since several changes can be monitored.

2. Cartosat-1 stereo imagery and DSM generation

The RPC provided with the Cartosat-1 Orthokit, exhibit a much lower absolute accuracy than the ground resolution of Cartosat-1 data (2.5 m). Correcting the RPC for high quality geolocation needs sub-pixel accurate Ground Control Points (GCPs). GCPs of this quality are usually derived from ground survey or high resolution orthorectified images and DSMs, which are not available in many cases. In this study, we use the Shuttle Radar Topography Mission (SRTM) DSM as absolute geolocation reference instead. SRTM is publicly available with a resolution of 3 arc seconds. The absolute horizontal (CE90) and vertical (LE90) accuracy of SRTM was estimated to be below 10 m for Europe (Rodriguez et al., 2005). Since this reference is used for the stereo data sets of both dates, the co-registration accuracy of the Cartosat-1 DSMs and orthoimages is between 2 m and 5 m (d’Angelo et al., 2008). 2.2. Semi-global matching (SGM) procedure First, high quality tie points between the stereo pairs are established using pyramidal local least squares matching. These are used to generate a quasi epipolar stereo pair by aligning the columns of the two stereo images. SGM is performed based on the generated epipolar stereo pairs. This matching procedure avoids using matching windows, and is thus able to reconstruct also sharp edges. Instead of strong local assumptions on the local surface shape, the matching step is cast into an energy minimization problem. It performs a semi-global optimization by aggregating costs from 16 directions, which finds an approximate solution to the global energy function E:

EðDÞ¼

X X X fCðp;Dp Þþ P1 T½jDp Dq j¼1þ P2 T½jDp Dq j>1g p

q2N p

ð1Þ

q2Np

This energy function (Eq. (1)) consists of a data term, measuring the similarity of possibly corresponding pixels in the two images, and a regularisation term, which favours similar disparities for neighbouring pixels, but also allows large jumps at discontinuities. Function C is the pixel wise matching cost between the image and the disparities D (for points p and q). A combined use of the mutual information and the census cost function are adopted as described by d’Angelo (2010). The second and third terms of energy function E are regularisation terms that penalize disparity changes in the neighbourhood NP for each pixel p. T is a conditional function, it equals to 1 if the argument is true, and to 0 otherwise. When the disparity changes by one pixel, a constant penalty P1 is added to the energy function. When a higher disparity changes occurs, a larger constant penalty P2 is added, but P2 is not a very high value and therefore also abrupt changes are allowed. P2 should always be bigger than P1, for processing of Cartosat-1 imagery we choose P1 = 300 and P2 = 800, with data costs C normalized to be between 0 and 1024. Finally, the disparity map is recomputed by choosing the disparities with the smallest E. Sub-pixel disparities are computed by fitting a parabola to the aggregated cost values next to the disparity with minimum aggregated cost. To detect blunders and occlusions, SGM is performed from the first to the second and from the second to the first image, and only consistent disparities are kept. Additionally, small isolated disparity regions are rejected as outliers as described in d’Angelo (2010). The core SGM algorithm is described in Hirschmüller (2008).

2.1. Input data 2.3. Image orientation and DSM generation Cartosat-1 stereo pairs in Orthokit format are used as input data for the DSM generation. The described process utilizes the Rational Polynomial Coefficients (RPCs) sensor model, which is used to transform the three-dimensional object-space coordinates into two-dimensional image-space coordinates (Grodecki et al., 2004).

A preliminary point cloud consisting of approximately 1 million points is generated from a thinned SGM disparity map by spatial intersection. This point cloud is aligned with the SRTM DSM using a 3D affine transformation. This aligned point cloud and the known

228

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

image locations are used in the image resection, yielding affine RPC correction parameters. The complete process is described in (d’Angelo et al., 2008). The DSM is generated by spatial intersection of the complete disparity map using the previously estimated RPC correction parameters and subsequent reprojection into a local coordinate system. 2.4. DSM refinement Due to matching errors in regions with weak contrast or occlusions, the generated DSMs will possibly contain a small amount of outliers. To eliminate large outliers, in this paper a check against the SRTM is performed. All points whose height deviates more than 200 m from the unfilled SRTM DSM are removed. Removing outliers will produce holes in the DSMs. Moreover, unmatched pixels also appear as holes. These holes are filled with SRTM data using the delta surface fill method by Grohman et al. (2006). The delta surface fill algorithm will effectively interpolate small holes and seamlessly fill large holes, e.g. created by clouds or water areas, with SRTM data. Remaining holes in areas where both generated DSM and the SRTM DSM do not provide data are filled by inverse distance weighted interpolation. 3. Change detection A two-step region-based change detection procedure is proposed according to the characteristics of Cartosat-1 images in this paper (as shown in Fig. 1). In the first step, after co-registration of the two DSMs and orthorectified images, segmentation on coregistered Cartosat-1 images is performed to obtain the initial regions. Then the regions from two dates (date1 and date2) are combined to get an initial segmentation map. To cover the over segmentation produced from the region combination, a similarity measure is proposed to reach a reasonable segmentation level. Through this similarity measure small segments are merged with their neighbours if appropriate. In the second step, multi-level change information from images and height change information from DSMs and images are combined with weighted change vector

analysis. A common resolution of all image data and DSMs is required for a combined data usage, thus in this paper the orthorectified Cartosat-1 images are resampled to 5 m.

3.1. Region preparation An additional 3D co-registration is needed to remove any remaining shifts in all three dimensions between the two generated DSMs (Tian et al., 2011). Since in the DSM generation procedure SRTM has been used for RPC correction in the DSM generation for both datasets, the shifts between these data sets are very small. In order to improve the computation efficiently, we use a simple version of least square matching, for which only linear shifting is considered, the shift distances in three dimensions (Xs, Ys, Zs) are estimated via iterative 3D shifts adjustment based on the minimization of the following equation: m X n X ðZ ref ðx; yÞ  ðZ 2nd ðx þ X s ; Y þ Y s Þ þ Z s ÞÞ2

ð2Þ

x¼1 y¼1

in which m and n are the column and rows of the DSMs, (x, y) are the plane coordinates, and Zref and Z2nd represent the height values from the two DSMs respectively, where we have to choose one DSM as reference DSM (ref) in this step. The shift-values Xs, Ys and Zs are iteratively adjusted to obtain a final shift result. Additionally an outlier rate is used to control the iterative procedure. The outlier rate can be adjusted according to the DSM quality. For instance, a larger outlier rate should be chosen when the two DSMs have higher spatial resolution or time distance. The pixels classified as outliers will not be included in the shift distance calculation procedure. After that, the horizontal shifts Xs and Ys are also applied to the orthorectified Cartosat-1 images of the same date. After co-registration, homogeneous regions are generated by applying segmentation separately to the images of both dates. Here, the objective of segmentation is to produce small units that have different spectral characteristics from the areas nearby. Many segmentation methods have been introduced in computer vision, like watershed, level-set, mean shift (Comaniciu and Meer, 2002)

Fig. 1. Flow chart of the proposed method.

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

229

Fig. 2. Region combination procedure. (a) regions from date1; (b) regions from date2; and (c) combination result.

and several more. However, automatic methods can hardly reach a proper level of segmentation for all land cover classes of interest, well-separated from each other, since different kinds of objects require different segmentation levels. In this paper, an over-segmentation is produced with mean shift from EDISON library (Comaniciu and Meer, 2002), to obtain firstly small units for further processing. This segmentation method could also be replaced by another segmentation method which can reach an over-segmentation. After segmentation two region maps are obtained from the two datasets. An intersection of the two region maps is performed to get all possible change regions. Fig. 2 shows the combination procedure. Fig. 2a shows the regions resulting from segmentation of date1 image (before change). Fig. 2b shows the regions from date2 (after change). The combination result is shown in Fig. 2c. In a region combination procedure (Bruzzone and Prieto, 2000; Bovolo, 2009), many small regions composed of only one or two pixels are produced, primarily at the edges of the segmentation results of the two dates (as shown in Fig. 3). Due to small remaining mis-registrations or radiometric differences (e.g. different shadow cast) caused by seasonal or sun angel changes, region edges of segment maps from two dates cannot be matched perfectly. Moreover, DSMs usually exhibit lower quality (height values are less accurate) at segment/object edges compared to the centre of segments/objects. Therefore, these small regions will introduce false alarm to the final results, and a proper region merging is required before performing further steps. A detailed requirement of initial regions for change detection was proposed by Bruzzone and Prieto (2000). However, how to reach this acquirement is not mentioned in that paper. In the following section a region merging procedure is proposed to solve this problem. 3.1.1. Greedy strategy for small region merging Each small region Gs that needs to be merged is surrounded by several regions, which are hereafter referred to as candidate regions Gi (i = 1, 2, . . ., k, k is the total number of the neighbours). The objective of this step is merging small regions into one of the candidate regions. To reach this the similarity between these small regions and these candidate regions is measured. Following the greedy strategy (Cormen et al., 2001) one candidate with the minimum heterogeneity increase after merging it with the small region is selected. An energy minimization model composed of four energy functions is proposed for this step. The similarities between the regions are measured based on the image intensity difference map D and height difference map H.

D ¼ jf1  f2 j

ð3Þ

H ¼ jh1  h2 j

ð4Þ

Fig. 3. Small regions at the edges resulting from region interaction.

f1 and f2 represent the pixel values of the panchromatic images from date1 and date2. And h1 and h2 are the two DSMs from date1 and date2. Merging Gs into one Gi will produce Gi0 . We define the following energy functions to measure the distance between Gi and Gi0 . The first function (Eq. (5)) is the area weighted mean value distance based on D.

E1 ¼ x  jlðGi0 Þ  lðGi Þj

ð5Þ

where l(Gi) is the average pixel value of region i, the pixel values are extracted from the image intensity difference map D. A weight x is multiplied to the energy term because the larger Gi is, the less likely it is to change after merging with Gs. In this paper, we use the size of Gi (total pixel number) as the weight x. We choose the second energy function for our model as the intensity range of each region.

E2 ¼ jMaximinðGi0 Þ  MaxminðGi Þj

ð6Þ

where Maximin(Gi) is the subtraction of maximum value and minimum value of region Gi. Also E2 is computed based on the D image. We compute similar energy functions based on H (shown in Eq. (4)), denoted as E3 and E4. E3 is the weighed mean value using average height value and E4 is the height range of each region. Next we rank the candidate regions four times according to the four similarity energy function separately, the smaller the value of the energy function, the higher rank O it gets. Oij represent the rank order of Gi according to energy feature j. A ranking matrix X is generated recording the ranks of the candidates.

230

8 1 9 > O1 O12 O13 O14 > > > > > > > < 2 = O1 O22 O23 O24 X¼ > >   > > > > > > : k ; O1 Ok2 Ok3 Ok4

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

ð7Þ

A ranking vector fOi1 ; Oi2 ; Oi3 ; Oi4 g is generated for each candidate i (i = 1, . . ., k). The candidate with the highest sum of ranking P Minð 4j¼1 Oij Þ is chosen as the object to ‘merge to’ (since a high rank results in a low value). To speed up the merging procedure, all of the neighbourhood candidates should have more pixels than the defined size of the small regions. For instance, when a region with less than 10 pixels has to be merged with one of its neighbours, only the neighbourhood regions larger than 10 pixels size are considered as candidates. 3.2. Feature extraction Pixel-based change detection operates on single pixels, while region/object based change detection methods operate on a group of pixels belonging to one region/object that are grouped together, change features should be thus extracted for each region. 3.2.1. Change information from DSMs Height change features are extracted by comparing the DSM from two different dates. Due to the potential outliers in the DSMs, the integral mean value will also be influenced if we use the mean value of each region directly. A histogram of H (Eq. (4)) inside one homogeneous region without any change is shown in Fig. 4. As can be seen, most of the reliable height values are clustered in the centre part of the histogram. The very high values (more than 10 m in Fig. 4) or very small values (less than 15 m in Fig. 4) have a higher probability to be outliers. Therefore, in this paper, only the middle part (96%) of the height histogram is analysed. The 2% highest values and the 2% lowest values are omitted. 3.2.2. Weighted multi-level change vector analysis Many change features were reviewed by Radke et al. (2005). The mean value difference of the grey values, which means the

intensity image difference is still an essential feature to change detection and easy to implement. The first step of this approach is to generate an absolute-value difference image map D (Eq. (3)) as we have already used in the region merging procedure. Change vector analysis (CVA) (Malila, 1980; Johnson and Kasischke, 1998) has been used for land cover change detection; it accepts any feature vectors as input. The multi-level feature sets (Celik, 2009) which are obtained by considering a series of neighbourhoods {S = 1; 3; 5; 7; 9} surrounding the corresponding pixels is given by

u¼

ss 1 X Dp s  s p¼1

ð8Þ

This can be also considered as the difference map at different resolution level. With this multi-level representation of changes, the neighbourhood of the corresponding pixels are included to enhance the behaviour of changes in various land cover classes. Hereafter, the region-based change vector (v) is described by using the mean value of the change features from DSMs and images (v ri is the ith change vector for region r; r = 1, 2, . . ., nr; nr is the total number of regions.). Since the Euclidean distance is adopted for the CVA in this paper, a standardization of the change vector is preferable when the change features in the change vectors are in the same range. All of the change vectors are normalized by dividing with their standard deviation (r(mi) vi is the ith change vectors of all regions; i = 1, 2, . . ., l; l is the total number of region-based change vectors). Moreover, weights qi are used to balance the influence of the change vector from the images and from the DSMs. In this paper, we use five region-based change vectors from the image difference map (multi-level features) and only one region-based change vector from DSM (see also Section 3.2.1). Weight q = 0.1 is used for the change vectors from images, and q = 0.5 for the change vector from DSM. This weight can also be manually adjusted, if the DSMs exhibit obviously lower or higher quality. The formulated weighted multilevel CVA for each region (r) is

Fig. 4. Height distributions inside one region.

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  r 2  r 2  r 2 vl v1 v2 r kDV k ¼ q1 þ q2 þ ::: þ ql rðv 1 Þ rðv 2 Þ rð v l Þ

ð9Þ

Here a larger value of kDV r k indicates a higher probability of change. 3.3. Accuracy assessment In order to evaluate the proposed change detection method, the reference data for the comparison were extracted manually by using both the images and the generated DSMs. The following measures are involved for computing the agreement of the change detection result and the reference change map. (1) True positive (TP): the number of changed pixels correctly detected as changed. (2) True Negative (TN): the number of unchanged pixels correctly detected as unchanged. (3) False positive (FP): the number of unchanged pixels incorrectly detected as changed. (4) False negative (FN): the number of changed pixels incorrectly detected as unchanged. (5) Overall accuracy: OA ¼ TPþTN  100%. N (6) Kappa Index of Agreement: KIA ¼ PrðaÞPrðeÞ . 1PrðeÞ

where PrðaÞ ¼ OA; PrðeÞ ¼

ðTP þ FPÞ  ðTP þ FNÞ þ ðFN þ TNÞ  ðFP þ TNÞ NN

.where N is the total number of pixels. These measures are only used to compare two binary masks. Thus, in addition to these measures, Receiver Operating Characteristic (ROC) curves (Hand and Till, 2001), which can be used to describe the relationship of true detection rate and false alarm rate, are generated to explore the quality of features and change maps. The area under the ROC Curve (AUC) is used for the ROC curve comparison, a large value of AUC indicates a higher agreement with the manual interpreted reference data. 4. Experiments and results 4.1. Study area and data Six study sites are chosen for the experiments. Three of them are located in forest areas in Bavaria, Germany, and the remaining three study sites are located in an industrial area in Istanbul, Turkey. The Cartosat-1 datasets of Bavaria were acquired in the years 2008 and 2009 respectively. Since they were acquired with a time difference of only one year, which is not enough for a considerable increase in tree height, we are mostly interested in deforestation in this area. The Cartosat-1 images in Istanbul were acquired with a time difference of 3 years, the first was captured on August 14, 2008 and the second was captures on May 08, 2011. In this industrial area our focus lies on newly constructed buildings. The original panchromatic images and the change reference data for three areas in each case are shown in Figs. 5 and 6. Except the second test area for forest area, which covers 1.5 km  1.5 km to contain more changes, the other five test areas cover an area of 1 km  1 km. 4.2. Results and evaluation 4.2.1. DSM quality assessment A quality assessment based on pixel level and region level was performed by comparing the generated DSM with first pulse LiDAR data, which exhibit a point density of 0.5 points/m2. The quality

231

assessment was carried out on three typical land cover classes: fields with low vegetation, forest area and residential area. Fig. 7 illustrates a visual comparison of these classes. It can be seen, that after co-registration Cartosat-DSMs share the same terrain trend as the LiDAR-DSMs. Forest areas and large houses are clearly visible in the Cartosat DSMs. However, the Cartosat-DSMs show fewer details in small structures than the LiDAR-DSMs, and the object boundaries are not very sharp. These restrictions are mainly caused by the resolution of the Cartosat data. This is the reason why we do not investigate residential areas with small buildings in this paper. To investigate the quantitative quality between Cartosat-DSM and LiDAR-DSM, four groups of pixels are chosen and compared in this test. The dashed lines of the profiles shown in Fig. 8 stand for height values from the Cartosat-DSM, while the solid lines correspond to heights extracted from the LiDAR-DSM. It can be noted that in field areas ‘‘noise’’ values of about 2–3 m exist in the DSM generated with the stereo imagery, which could be expected as it is within the range of the pixel size of Cartosat-1 data. A more detailed statistical comparison of the LiDAR-DSM and Cartosat-DSM is shown in Table 1. Although the existence of outliers leads to a relatively high maximum deviation in forest and residential areas, the medium and mean height values are around 2–3 m, which is less than the horizontal DSM resolution. Only in the case of forest the standard deviation is higher than the DSM pixel size, but here also the natural variation of height is quite large and LiDAR height data cannot be regarded as perfectly correct. As can be seen from Table 1, DSMs are less sensitive to outliers when using the region-based method. And the proposed region merge procedure can also largely decrease the calculation task. Fig. 9 shows the original intersected region map with 926 segments, including 253 regions having only one pixel inside. Fig. 10 shows the region size distribution of the original region map. From this distribution it becomes obvious that a large amount of regions have very small numbers of pixels, which are often generated from sub-pixel co-registration problems. After merging regions smaller than 50 pixels, only 1/5 of the regions remain (Fig. 9b and c). The minimum region size definition should be adjusted according to the size of the objects of interest. These small regions are generally assumed as noise from misregistration or might be due to DSM generation errors. Hence it is necessary to test the performance with and without the proposed region merge approach. For this, both the calculated speed and KIA are compared. As can be seen from Table 2, after the region merge procedure, only 25–30% of the amount of regions are left but the KIA is increasing. Moreover, the decline of the region numbers can also speed up the calculation procedure. 4.2.2. Change feature assessment In order to evaluate the proposed weighted change vector analysis, we compare it with other methods of change detection using different features, for instance, by using only image intensity, multi-level features, height values as well as normal change vector analysis. The same experiments are carried out for industrial and forest areas. The improvement of the multi-level features and weighed CVA is proven through ROC curve comparison, and is illustrated in Fig. 11. The AUC is used to evaluate the performance numerically. Table 3 shows the relative contribution of the input features. In general terms, this table reveals the importance of height features. It can also be seen that adding multi-level image features and height features in both forest and industrial areas leads to a noticeable increase in accuracy. The multi-level features also perform much better than single pixel values. This can be explained as the multi-level models are using the information in the neighbourhood of each pixel, and the combination of many window sizes leads to robust texture information.

232

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

Fig. 5. Cartosat-1 data in forest area. Figs. a1–a3 are the panchromatic images from year 2008; Figs. b1–b3 are the panchromatic images from year 2009; Figs. c1–c3 are the change reference data for each dataset.

4.2.3. Change detection result and assessment Change maps in forest and industrial areas are generated and displayed in Figs. 12 and 13 respectively. In achieving an over-segmentation based on mean-shift, we use bandwidths (hs, hr) = (7, 6.5) in the industrial area and (hs, hr) = (5, 5.5) in the forest area. For better view of the generated change mask, it is stretched to values from 0 to 1 using a sigmoid curve with the method described in Tian et al. (2013). Moreover, this stretch method can extend the range of the middle part, which helps in threshold value selection. The stretched change map can be considered as a change probability map, the grey value for each pixel represents its change probability. In the first step, we evaluate the change probability map for better visual comparison. In the second step, we investigate the agreement between the reference change map and the generated change map by using OA and KIA. We manually threshold this change probability map by using

T = 0.5 for industry areas and T = 0.9 for forest areas respectively. As forest areas display higher seasonal variability, a higher threshold value has to be chosen to detect the real forest changes. The comparison between the pixel-based and region-based change detection are displayed in Figs. 12 and 13. The first column shows the pixel-based change probability map, while the second column displays the region-based result. The produced change masks from the region-based method are shown in the third column. The datasets in forest area are quite complex because several change areas are quite small with irregular shapes. However, the result of our region-based change detection approach gives relatively good results in most of the test areas (shown in Fig. 12c1–c3). The datasets in industry areas are relatively straightforward due to large buildings, thus more precise change results are achieved. Tables 4 and 5 show the quantitative change detection results for all datasets. It is clear from the numerical results that the

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

233

Fig. 6. Cartosat-1 data in industrial areas. Figs. a1–a3 are the panchromatic images from year 2008; Figs. b1–b3 are the panchromatic images from year 2011; Figs. c1–c3 are the change reference data for each dataset.

region-based change detection approach leads to better results in both OA and KIA at 5 out of 6 test areas (about 10% improvement in KIA). Due to the high proportion of unchanged pixels in the overall scene the OA index is less sensitive, but it can be seen that especially for the industrial areas, the OA is significantly higher.

5. Summary and discussion In this paper, a novel approach to change detection using Cartosat-1 stereo images is presented. The approach is aiming to improve the change detection result by fusing all information that can be extracted from stereo images. Firstly, DSMs with 5 m resolution from two dates are generated with the semi-global matching method. Secondly, a region-based change detection workflow is proposed. Each region is characterized by a vector, which contains

the change features from both images and DSMs. The proposed method is quantitatively analyzed on six data sets, three forest areas and three industrial areas. The results show that, by using height information from stereo images, multi-level features and sharper boundaries from panchromatic images, our method is able to detect even changes in small areas in Cartosat-1 data correctly. The change detection results are more precise than using either DSMs or images alone. The proposed region extraction and merge procedure improves the detection speed and accuracy. The chosen change features and change detection methods suite well for the purpose. The outcome of this study is discussed in two steps, one is the generated DSM’s quality, and the other is the proposed change detection method. The DSM generation method explained in this paper is well suited for Cartasat-1. The quality of the DSM can be confirmed by the DSM quality assessment and by the change detection result

234

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

Fig. 7. DSM quality assessment test areas (field, forest and residential area).

assessment. As can be seen from Table 1, although some outliers still exist, the median and mean height errors in the generated DSMs are less than one pixel size. It indicates that the 5 m resolution Cartosat-1 DSMs can be used in small scale change analyses. Based on this accuracy, when only the change of high trees (height > 20 m) in forest areas is of interest, theoretically a change distribution can be obtained by DSMs subtraction alone. The DSM quality is also proven by the AUC comparison in Table 3. As shown there, height changes from DSMs have improved the change detection accuracy in all six test areas. It has to be noted, that the quality of the DSM can directly influence the change detection result. Although combined usage of images and DSMs can cope with some false height change alarms, when the area of outliers in DSM is larger than the defined regions

and the intensity value in the images of the area appear also different, all false alarms can be hardly prevented. As shown in ‘Test area 1’ of the forest areas, the deforestation is distributed in the middle of the forest, where the DSM quality is not optimal. In this case, our method can locate the changed places roughly, but the change boundaries cannot be defined exactly, and this also influences the region-based change detection accuracy (Table 4). The proposed two steps region-based change detection method is well suited for Cartosat-1 data. As it is an unsupervised change detection method, no initial training is required which makes this method very fast and flexible. Efficient and robust change detection can be performed in forest and industrial areas by combining height change with the spectral change from original pan-images. Although a pixel-based method is able to highlight most of the

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

235

Fig. 8. Profile analysis of three land covers (a) field, (b) forest and (c) residential area. (—: LiDAR-DSM, - - - -: Cartosat-DSM).

Table 1 Pixel and region based height difference between Cartosat-DSM and Lidar-DSM in metres. Method.

Max(abs)

Median

Mean

Std

Pixel based-Forest Pixel based-Ground Pixel based-Residential Region based-Forest Region based-Ground Region based-Residential

42 10 32.82 26.45 6.5 8.77

2 3 2.82 1.83 2.41 2.82

1.51 2.68 2.76 3.69 2.55 2.67

10.53 2.23 3.04 7.71 1.57 1.70

changes, it is more sensitive to noise in images and DSMs. Thus, in the pixel-based change detection results, some no-change areas also get high change probabilities, especially in forest areas. As illustrated in Fig. 12, this situation can be found at forest edges and some roads within the forest area. In addition, the advantage of a region-based method is also expressed in the extracted changed building boundaries. By using the region boundaries from panchromatic original images the changed building shape is often better preserved than in the DSM.

Fig. 9. Region merging procedure. (a) Original region map (926 segments); (b) after merging the small regions with less than 20 pixels (302 segments left); (c) after merging small regions with less than 50 pixels (171 segments left).

236

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

Fig. 10. Histogram of region size distribution for Fig. 9a.

Table 2 Comparison of the change map generated with and without the proposed region merge procedure (Higher accuracy is shown in bold font.). Forest area

Test area 1 Test area 2 Test area 3

Before region merge

After region merge

Region amount

KIA

Region amount

KIA

1359 1880 1151

0.4825 0.7698 0.6265

359 559 413

0.5273 0.7909 0.6936

Height change is helpful for more complicated/diverse areas. For example in forest area ‘Test area 2’, as can be seen in Fig. 6 a2 and b2, there is a lake in the lower-left corner. Due to season changes the colour of this lake is totally different in the two data sets. By including the height change information from DSMs, the AUC has improved from 0.9130 to 0.9874. In some cases, height information is indispensable for change detection, since it is often

very difficult or practically impossible to separate ground area and small trees using only intensity values especially within shadow areas (as shown in Fig. 14, marked with white circles). Also many crop areas share a similar intensity and texture information as some forest areas. In this case, without the multispectral information, height becomes the only information we can rely on. 6. Conclusion The main contribution of this research is to explore the applicability of Cartosat-1 stereo data for change detection even for objects like industrial buildings and minor forest regions. To achieve that, we have jointly used information from panchromatic images and height information coming from the DSMs derived from stereo imagery, and proposed a workflow which relies on a high quality DSM generation method, an efficient region merge procedure and a robust change detection method.

Fig. 11. ROC curve comparison in Test2-forest area.

237

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239 Table 3 Feature performance comparison (Higher accuracy is shown in bold font.). AUC (Industrial)

Image DSM Multi-level (image based) Multi-level + height (CVA) Multi-level + height (weighted CVA)

AUC (Forest)

Test 1

Test 2

Test 3

Test 1

Test 2

Test 3

0.9649 0.9567 0.9798 0.9878 0.9879

0.9560 0.8963 0.9728 0.9838 0.9860

0.9482 0.9185 0.9650 0.9721 0.9730

0.9399 0.8108 0.9801 0.9858 0.9845

0.9130 0.9361 0.9478 0.9830 0.9874

0.9809 0.8275 0.9962 0.9968 0.9954

Fig. 12. Change detection result in forest area. Figs. a1–a3: pixel-based change probability maps; Figs. b1–b3: region-based change probability maps; Figs. c1–c3: change detection masks (T = 0.90) (Blue: Change probability = 0; Red: Change probability = 1).

Without multi-spectral channels height information becomes especially important for Cartosat-1 change detection. Results based on the dense matching methodology, indicate that DSM generated

from Catosat-1 data can reach 5 m ground resolution, with a height accuracy of about 2–3 m. Since the DSM has its potential local errors, information from the original images like texture features

238

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

Fig. 13. Change detection result in industry area. Figs. a1–a3: the pixel-based change probability maps; Figs. b1–b3: region-based change probability maps; Figs. c1–c3: change detection masks (T = 0.50) (Blue: Change probability = 0; Red: Change probability = 1).

Table 4 Performance result for forest areas (Higher accuracy is shown in bold font.). Forest area

Test area 1 Test area 2 Test area 3

Pixel-based

Region-based

Table 5 Performance result for industrial areas (Higher accuracy is shown in bold font.). Industry area

OA

KIA

OA

KIA

0.9963 0.9950 0.9960

0.5721 0.6032 0.4580

0.9952 0.9973 0.9975

0.5273 0.7909 0.6936

and object boundaries, are used to overcome some false alarms. The proposed region-based method is a good choice in combining all of these features. Although the DSMs from Cartosat-1 are not precise, they can be a reliable source for achieving high

Test area 1 Test area 2 Test area 3

Pixel-based

Region-based

OA

KIA

OA

KIA

0.9803 0.9700 0.9735

0.7994 0.7640 0.6985

0.9854 0.9887 0.9845

0.8597 0.9131 0.8131

change detection accuracies for forest and industrial areas. Changes of these areas are very important for many applications, for instance, forest management, environmental evaluation and city planning.

J. Tian et al. / ISPRS Journal of Photogrammetry and Remote Sensing 79 (2013) 226–239

239

Fig. 14. Shadow in forest area. (a) Pancromatic image of date1; (b) pancromatic images of date2.

Furthermore deriving 3D change information from satellite imagery like Cartosat-1 is much more economical than from airborne LiDAR, InSAR or airborne stereo imagery, since the latter are by far more costly and need more logistics in data acquisition. The study confirms the high capability of the generated DSMs from Cartosat-1 data in assisting change detection. The method could be easily adapted for other satellite and airborne stereo imagery. The DSM’s quality can be improved when higher resolution stereo images are available, especially for building areas. If available, multi-spectral channels can also be used to enrich the change features for the weighted CVA. Thus, more precise results should be achieved, and changes in dense residential areas could potentially also be detected following this approach. This will be a matter of a further research study.

References Aguirre-Gutiérrez, J., Seijmonsbergen, A.C., Duivenvoorden, J.F., 2012. Optimizing land cover classification accuracy for change detection: a combined pixel-based and object-based approach in a mountainous area in Mexico. Applied Geography 34, 29–37. Berberoglu, S., Akin, A., 2009. Assessing different remote sensing techniques to detect land use/cover changes in the eastern Mediterranean. International Journal of Applied Earth Observation and Geoinformation 11 (1), 46–53. Blaschke, T., 2010. Object based image analysis for remote sensing. ISPRS Journal of Photogrammetry and Remote Sensing 65 (1), 2–16. Bovolo, F., 2009. A multilevel parcel-based approach to change detection in very high resolution multitemporal images. IEEE Geoscience and Remote Sensing Letters 6 (1), 33–37. Bruzzone, L., Prieto, D.F., 2000. An adaptive parcel-based technique for unsupervised change detection. International Journal of Remote Sensing 21 (4), 817–822. Celik, T., 2009. Unsupervised change detection in satellite images using principal component analysis and k-means clustering. IEEE Geoscience and Remote Sensing Letters 6 (4), 772–776. Chen, G., Hay, G.J., Carvalho, L.M.T., Wulder, M.A., 2012. Object-based change detection. International Journal of Remote Sensing 33 (14), 4434–4457. Comaniciu, D., Meer, P., 2002. Mean-shift: a robust approach toward feature space analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (5), 603–619. Coppin, P.R., Bauer, M.E., 1996. Digital change detection in forest ecosystems with remote sensing imagery. Remote Sensing Reviews 13 (3–4), 207–234. Cormen, T.H., Leiserson, C.E., Rivest, R.L., Stein, C., 2001. Introduction to Algorithms (2nd ed.). Greedy Algorithms, MIT Press and McGraw-Hill (Chapter 16). d’Angelo, P., 2010. Image matching and outlier removal for large scale DSM generation. Convergence in Geomatics. CGC & ISPRS. ISPRS Symposium Commission I, Calgary, Canada, 15–18 June. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 38 (Part 1), (on CDROM). d’Angelo, P., Lehner, M., Krauss, T., 2008. Towards automated DEM generation from high resolution stereo satellite images. International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences (Part B4) 37, 1137–1142. Desclée, B., Bogaert, P., Defourny, P., 2006. Forest change detection by statistical object-based method. Remote Sensing of Environment 102 (1–2), 1–11. Duveiller, G., Defourny, P., Desclée, B., Mayaux, P., 2008. Deforestation in Central Africa: estimates at regional, national and landscape levels by advanced

processing of systematically-distributed Landsat extracts. Remote Sensing of Environment 112 (5), 1969–1981. Grodecki, J., Dial, G., Lutes, J., 2004. Mathematical model for 3D feature extraction from multiple satellite images described by RPCs. In: Proc. ASPRS Annual Conference, Denver, Colorado, USA, pp. 23–28 May, (on CDROM). Grohman, G., Kroenung, G., Strebeck, J., 2006. Filling SRTM voids: the delta surface fill method. Photogrammetric Engineering and Remote Sensing 72 (3), 213–216. Hand, D.J., Till, R.J., 2001. A simple generalization of the area under the ROC curve to multiple class classification problems. Machine Learning 45 (2), 171–186. Hirschmüller, H., 2008. Stereo processing by semiglobal matching and mutual information. IEEE Transactions on Pattern Analyses and Machine Intelligence 30 (2), 1–14. Im, J., Jensen, J.R., Tullis, J.A., 2008. Object based change detection using correlation image analysis and image segmentation. International Journal of Remote Sensing 29 (2), 399–423. Johnson, R.D., Kasischke, E.S., 1998. Change vector analysis: a technique for the multispectral monitoring of land cover and condition. International Journal of Remote Sensing 19 (3), 411–426. Kamini, J., Jayanthi, S.C., Raghavswamy, V., 2006. Spatio-temporal analysis of land use in urban Mumbai – using multi-sensor satellite data and GIS techniques. Journal of Indian Society of Remote Sensing 34 (4), 385–396. Klonus, S., Ehlers, M., Tomowski, D., Michel, U., Reinartz, P., 2011. Detektion von zerstörten Gebäuden in Krisengebieten aus pan-chromatischen Daten. Photogrammetrie Fernerkundung Geoinformation 2011 (4), 219–231. Kumar, K.V., Martha, T.R., Roy, P.S., 2006. Mapping damage in the Jammu and Kashmir caused by 8 October 2005 Mw 7.3 earthquake from the Cartosat–1 and Resourcesat-1 imagery. International Journal of Remote Sensing 27 (20), 4449– 4459. Lu, D., Mausel, P., Brondizio, E., Moran, E., 2004. Change detection techniques. International Journal of Remote Sensing 25 (12), 2365–2407. Malila, W.A., 1980. Change vector analysis: an approach for detecting forest changes with Landsat. In: Proceedings of the 6th Annual Symposium on Machine Processing of Remotely Sensed Data, Purdue University, West Lafayette, Indiana, 03–06 June, pp. 326–335. Martha, T.R., Kerle, N., Jetten, V., van Westen, C.J., Kumar, K.V., 2010. Landslide volumetric analysis using Cartosat-1 derived DEMs. IEEE Geosciences and Remote Sensing Letters 7 (3), 582–586. Mas, J.-F., 1999. Monitoring land-cover changes: a comparison of change detection techniques. International Journal of Remote Sensing 20 (1), 139–152. Prabaharan, S., Srinivasa Raju, K., Lakshumanan, C., Ramalingam, M., 2010. Remote sensing and GIS application on change detection study in coastal zone using multi temporal satellite data. International Journal of Geomatics and Geosciences 1 (2), 159–166. Radke, R.J., Andra, S., Al-Kofahi, O., Roysam, B., 2005. Image change detection algorithms: a systematic survey. IEEE Transactions on Image Processing 14 (3), 294–307. Rao, C.V., Sathyanarayana, P., Jain, D.S., Manjunath, A.S., 2007. Topographic map updation using Cartosat-1 data. In: Proceedings of RSPSoc Annual Conference 2007, Newcastle upon Tyne, UK, 11–14 September, (on CDROM). Rodriguez, E., Morris, C.S., Belz, J.E., Chapin, E.C., Martin, J.M., Daffer, W., Hensley, S., 2005. An assessment of the SRTM topographic products. Technical Report JPL D31639, Jet Propulsion Laboratory, Pasadena, California. Tian, J., Leitloff, J., Krauss, T., Reinartz, P., 2011. Region based forest change detection from Cartosat-1 stereo imagery. ISPRS Hannover Workshop 2011: HighResolution Earth Imaging for Geospatial Information, Germany, Hannover, 14–17 June, International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 38 (Part 4/W19), 319–324. Tian, J., Cui, S., Reinartz, P., 2013. Building change detection based on satellite stereo imagery and Digital Surface Models. IEEE Transactions on Geosciences and Remote Sensing, in press (doi: (identifier) 10.1109/TGRS.2013.2240692). Walter, V., 2004. Object-based classification of remote sensing data for change detection. ISPRS Journal of Photogrammetry and Remote Sensing 58 (3–4), 225–238.