THE STUDY OF LAND COVER CHANGE USING CHANGE VECTOR APPROACH INTEGRATED WITH UNSUPERVISED CLASSIFICATION METHOD: A CASE IN DUY TIEN (VIETNAM)

. Investigating information on land cover changes is an indispensable task in studies related to the variation of the environment. Land cover changes can be monitored using multi-temporal satellite images at different scales. The commonly used method is the post-classification change detection which can figure out the replacement of a land cover by the others. However, the magnitude and dimension of the changes are not been always exploited. This study employs the mixture of categorical and radiometric change methods to investigate the relations between land cover classes and the change magnitude, the change direction of land covers. Applying the Change Vector Analysis (CVA) method and unsupervised classification for two Landsat images acquired at the same day of years in 2000 and in 2017 in Duy Tien district, the experimental results show that a low magnitude of change occurs in the largest area of direction I and direction IV regarding the increase of Normalized Difference Vegetation Index (NDVI), but the opposite trend of (Bare soil Index) BI in the rice field. Alternately, the high magnitude of change is seen in the build-up class which occupies the smallest area with 1700 ha. The characterized changes produced by the CVA method provide a picture of change dynamics of land cover over the period of 2000-2017 in the study area. The Study Of Land Cover Change Using Change Vector Approach Integrated With Unsupervised Classification Method: A Case In Duy Tien (Vietnam). Geography, 2, p. 175-184 Conflict of The


INTRODUCTION
Land cover, the layer of physical material at the earth's surface, is composed of water, trees, grass, bare ground, etc. Land cover change is one of the most important domains in environmental change research. Investigating information on land cover changes is an indispensable task in studies related to the variation of the environment (Chen et al. 2003;Dewi et al. 2017;limanova et al. 2017). Change detection using satellite images has been proved to be an efficient approach with the advantages of rapidly collecting information, and multitemporal, multi spatial resolution, which allows us to map land cover changes at different scales (Lambin and Strahlerf 1994;Chen et al. 2012). Various techniques for change detection using remote sensing data have been developed and broadly applied. These techniques can be separated into two groups: post-classification comparison or two-date classification comparison; and temporal radiometric change analysis (Ding et al. 1998;Chen et al. 2012). The post-classification comparison is based on the classification results to determine the variations, which can produce a detailed change matrix of land cover classes. However, this group method does not allow detecting subtle changes in the classes, and the mechanism of changes is not discovered. Alternatively, a temporal radiometric change method analyses the variation curve or trajectory of spectral reflectances achieved by remote sensing sensors for successive times. Several methods using temporal radiometric changes have been developed such as rationing, band differencing, vegetation indices, principal component analysis, regression analysis, and  (Chen et al. 2003). Similarly, to other temporal radiometric change techniques, the CVA one hand can minimize the errors accumulating from the individual classification of two or more satellite images, and it measures not only the intensity of change but also the trend of the change for each pixel on multitemporal satellite images (Dewi et al. 2017). In addition, the CVA method can process any number of spectral bands to generate a single band which represents the intensity of change between two or more satellite images (Molina et al. 2012). For these reasons, the CVA is a useful technique to investigate the information of land cover changes for various purposes as: detecting land-use/ landcover change in rural-urban fringe areas (He et al. 2011); assessment of wetland dynamic (Landmann et al. 2013); monitoring changes in fuzzy shorelines (Dewi et al. 2017); detecting forest changes (Malila 1980). On the other hand, applying the CVA technique may face some challenges composed of: the reliable radiometry of images over a period of time, the difficulties of effectively detecting thresholds for separating the change or no-change of the magnitude pixels, and the problem of discriminating phenomenological types of the changes due to the use of a great number of bands for analysis (Cheng et al. 2003). These challenges limit the potential of such an effective technique as CVA for the land cover change monitoring.
This study aims to overcome these challenges and applies the CVA technique to investigate the land cover change of a rural area at a local scale. The CVA method produces the direction and the magnitude of changes but the categories of land cover. Implementing the CVA method in a case study area of Duy Tien district, Vietnam, we propose a mixture approach of the temporal radiometric change and the categorical information of land covers to assess the effectiveness of the method and to figure out the relations between land cover classes and their change components.

Study area
The study area is in Duy Tien district of Ha Nam province, Vietnam, located in the center of the Red River delta with an area of 13,765 hectares composed of 19 communes and 2 towns (Fig. 1). Duy Tien was a purely agricultural district with the dominant land cover of the rice field. Recently, economic development required the exchange of a huge area of the agricultural land to industrial zones, dependent services, and residence expansion. Moreover, the very low valuation of rice cultivation led to the replacement of the rice fields with perennial trees or farms. The rapid change of land cover makes the need for monitoring not only the types of changes but also the direction of changes for environmental management of local government.

Date in use
To monitor the land cover changes of the study area, a Landsat 7 ETM+ satellite image acquired on 17/09/2000 and a Landsat 8 OLI image acquired on 17/09/2017 are used for analysis ( Table 1). The specifications of the Landsat 7ETM+ and the Landsat 8OLI satellite image are illustrated in Table 2. The two images are the no-cloud cover, and they are processed at the Landsat Surface Reflectance High-Level Data products (level L1TP) by the United States Geological Survey (USGS). The effects of the atmosphere and the distortion of the images related to topography are minimized at this processed level (US Geological Survey 2016). All images have been systematically imageto-image registered with the Root Mean Square Error (RMSE) less than 12m, and they are adjusted to the 30 m spatial resolution (30 m x 30 m pixel size) in the Universal Transverse Mercator (UTM) projection zone 48N (Zanter 2017). The images are acquired at the late-season stage (the ripening phase) of rice which presents identical color and tone. These two images were acquired on the same day of years 2000 and 2017 at similar local time, similar solar zenith and solar azimuth angle; hence, this provides very good conditions to apply the CVA method. Although the images were acquired in different paths/rows, the study area is fully covered by both scenes.

Methodology
CVA method. Ideally, the CVA method is described as a change vector with two components: a vector direction (the angle of the change) and a magnitude of the change between two or more images acquired in two dates (Malila 1980). The surface reflectance of an object measured on satellite images acquired in two dates corresponding to T1 and T2 are given by matrices T1 = (p 11 , p 12 , p 13 ,... p 1n ) and T2 = ( p 21 , p 22 , p 23 ,... p 2n ). Where p,n are the surface reflectance and the number of the bands, respectively. The change vector ∆ is presented as: The magnitude of changes (I∆I) between two dates are calculated from all the spectral reflectance change of the objects achieved on bands.
The direction of change can be investigated by calculating the angle of the change vector. The angle of the change vector represents the types of land cover change over a period of time. The individual band of satellite images can be replaced by a band index or a band ratio to calculate the magnitude and the type of changes.
In this study, we use the Normalized Difference Vegetation Index (NDVI) and the Bare soil Index (BI) calculated from Landsat images to quantify the change of the land covers in Duy Tien. The NDVI displays the relationship between the quantity of chlorophyll in leaves and the spectral reflectance in the wavelength of red and near-infrared, thus the NDVI image is popularly used to research vegetation as estimating biomass, plant productivity, fractional vegetation cover (Rouse 1974;Richardson 1977). The BI is used to distinguish an agricultural land and non-agricultural land (Jamalabad 2004). The NDVI and BI indexes are calculated by equations as: Where p Blue , p Red , p Nir , p Swir are the surface reflectance of the blue, red, near-infrared, and short wave infrared regions, respectively.
The CVA method is applied for the Landsat images captured in 2000 and 2017 to monitor the magnitude and the direction of changes using two components; the NDVI index and the BI index (Fig. 2). The direction of changes is clearly classified into four types of changes as type I, type II, type III and type IV corresponding to four possible corners of the change vector (Fig. 2B). The magnitude of change is then stratified into low change, medium change and high change using empirical thresholds detected by investigating the histograms of the image of the change magnitude combined with the knowledge of the study area (Fig. 3). As an example, the areas with rice no-change are checked with the values of the change magnitude image to define a threshold for the low change layer. The variations from rice fields into the other vegetation or infrastructures are used to define thresholds on the change magnitude image for the medium change or the high change layers, respectively (Fig. 3A, 3B).
The Landsat 8 OLI image captured on 17/9/2017 is performed classification to achieve a land cover map. To do the classification, we used the ISODATA method to cluster the images into 40 classes from running 20 iterations of all the bands of the Landsat 8 image. The unsupervised classification method is used in this study to maximize the possibility of distinguishing clusters based on spectral reflectance but it minimizes the impact of visual interpolation on results. A post-classification is implemented to combine the previous 40 classes into 8    (Table 3). The land cover map is overlaid with the change magnitude and the change direction maps to analyze the relationship between them. Field visiting was implemented on 14th September 2018 to observe the land covers in 47 locations in the study area. The information on field points is used to assess the accuracy of the land cover map.

RESULTS
The change magnitude and the change direction. The two dimensional scatter plots in Fig. 4A and Fig. 4B represent the spatial relation of NDVI and BI in 2000 and 2017. NDVI and BI are two dimensions in the change vector coordinate. The trend of data points (red line) in Fig.  4 shows a negative correlation between NDVI and BI. To apply successfully the CVA method, it must be sure that the variation of the vegetation index is opposite to the variation of the bare soil index. Thus, the more negative correlation between them is, the higher ability of change determination is. The negative correlation between the two indexes approves the suitability of using NDVI and BI for detecting the change vectors which produce the reliable maps of the change magnitude and the change dimension for further analysis.
As seen in Fig. 5, the change direction is discriminated by four types (so-called Direction) of the changes over

2000-2017.
Type I shows the increase in both of the NDVI and BI index, which covers only the area of rice field predominantly distributing in the study area. Type II represents the area with the increase of BI but the decrease of NDVI. This type of change is obviously seen in the changing area from the agricultural land to industrial zones or some roads as well. The degradation in both NDVI and BI values of Type III is found in the surrounding water surface such as the river of fishponds. The change in Type III may relate to the moisture variation of the land cover. The largest area of the change direction belongs to Type 4 defined by the increase of NDVI and the deduction of BI. Type 4 occupies all of the water surface and a vast area of vegetation. The distribution of Type 4 is complicated, thus, it needs to discuss more details when integrating with the other data such as the land cover map and the magnitude change map. Fig. 6 illustrates the magnitude of the land cover changes at three levels of low change, medium change, and high change corresponding to three threshold ranges as 0.0004-0.1428, 0.1428-0.2709, 0.2709-1.235 of the change magnitude. These thresholds are empirically defined by the correlation between the land cover changes and the mean, the standard variation of the histogram of the change magnitude image (showed in Fig. 3). Low magnitude change represents the variation of the spectral reflectance at an insignificant level that is not strong enough for indicating a replacement of a land cover class over the period 2000 -2017. Hence, the low change level covers almost all of the study area. Medium  However, the medium change of the river is an exception due to the variation of turbidity of water at the two dates. High change magnitude shows the replacement of a class by the other classes with different physical characteristics such as an exchange of rice by build-up area or bare soil by vegetation.
The map of land covers and accuracy assessment. Fig. 7 represents the distribution of 8 land cover classes in the study area on 17th September 2017. Although some vast industrial zones have been constructed in Duy Tien, the agricultural land occupies most of the area with 45% of rice, 9% of perennial trees and 6% of annual agriculture. Clear water is fishponds and canals, which are discriminated from turbid water by the density of suspended sediment.
The extremely high turbidity of up to 1 kg of sediment in 1 m3 water of the Red river, as reported by (Hoa 2001) leads to very high reflectance of visible wavelengths, thus the turbid water is clearly seen in the multispectral satellite image. Clear water and turbid water cover 11% and 3% area, respectively. Classes without vegetation as build-up area, rural residence, bare land occupy an amount of 26% area in which the bare land is some sandbars along Red river, cover only 1% of the study area. Table 4 shows the user's accuracy and the producer's accuracy of each class of the land cover map built from Landsat 8 image acquired in September 2017. Totally, 47 field points have been used as the reference data for the accuracy assessment. The highest accuracy is found in the bare land and the turbid water with 100% of the right pixels because of the clear presentation of these features on the satellite image. The lowest accuracy is seen in the annual vegetation class with only 29% of the user's accuracy. This number of accuracy is unexpected but it represents the difficulty to extract such a mixed class as annual vegetation in the agricultural rural area. Above all, the general accuracy of the map is 68%. The land cover map with this accuracy is acceptable in this experiment. The assessment of temporal radiometric changes and the categorical information of land covers. As seen in Fig.  8A, the low magnitude of change occurs at the largest area in the direction I (Type I) and direction IV corresponding to 4200 ha (32% of total area) and 2600 ha (20%). Direction I and direction IV are both increases of NDVI but the opposite trend of BI. Hence, these directions are related to vegetation such as rice with 23% in direction I, 16% in direction IV or perennial trees with 7% in direction IV (Fig.  8B). The medium magnitude of change covers the largest area of direction IV with 11% of the total area. However, this level of change does not emphasize a specific land   cover class. The area of medium change is relatively equal for all the classes (Fig. 8C). Alternately, a high magnitude of change occupies the smallest area with 1700 ha (12% of total area), but it almost covers the build-up area with 3% and 4% of the clear water area. High change can obviously be observed in the satellite images due to the conversion of vegetation to the concrete surface which makes the absolute change of spectral reflectance. An amount of 519 ha (4%) of clear water felt in the high change area may be caused by the appearance of the vanishing of vegetation on the clear water surface. In addition, high change occurs mostly in direction IV area (6%) where the NDVI increases over the period 2000-2017. In sum, the results of the experiment show the correlation of the change magnitude, the change direction, and land cover classes. The characterized changes in land covers exploit the nature of change.

DISCUSSION
The CVA method is based on the radiometric change, which is an efficient approach by exploiting the hidden information of the changes through the change magnitude and the change direction. However, the challenges related to using spectral reflectance and the determination of thresholds of change may impact the accuracy of applying the CVA method. Thus, we discuss more of these issues.
As mentioned by (Chen et al. 2003), the radiometric changes of the images acquired at the two dates can be affected by disturbing factors caused by the temporal variation such as the different atmospheric conditions, solar angle, soil moisture, and vegetation phenology, etc. These quick variation conditions may contribute to the bias of the results more critical than the spectral reflectance change detected by the CVA method. To minimize the error coming from temporal variation, it is necessary to use the pair of satellite images acquired at as much similar as possible conditions. In this study, we use two images acquired on the same day of the year on 17th September in 2000 and 2017. The difference of several seconds of the acquisition times (Table 1) provides a radiometric signal at a similar solar zenith, horizontal angle. In addition, the study area suffers the sub-tropical climate with a cold season, thus, the time of cultivation is relatively stable over the years. The images acquired on the same day of year measures the spectral reflectance from the similar phenological stage of rice. Moreover, these two images are captured by the same family of satellite sensor (Landsat family), the same way of atmospheric correction, the same spectral resolution, and the same range of the accuracy of geometric correction. Thus, the satellite images used in this case are wonderful data to apply the CVA method to achieve the highest accuracy of results.
In the concept of CVA, this method is applicable to any number of spectral bands, even though any measurement scale of radiance is used (Malila 1980). However, the use of multispectral images can make the confusion and the difficulties to discriminate the categories of changes when analyzing the change directions. The more the number of image bands in use, the analysis of change direction is more complicated, especially in the study area with numerous types of land covers. Three methods to face this issue have been developed (Chen et al. 2003): trigonometric function of vector angle in two spectral dimensions (Malila 1980); principal component analysis in multi-temporal space (Lambin and Strahler 1994), and sector coding in more than two spectral dimensions (Virag and Colwell 1987). The use of the sector coding method to separate the change direct into 4 categories in this study is the appropriate option for the study area. Each category represents a type of the change direction corresponding to two correlated components NDVI and BI index. The stratification scheme of the change direction can be used to figure out some causes of the changes which may not be detected by visual interpolation.
One of the most difficult tasks of the CVA application is the determination of appropriate thresholds to distinguish the change-magnitude pixels with change or no-change. However, there is no automatic or semi-automatic method to define the change thresholds because it is necessary to consider some ecological and spectral conditions corresponding to threshold selection and overall change sensitivity. As an example, a lower change threshold value may allow the inclusion of slightly changed wetlands into the change analyses, while a high threshold value may only include the locations of significantly changed areas (Baker et al. 2007). Hence, using the remote sensing analyst's expert knowledge of the study area is the best method to detect the change thresholds (Jano et al. 1998). In this study, we determine the change thresholds based on the mean and standard deviation of the histogram of the change magnitude image combining with the knowledge of land cover changes. The results proved that this approach is effective and it produces reliable maps of land cover changes in such a small known area as in Duy Tien. However, for a large scale area with the complication of land cover classes, the use of the CVA method with this method of threshold determination will face numerous challenges and it needs to be studied in future works.

CONCLUSIONS
The application of CVA using two components NDVI and BI index is an effective method to discover the dynamic of land cover/ land use changes. The CVA method does not only exploit the magnitude element but also the direction element of the changes which represent the nature of land cover changes using remote sensing data. On the other side, the results of the CVA application make it difficult to indicate the specific land covers that can obviously be obtained using a single classification of the satellite image. This study successfully applies both CVA and unsupervised methods to investigate the change information of land cover. The study shows that the CVA method using NDVI and BI index is suitable for Landsat images to monitor land cover change in such a purely agricultural area as Duy Tien district. The key step of the CVA method is the method to detect appropriate thresholds for stratifying pixels with change or no-change, for which the knowledge of the study area combined with the spectral reflectance of surface features is the most accurate approach. The change magnitude or change direction obtained from the CVA application can be disturbed by the temporal variation of the environment, the different sources of satellite images and the radiometric, geometric processing as well. Thus, the used satellite images must be carefully selected and processed to minimize the impacts of those disturbing factors. Overall, the CVA with the ability to measure the magnitude of land cover change is a promising method to monitor the growth of the plants from that we can estimate the production of crops. This application is necessary for an agricultural country like Viet Nam, and it may be employed for future studies.