Change Detection in Boreal Forests Using Bi-Temporal Aerial Photographs

Increased need for timely forest information is leading to continuous updating of stand databases. In continuous updating, stand attributes are estimated in the field after an operation and stored in databases. To find the changes caused by operations and forest damage, a semi-automatic method based on bi-temporal aerial photographs was developed. The test data were classified into three classes: No-change (952 stands), Moderate-change (163 stands) and Considerable-change (44 stands). The aerial photographs were acquired in years 2001 and 2004 with almost the same image specifications. Altogether 110 features at stand level were extracted and used in change detection analysis. The test data were classified with stepwise discriminant analysis. The overall accuracy of classification varied between 75.3 and 84.7%. The considerable changes were found without error, whereas the Moderate-change and No-change classes were often confused. However, 84.2% of thinned stands were classified correctly. The best accuracy in classification was obtained by using the histogram and textural features extracted from the original, uncorrected images. Radiometric correction did not improve the accuracy of classification. Soil type, characteristics of the growing stock and the location of a stand in an image were found to affect the change detection. Before the method can be applied operationally, issues related to, e.g., confusion between No-change and Moderate-change must be solved.


Introduction
In Finland, stand attributes have traditionally been gathered by periodical field inventories in which the inventory cycle has been 10-15 years.Increased need for timely forest information is, however, leading to continuous updating of stand databases.In continuous updating, a stand database is kept up-to-date computationally using statistical growth models.After an operation (i.e.cutting or silvicultural treatment), stand attributes are estimated in the field.
An operation can be registered during the work, but, e.g., forest damage must be found some other way.There are also errors that should be controlled in the databases.Medium-resolution satellite images (e.g.Landsat TM images) have been successfully applied in detecting considerable changes like clear cuttings, removals of hold-over trees, soil preparations or drastic damage, while the detection of moderate changes like thinnings, preparatory cuttings or slight damage, has been difficult (e.g.Häme 1991, Varjo 1996, Holmgren and Thuresson 1998, Wilson and Sader 2002, Saksa et al. 2003, Heikkonen and Varjo 2004).The reason for this is that typical thinnings, where about 20-40% of basal area is removed, cause only subtle change in reflectance (Häme 1991, Olsson 1994).
The reflectance of a single pixel with mediumresolution images is an average of the reflectance from an area of more than about 100 m 2 on the ground.With high-resolution remote sensing material, even disappearance of individual trees can be detected.For example, using airborne laser scanning (ALS) data, Yu et al. (2004) found 61 of 83 harvested trees.Compared to aerial photography, however, ALS is very expensive.
Operational high-resolution applications for change detection of vegetation cover utilize visual interpretation of aerial photographs (see e.g.Fensham and Fairfax 2002), although more automatic methods have also been proposed.Hudak and Wessman (1998) investigated a transition from grassland to shrubland using historical aerial photographs.They classified the images using variograms that characterized the image texture.Changes were determined by post-classification comparison.Kadmon and Harari-Kremer (1999) also studied vegetation dynamics.Their approach comprised pixel-level spectral classification, averaging the results to larger cells and differencing the cell values.In both studies, good results using automatic change detection were achieved; but compared to the context of this study, the timeintervals were much wider and the changes more radical.In the study of Saksa et al. (2003), clear cuttings were detected by using three approaches: pixel-by-pixel differencing and segmentation, pixel block-level differencing and thresholding, and presegmentation and unsupervised classification.All the methods were found to be suitable for operational use.Hyppänen (1999) applied image differencing to bi-temporal aerial photographs.In his study, considerable changes were detected while moderate changes were not.Consequently, the problem of how to detect moderate changes automatically using aerial photography has remained unsolved.
In automatic change detection, the effect of disturbing factors should be minimized.These factors include differences in atmospheric conditions, sun angle, soil moisture, vegetation phenology and, in the case of aerial photographs, the differences in viewing angles (Singh 1989, Hyppänen 1999).Especially factors that might cause differences between stands in different parts of the images should be eliminated.By using the present satellite positioning systems, aerial photographs can be shot very close to each other in different years.
This study had two objectives.The first was to investigate whether bi-temporal aerial photographs taken with similar image specification are beneficial in change classification, especially in finding moderate changes like thinnings.The second was to test the effect of the following factors on the accuracy of classification: image resolution, radiometric correction, extracted features and stand properties.

Field Data
The study area is located in Western Finland near the town of Kauhajoki (22°18´E, 62°24´N).The area is mainly owned by private ownerships.The forest holdings are formed in so called "stripsharing", and therefore the shape of these holdings is usually like long and narrow rectangle.The landscape of the area is dominated by flat terrain.The elevation varies between 125 and 185 m above sea level.The main tree species are pine (Pinus sylvestris L.) and spruce (Picea abies (L.) Karst.).Main broadleaved species are silver birch (Betula pendula Roth.) and pubescent birch (Betula pubescens Ehrh.)There are also sparse stands in places.Main site types are fresh (26% of the area), dryish (42%) and dry sites (25%).Site types are according to Kuusela and Salminen (1969) and the fertility of peatlands is determined and classified using the same system as that for mineral soils (Laine and Vasander 1993).
The field data consisted of 2361 forest stands (Table 1), for which the stand attributes were measured during summer 2002 based on a visual stand-level inventory system used in Finland (Solmu metsäsuunnittelun… 2000).In this system the stands are first pre-delineated by visual interpretation of aerial photographs.Secondly, the forest characteristics of the stands are estimated visually with the aid of some measurements in the field, and the delineation is confirmed.
The changes in stands between 2001 and 2004 were collected from different databases of the Regional Forest Center and all stands were checked visually using bi-temporal aerial photographs.The uncertain cases were checked in the field.Subsequently, all stands were divided into three classes according to the type of change.
No-change class included stands with no changes other than growth (1913).Moderate-change class included stands with moderate changes (350), i.e. thinning, seed tree felling, tending of seedling stand, improvement of young stands, removal of hold-over trees, soil preparation, slight storm damage, partly operated stand and forest road building.Slight storm damage means that the storm has felled only few trees on the stand.Considerable-change class included stands with considerable changes, i.e. clear cutting and severe storm damage (98).Severe storm damage means that the storm has felled several trees or group of trees on the stand.

Image Data
The image data consisted of two aerial photographs covering the study area.The idea was to take the latter photograph as close to the former one as possible with respect to time, date and location (Table 2).The photographs were acquired using Leica RC30 camera and an antivignetting and an infrared radiation filter.The nominal scale of the color-infrared photographs was 1:30 000.The photographs were first scanned with a photogrammetric scanner at 14 μm resolution into RGB-images and no tone adjustment was done.Secondly, the images were ortho-rectified to a spatial resolution of 0.5 m by using 11 (1422) and 13 (755) control points that were located from digital base maps and digital elevation model.The resolution of the digital elevation model was 25 m.The root mean square errors (RMSEs) of the rectification were 5.1 m (1422) and 3.2 m (755).These photographs were used as original photographs in feature extraction.Thirdly, because the other radiometric correction technique was computationally heavy (discussed closer in the following sub-section), the photographs were resampled to the spatial resolution of 3 m.This has been shown to be suitable resolution for numerical interpretation of aerial photographs in boreal forests (Hyppänen 1996).

Radiometric Correction
Aerial photographs are acquired from low altitudes and using wide-angle lenses.This causes strong bi-directional reflectance effect which can be seen within an image in the form of variation in brightness that is not due to variation in stand properties.In order to decrease the effect of varying bi-directional reflectance, the aerial photographs were corrected radiometrically by two different methods: by rationing and by an empirical correction chain.With aerial images, good results have been achieved by rationing (e.g.King 1991, Holopainen andWang 1998).When channel ratios are used, it is assumed that the sensitivity of the sensor, or the film and the scanner, is similar on different channels and that the ratio is similar on the different parts of the image.Three new images per date were obtained by calculating pixel-by-pixel the following ratios: IR/R, IR/G and R/G, where IR is the near-infrared channel, R the red channel and G the green channel.
In the empirical correction chain, the radio-metric correction for every pixel in a photograph is calculated by the functions depending on the pixel's location on the photograph, solar elevation and azimuth angles, flight course and altitude and scattering angles (Pellikka 1998, Pellikka et al. 2000).The first step was to reduce the light falloff effect.First, the zenith viewing angle was calculated.Secondly, because the calibration method is sensitive to vegetation type, a reference vegetation type was selected (Pellikka 1998).The dominant tree species in the study area was pine, and most of the interest was focused on thinned stands.Therefore, the reference areas were selected from middle-aged, pine-dominated stands.The reference areas were situated in the central part and in the corners of the aerial images.The spectral means of the reference areas were calculated, and the correction factor n, n-factor, for each channel was determined: where DN cio is the corrected spectral mean on channel c in reference area i, DN ciq is the spectral mean on channel c in reference area i, cn is the correction factor on channel c and q i is the zenith viewing angle in reference area i.
The value of n-factor for each channel was determined in such a way that the correlation coefficient between corrected spectral means and the zenith viewing angles of the reference areas were zero.The n-factor is related to the reduction in film exposure for an off-axis point (Lillesand et al. 2004p. 68, Pellikka 1998).The n-factor and zenith viewing angle were then used to correct the light falloff effect of the aerial photographs.
The next step was to reduce the bi-directional reflectance effect.First, the scattering angle was calculated.Secondly, six reference areas overlapping three different aerial photographs were selected, and in order to determine the correction factor f, f-factor, for each channel, the spectral means of the reference areas in every aerial photograph and channel were calculated: where DN* ci is the corrected spectral mean on channel c in reference area i, DN ci is the spectral mean on channel c in reference area i, f c is the correction factor on channel c and g i is the scattering angle in reference area i.The value of f-factor for each channel was determined in such a way that the correlation coefficient between corrected spectral means and the scattering angles of the reference areas were zero.
The f-factor is related to the Mie and Rayleigh scattering and is used to reduce the effect of the correction (Pellikka 1996).Thirdly, the scattering angle and the f-factor were used to decrease the effect of bi-directional reflectance.

Feature Extraction
Due to the manual stand delineation, the stand borders were often displaced from their "true" locations by several meters.Therefore, the borders were buffered by 5 m to decrease the effect of location errors.Spectral and textural features were extracted for these buffered stands resulting in 11 different data sets (Table 3).The median of the standwise DNs was chosen to represent the spectral characteristics of each channel, which further reduces the effect of displaced stand borders (Anttila 2002).
The following features were calculated to describe the histograms of DNs: median, standard deviation, skewness and kurtosis.In addition, an error index (modified from Reynolds et al. 1988), measuring the difference of histograms at the two points of time, was derived: where f i01 is the frequency of pixels with DN value i within the stand in 2001, f i04 is the corresponding frequency in 2004, and # pix is the number of pixels within the stand.
The histogram features are first order statistics that also describe, to some extent, the texture of a stand.To further characterize the texture, 11 texture measures based on gray-tone spatial dependencies were extracted on each channel.The measures were Angular Second Moment, Contrast, Correlation, Sum of Squares, Inverse Difference Moment, Sum Average, Sum Variance, Sum Entropy, Entropy, Difference Variance and Difference Entropy (Haralick et al. 1973).These measures were extracted only on the original channels, because in a preliminary examination the resolution of 3 m was found to be too coarse.Based on the results of Tuominen and Pekkarinen (2005), the number of requantification classes was set to 16 and the lag value to 3 m.The textural features were calculated omnidirectionally.
Finally, it was tested whether the results of change detection could be improved with old inventory data from 2002.Standwise soil and

Change Detection
The field data were randomly divided into training and test data (Table 1).Due to the large number of potential predictors, i.e. extracted features from the original images and the image conversions on three channels in 2001 and 2004, stepwise discriminant analysis was used in change classification (Tabachnick andFidell 1989, SPSS 2005).All the calculated features in the data sets (Table 3) were included when the classification functions were estimated for the three change classes in the training data.The prior probabilities were set to 1/3 for each change class.The estimated classification functions were subsequently applied to the test data.
The accuracy of the classification results was evaluated by means of confusion matrices (e.g.Campbell 1987).Overall accuracy (OA), producer's accuracy (PA) and user's accuracy (UA) were calculated (e.g.Stehman and Czaplewski 1998).The OA gives too optimistic results if the proportion of one class is high (Campbell 1987); thus the k (kappa) coefficient was calculated (Rosenfeld and Fitzpatrick-Lins 1986).The statistical significance between k coefficients was tested as in Congalton and Mead (1983), and the p-values were corrected by Bonferron correction.

Effect of Stand Properties on the Change Detection
The factors that affect the change detection were tested by logistic regression.Because there were many variables, and their possible interactions and nonlinear effects were also of interest, Classification and Regression Trees method (CART, Breiman et al. 1984) was utilized first.Interesting factors included variables describing the size of the stand, position of the stand in the image, soil type (mineral soil / peatland), the density and size of the growing stock, and tree species composition.Secondly, the statistical significance of the factors and their interactions revealed by CART were tested by logistic regression analysis.

Results
Tables 4 and 5 show the classification results with the different data sets used in classification.
When the k values of data sets Orig* and R3m* are compared, it can be noted that the resolution of 0.5 m performed equally well or slightly better than the resolution of 3 m.The differences were not, however, statistically significant.The radiometric correction did not significantly affect the accuracy of classification (compare data sets R3m*, Ratio* and Emp* in Tables 4 and 5).Adding features to the classification procedure resulted in higher accuracies.For example, with the Orig* data sets, the k values were 0.43, 0.54, 0.56 and 0.58 when median values, histogram features, texture, and both histogram features and texture, respectively, were used in classification.The k value of data set OrigHis was significantly higher than the k of data set OrigMd (p < 0.05).Unexpectedly, the use of old inventory data lowered the accuracy, but the difference was not statistically significant.
With the best data set, OrigHisTex, only 17 features of a total of the 93 possible features were selected in the stepwise discriminant analysis.The selected features were the spectral medians on all the channels, rei on R and G, Inverse Difference Moment on R in 2001and on G in 2001and 2004, Sum and Difference Entropy on G in 2001, Entropy, Contrast and Sum Variance on R in 2004and Sum Average on G in 2004.When both the Moderate-and Considerablechange classes were considered, the comission error varied between 14.3% (RatioHis) and 24.9% (EmpMd) and the omission error between 15.5% (EmpHis) and 26.1% (OrigMd).In the Moderatechange class, 84.2% of the thinned stands were correctly classified (OrigHisTex) and 15.8% were classified into the No-change class.
Because change detection was not problematic in the Considerable-change class, the effects of  In the No-change class, peatland stands were more often classified as unchanged than were stands on mineral soils (Fig. 1).The same holds true for stands where the volume of birch is high.There was a tendency to missclassify unchanged tall stands far from the nadir as well as frontlighted stands far from the nadir as changed.
In the Moderate-change class, the classification failed most often in very sparse and very well-stocked stands (Fig. 2).Dense stands near the nadir and back-lighted dense stands were especially difficult to classify.Some types of operation were so rare that the use of asymptotic statistical tests would have been questionable.The proportion of correctly classified stands by operation type is, however, illustrated in Fig. 3.

Discussion
A semi-automatic method for change detection in boreal forests using bi-temporal aerial photographs was developed.The photographs were   taken as close to each other as possible with respect to time, date and location.Furthermore, the effects of image resolution, radiometric correction, extracted features and stand properties on the change detection were tested.The radiometric correction did not improve the accuracy of classification compared to the accuracy with the original photographs.The obvious reason for this was that the aerial photographs were acquired with almost the same image specifications.Thus, the same unchanged stand looks quite similar on different aerial photographs without any image correction or calibration.Furthermore, in the empiric correction chain different types of vegetation should be corrected separately because the correction chain is sensitive to vegetation type.In this study the whole image area was corrected with the information of the same reference areas.However, these areas were selected from pine-dominated stands where the main interest was focused.
The accuracy of classification was in general at the same level as in some earlier studies (e.g.Saksa et al. 2003, Haapanen and Pekkarinen 2000, Heikkonen and Varjo 2005), but also some improvement was achieved.The considerable changes were classified with high accuracy and, with the best data sets, without error.Meanwhile, the Moderate-change and No-change classes were often confused.With the best data set, however, 84.2% of thinned stands in the Moderate-change class were classified correctly.Furthermore, the omission error in the Moderate-change class was quite low.For example, Hyppänen (1999) found no thinnings and Heikkonen and Varjo (2005) found only 40% of the thinned stands.Thus, the high spatial resolution of aerial photographs clearly helps to detect moderate changes compared to medium-resolution satellite images (cf.Heikkonen and Varjo 2005).
On closer visual inspection, it was found that the visible changes on the images of the undetected stands were mostly very small.For example, only few trees were missing in the later image or only a small part of a stand was thinned.The visual inspection revealed the following common reasons for overclassifying: 1) A road or a clear cutting on the border of an unchanged stand combined with the displacement due to the fact that the images were not registrated.2) In very sparse stands the shadows of the trees were clearly visible on the ground.A minor change in the direction of sunlight caused locally large changes in the DNs when the shadows moved slightly.3) There were still some stands with slight changes that were classified as No-change.Depending on the interpreter, these stands could be classified as No-change or Moderate-change in visual checking because it is hard to make clear line between these classes.Thus, the commission error was actually smaller than shown in results.On the other hand, the large number of unchanged stands classified as moderate changes can make the method unsuitable for practical forestry.However, the unreliably classified stands could be checked visually from bi-temporal images, if some reliability-index was associated to the stands in the classification procedure.
Although some of the stands were very small (0.1 ha), the area of a stand was not found to affect the change detection.In the No-change class the accuracy of classification declined as the distance to nadir increased.This can be explained by the fact that far from nadir even little change in the sun angle causes relatively large change to the positions of tree shadows, which reflects to stand level features.In the Moderate-change class the effect of the distance to nadir was not as clear as in the No-change class which can be seen by comparing Figs. 1 and 2. According to the results in the No-change class, however, it might be wiser to restrict the change detection only to the central parts of the images.
Surprisingly, in the No-change class, the classification results were better on peatlands than on mineral soils.Due, e.g., to changes in water-level, in previous studies the accuracy of change detection on peatlands has been lower than on mineral soils (e.g.Varjo 1996).In the Moderate-change class, soil type had no effect.The effect of soil type was also examined in the training data before the final classifications of the test data.When the mineral soils and the peatlands were classified separately, the accuracy was lower than when they were classified together.The accuracy of classification in this examination was higher in the peatlands than in the mineral soils, as it was also in the final classifications.
The influences of basal area and mean height on the accuracy were significant, but no good explanation for this can be given.The classifica-tion probably fails most often in very sparse and dense stands and in tall stands far from the nadir.In the sparse stands, as noted above, one possible reason is the shadows.High volume of birch increased the probability of correct classification in the No-change class.It was anticipated that this would happen in the Moderate-change class, because deciduous species are often removed in thinnings, but there was no evidence of that.Normal thinnings were quite well recognized compared to other moderate operations and slight storm damage.
The aerial photographs were not temporally registrated, which caused displacement of objects between the ortho-rectified images.The differences in locations varied from 0 to 7 m.At stand level, this was not considered to be a major problem.However, to improve detection of moderate changes, stand-level examination is not adequate.For sub-stand or even tree-level studies the images need to be accurately registered.
In practice, obtaining aerial photographs with the same spatial specification would not be a problem.In Finland, however, the weather conditions limit the number of suitable days for aerial photography.This may complicate obtaining photographs with similar temporal specifications and thus limit the operational use of the method.In this study, photographs with similar specifications were obtained, which enabled exploration of how reliable change detection is by using the bi-temporal aerial photographs in nearly optimal conditions.
There are some issues that need further study.Firstly, the change detection should be carried out with temporally registrated images.Secondly, one image-pair may not be enough to cover the area of interest.Mosaicking of images introduces error, which must be taken into account in further studies.Thirdly, another source of error would be the practical training data from operational databases which, as noted in this study, contain incorrect operations.Finally, the effect of different methods of classification on the accuracy should be studied.In the study of Heikkonen and Varjo (2005), nonparametric classification methods worked better than the parametric method tested.With a maximum likelihood classifier there was a strong tendency to overclassify the Moderate-change class, as was also the case in this study.
stand properties were analyzed in the No-change and Moderate-change classes.Based on CART, potential factors affecting the change detection could be identified.These factors were used as covariates in stepwise logistic regression.Statistically significant (p < 0.05) covariates are depicted in Figs. 1 and 2.

Fig. 1 .
Fig. 1.Stand properties affecting change detection with data set OrigHisTex in the No-change class.The predicted probability of a stand to be correctly classified as unchanged is given on the vertical axis.On each sub-figure, smoothed curves clarify the effect.

Fig. 2 .
Fig. 2. Stand properties affecting change detection with data set OrigHisTex in the Moderate-change class.The predicted probability of a stand to be correctly classified as changed is given on the vertical axis.On each sub-figure, smoothed curves clarify the effect.

Fig. 3 .
Fig. 3. Proportion of correctly classified stands by operation type with data set OrigHisTex.Operation types: 1 = thinning or seed tree felling, 2 = removal of hold-over trees, 3 = tending of young stand, 4 = tending of seedling stand, 5 = slight storm damage, 6 = other operation (soil preparation, partly operated stand, forest road building).The number of stands is given above each bar.

Table 1 .
Main forest characteristics of the training and test data.

Table 2 .
The metadata of the aerial photographs.

Table 3 .
The 11 data sets used in the study.The resolution of the original channels was 0.5 m and that of the other channels 3 m.

Table 4 .
Confusion matrices (3 × 3) derived from 15 classifications of the test data.Each row corresponds to the ground truth and each column to the classification results.The first class is No-change (952 stands), the second Moderatechange (163 stands) and the third Considerable-change (44 stands).The data sets are explained in Table3.

Table 5 .
k, overall, producer's and user's accuracies (OA, PA and UA, respectively) in the test data.OA, PA and UA are in percentages.C1 means No-change, C2 Moderate-change and C3 Considerable-change.The data sets are explained in Table3.