3.1.

Data Collection and Preprocessing

Black-and-white aerial photographs with a paper size of $18\times 18{\mathrm{cm}}^2$ each, which were acquired in August 1962 and August 1981, were used in this research. The 1962 photographs have a scale of $1:6,000$, and were acquired at a flight altitude of 6,450 m. The 1981 photographs have a scale of $1:8,000$, and were acquired at a flight altitude of 5,575 m. All aerial photographs were scanned and the scanned products have a pixel size of 0.5 m. The processing of aerial photographs includes inner orientation, external orientation, linage control points between aerial photographs, orthorectification, and image output. Ground control points were collected from topographic maps with a scale of $1:25,000$ and the 2006 QuickBird image. The root mean square error was less than 1 pixel. Figure 2 illustrates the orthorectified aerial photographs in 1962 and 1981. A detailed description of aerial photograph processing is provided in Dai’s dissertation.

Fig. 2

The 1962 (a) and 1981 (b) aerial photographs after orthorectification.

An orthorectified QuickBird image was acquired in August 2006. QuickBird imagery has four multispectral bands covering visible and near-infrared wavelengths with 2.4 m spatial resolution and one panchromatic band with 0.6 m spatial resolution. In order to make full use of the multispectral bands and panchromatic band with their different spectral and spatial resolutions, principal component analysis was used to integrate both multispectral and panchromatic data into a new dataset with a spatial resolution of 0.5 m for matching the pixel size of aerial photographs.

The stereo photographs were used to produce the 10 m interval contours using MapMatrix software. The contour data were then interpolated and output to produce digital elevation model (DEM) data with a cell size of 2 m.²¹ The MapMatrix software package is a comprehensive platform for processing multiple source data to generate DEM data with different spatial resolutions and digital orthorectified images based on aerial photographs and/or satellite images.²¹ Through automatic control and manual editing of the 10 m contour data, elevation error can be controlled to less than 15 cm.

3.2.

Vegetation Classification

Before performing a land cover classification, it is necessary to design a suitable classification system based on the characteristics of the study area and research objective. A classification system with three land-cover types–Schrenkiana, Sabina, and others were designed for this research. The Schrenkiana class is defined as a single tree or forest having tree height of over 2 m and crown size of greater than $1.5{\mathrm{m}}^2$. The Sabina class is defined as the vegetation having coverage of greater than 10% within an area of $20{\mathrm{m}}^2$ and a height of greater than 0.5 m. The others class is any other land covers except Schrenkiana and Sabina.

Per-pixel based classification algorithms are the most common methods for land cover classification,²² but often generate relatively poor classifications when very high spatial resolution remote sensing data are used.²³ Utilizing the spatial features inherent in high spatial resolution images, object-oriented classification algorithms have been proven to provide better classification than the per-pixel based algorithms^23,24 and thus the object-oriented classification using the Definiens eCognition software²⁵ was used in this research. The object-oriented classification method is based on two critical steps—image segmentation and classification. A good segmentation result is fundamental for this method. One critical step is to identify suitable parameters to generate the segment image. The determination of segmental scale is required, which is often based on multiple scales, depending on the complexity of landscape and spatial resolution of remote sensing image. In general, some land covers such as bare soil, desert, and dense crops are composed of relatively large, homogeneous patches, and thus a relatively large segment scale is required; while other land covers such as residential and grass in urban landscape have relatively small patch sizes and complex boundaries, thus a relatively small segment scale is necessary. In addition to the use of scale parameter, other parameters—color, shape, smoothness, and compactness—must also be evaluated and defined. These parameters are used to determine the fitness between the segmentation boundary and land cover boundaries, and they are modified through weights of each parameter. The procedure to identify the optimal parameters is often dependent on the complexity of the study area and the nature of the remotely sensing data being classified. After segmentation, the classification is conducted on the segment images. Different supervised classification algorithms such as k-nearest neighbor and decision tree can be used, which are based on training samples selected from the segment images respectively. A detailed description of the object-oriented classification is provided in the user guide of eCognition Developer. In this research, we used this software to develop vegetation classifications using the 2006 QuickBird image as well as scanned aerial photographs, based on the optimization of the parameters used in eCognition for segmentation and the k-nearest neighbor algorithm for vegetation classification.

In the 1981 photographs, most of the Sabina class can be accurately classified, but in the 1962 photographs, Sabina cannot be separated from other land cover due to the photograph quality. Thus, only two classes—Schrenkiana and others were classified in the 1962 aerial photographs. Compared with the QuickBird image, the aerial photographs have relatively poor data quality, thus further post-processing was required. For example, the Schrenkiana canopy-cast shadows are confused with grass in the black-and-white photographs. This problem is solved by establishing the relationship between canopy and neighboring tree cast shadows. The relatively large shadows in the forest edge often result in forest misclassification in segmentation process due to the impacts of terrain slope and sun elevation angle. This problem is solved using the standard deviation of pixels after segmentation.

Field work was conducted in September 2010 and September 2011. Ninety-two plots in the reachable slope areas were documented. Each plot was a circle with a semi-diameter of 5 m. The area percent of Schrenkiana and Sabina were measured within each plot. These plots were used to evaluate the classification result from the 2006 QuickBird image. Although the date for ground truth data collection has a 4 to 5 years difference with the QuickBird image acquisition date, the vegetation type can be assumed unchanged in this study area due to the low-level of human activity and very slow growth rate at the upper area of the mountainous region. For the results from the 1962 and 1981 aerial photographs, since no field measurements were available, the test samples were randomly collected by visual interpretation from the aerial photographs. In this study, a total of 90 plots were randomly generated from each date of photographs, with a circle of 5-m semi-diameter for measuring the coverage of each class. For example, Schrenkiana, Sabina and others have coverages of 30%, 20%, and 50%. The error matrix approach was used to assess the classification result for each date. Overall classification accuracy and kappa coefficient for each classification image, and producer’s accuracy and user’s accuracy for each land cover type were calculated from each error matrix.^26,27

Topographic factors, such as elevation, slope, and aspect are important factors influencing sun illumination, precipitation, and temperature, as well as soil type and depth in mountainous regions. These factors comprehensively affect vegetation species distribution and vegetation growth rates.^6,28,29 In this research, the analysis of vegetation distribution and topographic factors was based on the 2006 QuickBird classification image and the DEM data. Considering the data availability, a total area of $33.29{\mathrm{km}}^2$ [see subregions 1 and 2 in Fig. 1(b)], which coincided with both the QuickBird coverage and DEM data, was selected. The elevation range of approximately 2,300 to 3,300 m in this study area is a typical site for examining tree line change in an arid ecosystem.

3.3.

Vegetation Change Detection

The vegetation change analysis was based on the overlap areas of aerial photographs in 1961 and 1981, the 2006 QuickBird image, and the 2-m DEM data. Two subregions as illustrated in Fig. 1 were selected to examine the impacts of topographic factors on vegetation expansion. Considering the vegetation classification results, the analysis of Schrenkiana change was based on both periods from 1962 to 1981 and 1981 to 2006 for both subregions, and the analysis of Sabina change was only based on the period from 1981 to 2006 because of the difficulty in separating Sabina from other land covers based on the 1961 aerial photographs.

The postclassification comparison approach is often used to examine land cover change trajectories and this method was used to examine vegetation expansion in this research. The vegetation change result was then overlaid on the topographic layers (e.g., elevation, slope, and aspect) to examine how different terrain conditions affect vegetation change, especially the dynamic change of tree lines over the past four decades.

4.1.

Analysis of the Relationship between Vegetation Distribution and Topographic Factors

Classification accuracies for the three dates are summarized in Table 1. The 2006 QuickBird image produced a higher overall classification accuracy than the 1961 and 1982 aerial photographs due to its high-quality multispectral image compared with only one panchromatic band in the aerial photograph. The relatively high accuracy result in 1962 compared with the result in 1981 may be due to the different numbers of classes, that is, only two classes in 1962 compared with three classes in 1981. As an example, Fig. 3 shows how vegetation distribution relates to elevation by overlaying contour lines on the 2006 QuickBird vegetation classification image. It indicates that Schrenkiana is mainly distributed in the regions with elevation values of less than 2,600 m, and most of Sabina is distributed in the regions with elevation of less than 2,900 m but greater than 2,600 m. This implies that elevation may be an important factor influencing vegetation distribution.

Table 1

Accuracy assessment results for the three dates of classification images.

Fig. 3

Vegetation distribution within the study area by overlaying contour lines.

The statistical analysis of vegetation distribution along elevation and slope ranges can better show the relationship of vegetation distribution and topographic factors (see Fig. 4). In this case study, Schrenkiana is mainly distributed within the elevation range of 2,300 to 2,650 m with a peak value of 2,400 m. The Schrenkiana area decreases gradually from the maximum value at the elevation of 2,400 m to the minimum value at the elevation value of 2,850 m. Sabina is mainly distributed in the elevation range of 2,600 to 2,900 m with the peak at 2,750 to 2,800 m. The Sabina area decreases gradually from the peak at 2,800 m to the minimum at the elevation of 3,250 m [see Fig. 4(a)]. In this study area, Schrenkiana is the dominant vegetation in the elevation range between 2,300 and 2,650 m; Schrenkiana and Sabina become evenly distributed in the elevation range between 2,650 and 2,700 m; and Sabina areas become dominant when the elevation value is greater than 2,700 m. However, after the elevation reaches over 3,250 m, grass and bare lands become dominant. The relationship between vegetation distribution and elevation fit the alpine vegetation distribution well.

Fig. 4

Relationship of vegetation areas with elevation (a) and slope (b).

Both Schrenkiana and Sabina have similar relationships to slope, and are mainly distributed within the slopes between 25 and 36 deg [Fig. 4(b)]. However, Schrenkiana area reached a peak at slopes of about 30 deg, and Sabina reached a peak at slopes of about 33 deg. This is because Sabina is often distributed in higher elevation areas with steeper slopes than Schrenkiana [Fig. 4(a)].

The statistical analysis of vegetation distribution along slope aspects indicates that Schrenkiana and Sabina distribution varied significantly in aspects [Fig. 5]. Schrenkiana is mainly distributed along the northern aspect (about 300 to 60 deg), and its area gradually decreases from both northwest (from about 300 to 180 deg) and northeast to south (from about 60 to 180 deg) [Fig. 5(a)]. On the other hand, Sabina distribution is mainly distributed in the southwest (about 200 to 300 deg) and northeast (about 40 to 100 deg) and had limited distribution in the north, northwest and northeast (about 300 to 40 deg). The study area has very limited area distribution in the southern aspect (about 140 to 190 deg). Figure 5 also indicates that Schrenkiana is mainly distributed in the areas away from the sun, and Sabina is mainly distributed in areas facing the Sun. This is related to the physiological property of both Schrenkiana and Sabina. Schrenkiana is a shade-loving plant at young stages and sensitive to aspects, while Sabina is sunny-loving plant, less sensitive to aspects, and has a strong capability in different biophysical environments, such as dry and poor soil conditions, thus Sabina is more commonly distributed in the aspects facing the Sun than Schrenkiana.

Fig. 5

Relationship between vegetation distribution and aspects displayed by using (a) scatter plots and (b) 360-deg circles.

4.2.

Analysis of the Relationships between Vegetation Change and Topographic Factors

Analysis of Schrenkiana change at both subregions as illustrated in Fig. 6 indicates that the increased Schrenkiana areas are mainly distributed at the forest frontiers, especially in the relatively high elevation areas, indicating the tree line change. Table 2 indicates that the increased areas in both subregions are higher in the period from 1981 to 2006 than in the period from 1962 to 1981. Examination of the changed Schrenkiana areas and corresponding elevation data showed that the elevation of Schrenkiana distribution increased by 3 m in subregion 1, and increased by 2 m in subregion 2 between 1962 and 2006. The increased elevation of Schrenkiana distribution in the past 44 years was not much due to its slow growth rate. If the Schrenkiana distribution in 1962 is used as a base line, the increased areas vary at different elevation ranges. At subregion 1, the increasing trend is similar in both periods of 1981 to 2006 and 1962 to 1981 (see Fig. 7), and the peak is in elevation value of 2,650 m. When the elevation falls between 2,300 and 2,500 m, Schrenkiana has the same increase area values in 1981 and 2006 [see Fig. 7(a)], implying that climate has little impact on vegetation change at this region. However, when elevation reaches over 2,650 m, the increased Schrenkiana area in 2006 was obviously higher than in 1981, implying the increasing trend of Schrenkiana area in higher elevations. In the subregion 2, increased Schrenkiana area is obvious in higher elevation areas. The peak of increased area in 2006 has a higher elevation value than the peak value in 1981, and the highest elevation value of increased areas in 2006 is also higher than in 1981 [see Fig. 7(b)], implying that Schrenkiana expanded to higher elevation areas.

Fig. 6

Distribution of Schrenkiana dynamic change in subregion 1 (a: 1962 to 1981; b: 1981 to 2006) and subregion 2 (c: 1962 to 1981; d: 1981 to 2006).

Table 2

A summary of vegetation areas among three dates between two subregions.

Fig. 7

Relationships between increased Schrenkiana areas and elevation for (a) subregion 1 and (b) subregion 2.

In subregion 1, Schrenkiana expansion areas within the slopes of 25 to 40 deg account for 77.5% of the total expansion areas and have similar slopes to the peak of expansion in 2006 and 1981, but the peak value in 2006 is much higher than that in 1981, as illustrated in Fig. 8(a). Another relatively small peak of expansion areas is in slopes of 10 to 15 deg. In subregion 2, Schrenkiana distribution has two peaks in 1962, and the area expansion in 1981 has a similar trend as the year of 1962, but only one peak in 2006 [see Fig. 8(b)]. The expansion peaks in 2006 and 1981 are in the slope ranges of 25 to 40 deg, accounting for 74.5% of the total expansion areas. Although Schrenkiana distribution has different slope features in both subregions, the expansion areas falls in the same slope ranges, i.e., 25 to 40 deg, implying that topography is an important factor influencing Schrenkiana distribution and expansion.

Fig. 8

Relationship between Schrenkiana expansion areas and slopes for (a) subregion 1 and (b) subregion 2.

Considering slope aspects, the expansion trends in 2006 and 1981 are similar, as illustrated in Fig. 9. In subregion 1, Schrenkiana expansion in both periods varies in aspects compared with the Schrenkiana distribution in 1962. The peaks of the 1962 Schrenkiana distribution fall in the slope aspects of 0 to 45 and 270 to 315 deg, respectively [see Fig. 9(a)], while the expansion peak is also included in the slope aspect of 45 to 90, in addition to 270 to 315 deg. In subregion 2, most of the Schrenkiana expansion is distributed in semi-northward and northward slopes [see Fig. 9(b)], but has some in sunny-slopes, especially in 2006 with less than 3% of the total expansion areas.

Fig. 9

The relationship of Schrenkiana expansion areas and slope aspects for (a) subregion 1 and (b) subregion 2.

The distribution of Sabina expansion areas between 1981 and 2006 can be explored by overlaying contour lines (Fig. 10) and indicates that the increased Sabina areas in both subregions are distributed at high elevation areas. Analysis of changed Sabina area and the corresponding elevation data shows that the elevation of Sabina distribution increased by 7 m in subregion 1 and increased by 12 m in subregion 2. The increased Sabina areas are mainly distributed along the edges of original shrubs, especially in subregion 2, although limited areas increased in the inner open space.

Fig. 10

Distribution of Sabina expansion areas between 1981 and 2006 for (a) subregion 1 and (b) subregion 2.

Sabina expansion areas along with elevation change varied considerably between the subregions, as illustrated in Fig. 11. In subregion 1, the peak of Sabina expansion area varies with the peak of the 1981 Sabina distribution, and the elevation of peak becomes smaller [Fig. 11(a)]. The Sabina distribution in 1981 has only one peak, but the expansion between 1981 and 2006 has three peaks with the peak values become smaller as elevation becomes higher. In subregion 2, the peak of Sabina expansion area is similar to the 1981 Sabina distribution [Fig. 11(b)]. However, the highest elevation of Sabina expansion in subregion 2 is obviously lower than that in subregion 1; in contrast, the Sabina expansion in subregion 2 at lower elevation areas is obviously higher than in higher elevation areas, and the expansion trend in the high elevation area in subregion 2 is smaller than in subregion 1, implying that the Sabina expansion in subregion 1 is more prone to occur at higher elevation areas.

Fig. 11

Relationship of Sabina expansion areas and elevation for (a) subregion 1 and (b) subregion 2.