OBJECT-BASED FOREST CHANGE DETECTION USING HIGH RESOLUTION SATELLITE IMAGES

An object-based approach for forest disaster change detection using High Resolution (HR) satellite images is proposed. An automatic feature selection process is used to optimize image segmentation via an original calibration-like procedure. A multitemporal classification then enables the separation of wind-fall from intact areas based on a new descriptor that depends on the level of fragmentation of the detected regions. The mean shift algorithm was used in both the segmentation and the classification processes. The method was tested on a high resolution Formosat-2 multispectral satellite image pair acquired before and after the Klaus storm. The obtained results are encouraging and the contribution of high resolution images for forest disaster mapping is discussed.


INTRODUCTION
In a climate changing context, wind storms have become more and more frequent.Wind-fall damages have to be quickly mapped to prevent fire risks, call for financial compensation and to update the national forest inventory.While ground investigations are complex due to fallen trees, remote-sensing techniques enable fast monitoring of large and unreachable areas.Their use have widely spread through the world for disaster change detection, especially with the growing spatial and temporal resolutions of new satellites.The objective of the present study is to provide a binary mapping discriminating damaged and intact areas using HR satellite imagery consisting of a pair of multispectral Formosat-2 images of 8 m resolution.The images were taken before and after the Klaus Storm that happened the 24 th January 2009 in the South West of France.In literature, the former works in forestry produced low scale maps, near to the hectare and studied essentially the clear-cuts.Few works have dealt with smaller structural changes such as wind-fall damages.Recently, object-based classifications were used for change detection in forestry (Desclee et al., 2006, Conchedda et al., 2008).These methods are based on a segmentation process that combines spatial and spectral information to group pixels into homogeneous regions before their classification using new object descriptors.The latter can be geometrical and textural (Fraser et al., 2005) or temporal (Desclee et al., 2006).Besides, change detection can be either based on a comparison of before and after storm classifications (post-classification), or directly processed on multitemporal images (joint-classification).Post-classification approaches are robust to radiometric differences between images, and provide an accurate "from-to" change information (Im and Jensen, 2005) but suffer from errors propagation.Joint-classification approaches provide more information to classify small changes.They can either rely on machine learning algorithms using a training set (Im and Jensen, 2005) or use thresholding which involves a parametric statistical test (Desclee et al., 2006) or an expert knowledge (Fraser et al., 2005).
In this study, the proposed method is an object-based, multitemporal classification that maps storm damages at a fine spatial scale.We propose a nearly automatic method requiring very limited data for rapid mapping at a regional scale.An unsupervised multitemporal classification was preferred to a post-classification scheme since, in our context, changes are subtle and hard to detect in a single after-storm image.The algorithm is based on the mean shift segmentation that will be detailed in section 2.An automatic feature selection process is used to optimize image segmentation via an original calibration-like procedure and will be presented in section 3. Section 4 presents the binary multitemporal classification which is based on the mean shift algorithm and uses a new descriptor that depends on the level of fragmentation of the detected regions.Experimental results are shown and discussed in section 5 and finally conclusions are drawn.

MEAN SHIFT SEGMENTATION
The Mean Shift (MS) algorithm is a non-parametric feature-space analysis technique.(Fukunaga andHostetler, 1975, Comaniciu andMeer, 2002) showed excellent results in clustering and object delineation in color images.It is based on a density mode searching and clustering technique.The feature space is considered as the empirical probability density function (p.d.f.) of the input features.The algorithm proposes a filtering step that associates each pixel in the image with the closest local mode in the density distribution of the feature space.The MS procedure actually locates theses modes without estimating the global density.The segmentation into a piecewise constant structure requires one more step, the fusion of the regions associated with nearby modes.The implementation of (Comaniciu and Meer, 2002) searches for local modes in the joint feature and spatial domain of n + 2 dimensions, where n is the number of considered features.An iterative procedure of mode seeking consists in shifting the n+2 dimensional window to a local mode.The search window involves two user-defined inputs that can be deduced from desired object sizes or physical properties.A radiometric range (hr) corresponds to the unique spectral radius in the n-dimensions search window and a spatial bandwidth (hs) corresponds to the spatial radius of the window.
Figure 1 shows the impact of both parameters on segmentation results.
Figure 1: Mean shift segmentations of a 4-channel image (B,G,R,NIR) using different parameters (hs, hr).When hr increases, only highly contrasted and homogeneous regions remain: intact young stand on the left, blue farming area, medium aged stand on the right.If hs increases, only larger regions remain : shadow areas and standing tree groups disappear within the central, strongly damaged area.
Unlike hierarchical methods used in (Desclee et al., 2006, Conchedda et al., 2008), the Mean Shift does not use any heterogeneity measure which is more appropriate in forestry where the objects of interest, i.e. tree stands, include heterogeneous pixels such as vegetation, ground and shadows.To our knowledge, the MS algorithm has not been used for forestry mapping.In our method, the Mean Shift procedure is used in both segmentation and classification steps using the joint spatio-spectral domain and solely the spectral domain respectively (cf.Section 4).

Feature selection using test frames
In literature, the segmentation is usually processed on all available before and after image bands.In this study, in order to determine the most relevant features for segmentation, input feature selection was carried out through an original generic calibrationlike procedure using a test frame.Moreover, this process aims to automatically optimize the segmentation parameters.The test frames are constructed with n small image samples corresponding to n classes, yielding as much test frames as input features, thus resulting in a multi-band test frame, at the core of our feature selection method.The used data is not included in the validation dataset.The test frames are then segmented by the MS algorithm using either one input feature (single test frame) or multiple normalized features (multi-band test frame) while testing various parameterizations.The feature Segmentation Performance (SP) is defined as: where A is the area, nc the actual number of classes in the test frame, Ri the actual regions and RS j the segmented regions.The segmentation performance depends on the regions overlap percentage.Figure 2 illustrates an example of a test frame with three regions and the SP computation.
The highest SP value leads to the best feature or set of features.

Proposed input features
Forest features can be grouped into three categories: spectral, textural and temporal features.Spectral features include raw or corrected spectral bands and derived bands such as the NDVI (Normalized Difference Enhanced Vegetation Index) (Conchedda et al., 2008), the PCI (Principal Components Inversion) (Inglada, 2001), etc.In this study, the four Formosat-2 spectral bands (blue, green, red and infrared) were used as well as two vegetation indexes, the Normalized Difference and Soil Adjusted Vegetation Indexes (NDVI and SAVI).First order statistics such as mean or variance of reflectance are also used.As for textural features, the more common ones are Haralick features (Haralick et al., 1973) as used in (Ruiz et al., 2004, Kayitakire et al., 2006, St-Louis et al., 2006, Trias-Sanz et al., 2008, Tuominen and Pekkarinen, 2005).Finally, three common temporal features were computed: mean correlation, difference and ratio between both images.More temporal features can be found in literature such as the pixel-wise Magnitude of the Temporal Change Vector (Fraser et al., 2005), the multi-temporal PCI (Inglada, 2001), the Neighborhood Correlation Index (Im and Jensen, 2005) and will be investigated in a future work.Temporal and textural features were then processed for each spectral feature.

UNSUPERVISED OBJECT-BASED MULTITEMPORAL CLASSIFICATION
The global unsupervised object-based multitemporal change detection scheme is depicted in Figure 3.The before and afterstorm images are segmented independently using the MS algorithm and the feature selection process as explained in Section 3. The before-storm segments correspond to homogeneous structural regions, i.e belonging to the same age class.The after-storm segments reflect the change degree.

Mean Shift spectral classification
The

Fragmentation rate
Damaged areas are heterogeneous and therefore appear over-segmented in the after-storm image.Conversely, intact areas correspond to larger regions and have similar delimitations in both images.The Fragmentation Rate (FR), characterizes the beforestorm regions and reflects their over-segmentation in the afterstorm image.It is computed by a comparison between before and corresponding after-storm regions as follows: where A is the area, R is the before-storm region and R a j the after-storm regions that are included (partly or entirely) in the before-storm one.
The average FR is then computed for the change clusters afterstorm.The change class fragmentation rate is then defined as : where RC is the change cluster, p and N are respectively the pixels and the number of pixels of the change class C, Ri are the before-storm regions.The clusters are finally divided into two intact and damaged classes, based on their fragmentation rate and using the unsupervised Otsu thresholding (Otsu, 1979) which minimizes intra-class variances.

Study area and data set
The Nezer forest is located on the French Atlantic coast.It is made up of rectangular stands of pine trees that have the same age and height.Available images of the area are a set of orthorectified, geo-referenced multispectral Formosat-2 images before and after the Klaus storm, acquired on 22/12/08 and 04/02/09, respectively (cf. Figure 4).The images have a 8m spatial resolution and four spectral bands (Blue, Green, Red and Near Infrared).Ancillary data include tree stand delimitations and ages GIS layer and 100 reference areas identified on orthophotos of 15 cm resolution dating back to 26/02/09.

Feature Selection for segmentation and classification
Table 1 outlines the used input features, separated into three groups: spectral, textural and temporal.A total of 84 features are used: 6 spectral features, 10 textural and 3 temporal features all computed on the six spectral features respectively.For the textural features, the neighborhood radius and the directional vector offsets were both experimentally set to 1 pixel, the displacement vector being horizontal.
These features were used in the feature selection process, individually or combined in a multi-feature test frame after their normalization.A test frame composed of four regions was used.The spatial radius hs of the MS segmentation was set using a prior thematic knowledge on the desired objects size.A 3-pixel area can reasonably be considered as a regular forest pattern, so hs was set to 3. To set the radiometric range hr, all features were individually rescaled between 0 and 255 and spectral ranges from 2 to 60 were experimented via our automated Mean-Shift parameter optimization procedure.
In the MS spectral classification, the smallest spectral range hr=2 was used in order to detect subtle changes.Besides, unlike the MS segmentation, this spectral-based clustering stage involves only object temporal features.They derived from averaging the temporal features over the segmented after-storm regions.Tables 2 and 3 present the best features and spectral ranges for image segmentation and classification respectively.
After In this study, the optimal features are essentially spectral or temporal.Textural features should give better results on panchromatic images with a spatial resolution of 2m.Indeed, the relatively low spatial resolution (8m) of the Formosat-2 multispectral images turns out insufficient to properly exploit the textural information, particularly of significance in the context of forestry.
Besides, the Haralick features were processed on small windows (3 × 3) which may be insufficient to capture the stand texture.Moreover, in some cases, the broken stems can be oriented in similar directions which might lead to maxima on texture features.Therefore, texture features should be more useful on higher resolution and with an optimization of Haralick parameters.
For the after-storm segmentation process, the best feature is temporal.Indeed, the stand structure (which depends on the tree age) before the storm helps to determine the change degree after the storm.In addition, one can observe that the temporal features lead to higher spectral ranges than the original image bands.
Temporal features present a Gaussian distribution, a large spectral range allows to segment the image into different change degrees, whereas the initial image bands present a higher variability, a finer spectral range is necessary to segment subtle changes.
Multiple feature segmentations were also tested.In our experiments, the segmentation performances were better using one feature only.This rather unexpected result can be explained by two reasons.First, the used implementation of the Mean Shift (Comaniciu and Meer, 2002), involves one unique spectral range hr for multiple features.Indeed, it was initially proposed for grayscale and color image segmentation.An adaptive spectral range per feature, should enhance the results as it is more appropriate for remote sensing images where spectral distributions of image bands are different.Secondly, the forest canopies are very complex and their variability depends on various parameters or features that are not correlated to the damage degree.For instance, the classical use of the eight initial bands (i.e before and afterstorm images) for the segmentation process led to a decrease in the global classification accuracy of 5% (Orny et al., 2010).

Segmentation and fragmentation rate
Figure 5 depicts the segmentation results before and after the storm using the respective best features.One can visually distinguish intact and damaged areas.Intact areas are larger and have the same delimitations in both segmentations, whereas the damaged areas are more heterogeneous leading to an over-segmentation into many small regions.Figure 6 shows the fragmentation rate (FR) of before-storm segments in gray levels.The lighter the regions, the more damaged  they are.This result matches visually the tree stand age map where the older stands appear to have more damage extent than the younger stands.In fact, among numerous factors, the tree height influences the sensitivity of the tree to the wind.The young stands are dense with small trees which make them more robust to the wind.On the contrary, older stands are less dense, more heterogeneous due to sylvicultural practices, and present higher trees that are more vulnerable and likely to be damaged by storms.
Figure 7 shows the obtained change clusters using the MS spectral classifier.About 30 change clusters are obtained.Given the histogram complexity of this image, it is difficult to separate clusters into two classes automatically.On the contrary, one can observe that after characterizing these regions by their average fragmentation rate, the damaged areas are better discriminated with respect to the intact areas.

Map validation
The figure 8 illustrates the final binary map for forest disaster detection with intact and damaged areas.Reference data were collected from orthophotos of 15 cm resolution, that were available after the storm.
The Overall accuracy was of 87.2%.Some small changes or intact areas were not extracted due to the limited spatial resolution of Formosat-2.Besides, in our method, the shadows are not taken into account, therefore the segmentation was not robust to shadow changes.
The INRA 1 inventory data layer of the Nezer site references ages of all forest stands.It was confronted to the obtained classification.The overall accuracy and the percentage of damaged pixels (damaged rate) were computed for all available ages.The obtained results are depicted in figure 9.
One can observe that the classification accuracy increases significantly with the class age.The tree height influences its sensitivity to the wind (Cucchi et al., 2005).The damage rate is higher for the older stands and reaches 70% for stands whose age is superior to 25 years.
The detection precision of stands aged from 14 to 39 years (four intermediate classes) is high and ranges between 93.3 and 99.4%.
1 Institut National de Recherches Agronomiques However, the two youngest and the oldest stands have a lower detection rate combined to high omission and commission errors (cf.Table 4).The confusion happens between damaged areas and old stands that are heterogeneous and sparse.In addition, small damaged areas or areas with a low intensity damage (leaning trees) are difficult to detect in young dense stands.
These results can be compared to the results of a similar work obtained in (Schwarz et al., 2001) from 10m high resolution multispectral images using a supervised object-oriented approach.Using an automatic object-based approach, our overall accuracy is inferior (87.2%) to the one obtained by (Schwarz et al., 2001) (96%).However, the classification accuracy obtained on trees, aged from 14 to 39 years, is slightly better (96% versus 95%).Good results are obtained on intermediate and older classes.Some problems persist on younger classes as explained above.However, our method requires only a few samples to construct the test frame and only two parameters to tune.Reducing the training and tuning time is essential for emergency mapping.

CONCLUSION
We presented, in this paper, an object-based multitemporal change detection method, well-suited for emergency mapping.Our contribution provides two main novelties with respect to similar works in forestry.Firstly, an automatic feature selection process, applied to both the segmentation and classification steps, is introduced, whose originality lies on the use of test frames (single or multi-bands) of adequate forest samples.This innovative feature selection process, inspired by camera calibration procedures, allows a rapid evaluation of hundreds of features and combined features.It is applicable to the optimization of any other segmentation algorithm and in the context of any application related to forestry or not.The second originality of our approach is a relevant fragmentation rate, dedicated to storm damages, that ensures an automatic thresholding of the change clusters into a binary map.Our method gives good results.It has the appealing property of requiring only a few samples to construct the test frames, leading to a slight supervision, compared to more traditional supervised methods, hence categorizing it as an unsupervised approach.Moreover, our method involves only two user-defined parameters and does not need any statistical assumption, thanks to the powerful Mean-Shift clustering algorithm, at the core of our change detection scheme.
However, the Formosat-2 multispectral images resolution appears to be not well-suited to detect scattered small damages, due to the underlying relatively low spatial resolution (8m).Indeed, this resolution turns out insufficient to properly exploit the textural information, particularly significant in the context of forestry.Panchromatic images, of 2m resolution, have the potential to enhance the current mapping performances.Finally, the confrontation of the resulting binary (intact/damaged) classification to the age class map was consistent and confirmed the increase in sensitivity of the tree to the wind with the age.
Future work will be devoted to the application of our change detection scheme to panchromatic and very high resolution images of the Nezer site.A higher spatial resolution would be certainly more appropriate to capture textural forest details and improve the detection of small changes.A hierarchical approach could also enhance the detection of subtle changes in young dense stands.In addition, in the Mean Shift algorithm, the spectral range hr should be adapted to each feature, to take advantage of the spectral distribution variability and the combination of multiple features.
after-storm segmented regions are represented by object mean temporal descriptors.The automatic feature selection process (cf.Section 3.1) provides the input features that optimize the MS spectral classifier.The segmented regions are then clustered automatically into change classes using this optimized mean shift spectral classifier.Unlike the MS segmentation, this modified version is independent of pixel positions and involves solely the spectral domain which allows to cluster similarly damaged regions that are spatially distant into the same change class.The MS classifier has a single parameter, hr, which is easy to specify.The MS spectral classification leads automatically to many change clusters.No reference data were available to validate the obtained change degrees.Consequently, our objective was limited to the production of a binary change map even if the MS spectral classification provides multiple change classes.In order to group the change classes into intact and damaged classes, the clusters were characterized by an innovative temporal feature : the fragmentation rate.
Figure 4: Formosat-2 multispectral images acquired before and after the Klaus storm

Figure 5 :
Figure 5: Segmentation of before and after-storm images using the best features, i.e red band and red band ratio respectively.

Figure 6 :
Figure 6: Comparison between the fragmentation rate and the tree stand ages.
Figure 7: MS spectral classification and fragmentation rate

Figure 8 :
Figure 8: Binary map of wind storm damaged areas.

Figure 9 :
Figure 9: Classification accuracy and damaged pixel rate with respect to age classes.

Table 1 :
Input features

Table 2 :
-storm segmentation Before-storm segmentation Feature hr SP (%) Feature hr SP (%) Optimal features and spectral ranges hr for the afterstorm and the before-storm segmentations using a test frame of 4 samples.SP is the segmentation performance (cf.Eq. 1).

Table 3 :
Optimal features for the binary classification using a test frame of 4 samples with hr = 2.
table 4 shows the global confusion matrix, obtained by comparing pixel values between the classification results and the reference data.