Land Cover Change Detection Using Multiple Shape Parameters of Spectral and NDVI Curves

: Spectral and NDVI values have been used to calculate the change magnitudes of land cover, but may result in many pseudo-changes because of inter-class variance. Recently, the shape information of spectral or NDVI curves such as direction, angle, gradient, or other mathematical indicators have been used to improve the accuracy of land cover change detection. However, these measurements, in terms of the single shape features, can hardly capture the complete trends of curves affected by the unsynchronized phenology. Therefore, the calculated change magnitudes are indistinct such that changes and no-changes have a low contrast. This problem has prevented traditional change detection methods from achieving a higher accuracy using bi-temporal images or NDVI time series. In this paper, a multiple shape parameters-based change detection method is proposed by combining the spectral correlation operator and the shape features of NDVI temporal curves (phase angle cumulant, baseline cumulant, relative cumulation rate, and zero-crossing rate). The change magnitude is derived by integrating all the inter-annual differences of these shape parameters. The change regions are discriminated by an automated threshold selection method known as histogram concavity analysis. The results showed that the mean differences in the change magnitudes of the proposed method between 2100 changed and 2523 unchanged pixels was 32%, the overall accuracy was approximately 88%, and the kappa coefﬁcient was 0.76. A comparative analysis was conducted with bi-temporal image-based methods and NDVI time series-based methods, and we demonstrate that the proposed method is more effective and robust than traditional methods in achieving high-contrast change magnitudes and accuracy.


Introduction
Land cover dynamics and their spatiotemporal information play a vital role in many scientific fields, such as natural resources management, environmental research, climate modeling, and earth biogeochemistry studies [1][2][3]. In the past decades, a variety of bi-temporal change detection methods have been developed to monitor land cover conversions between two years based onLandsat images or NDVI time series [4][5][6][7]. However, the method of comparing the responses of a single image band or NDVI which is directly accepted by some previous algorithms (e.g., image differencing and image ratioing) may lead to pseudo-changes, which hampers the monitoring process from local to regional scales [8][9][10]. Individual values in the spectrum or temporal curve of land cover type are more vulnerable to environmental factors and thus result in data uncertainty [11,12]. If change magnitudes were calculated based on contaminated values of the spectrum or temporal curves, pseudo-changes

Methodology
The main idea of the multiple shape parameter-based (MSP) change detection method is the utilization of shape features from both the spectrum and the NDVI temporal curve to calculate a high-contrast change magnitude that can mostly discriminate the change region from various change magnitude types. It includes three major components: change magnitude calculation (see Section 2.1), change magnitude integration (see Section 2.2), and change region determination (see Section 2.3). Figure 1 gives the framework of the MSP method. In the change magnitude calculation part, five shape parameters (correlation operator for spectral curves, phase angle, baseline cumulant, relative cumulant, and zero-crossing rate for NDVI temporal curves) are adopted to describe the shape features of the spectrum and the NDVI temporal curve and perform measurements from these different parameters. The differences of bi-temporal multispectral image pairs or NDVI time-series pairs are taken to calculate the change magnitudes independently using these shape parameters. In the change magnitude integration part, the units of different change magnitudes are assimilated by standard normalization. These unified change magnitudes are integrated with a weighting function. Finally, the changed and unchanged regions are acquired by comparing the change magnitude with a threshold.
Remote Sens. 2018, 10, x FOR PEER REVIEW 3 of 20 parameters. For instance, some of the parameters are the beginning of the season, the end of the season, amplitude, and the rate of increase [19]. The phenological variability of land cover could be measured by these metrics [29]. Therefore, using these parameters or metrics to model the shape feature of the NDVI temporal curve is feasible. Considering the different spectrum of the same land cover type, the correlation operator of the spectrum is more appropriate to ensure that the shape information of the spectrum is utilized but not overexploited. The integration of advantages from the spectra and NDVI temporal curve could further reduce the pseudo-changes and allow easier discrimination of real changes from non-changes.
In this study, a multiple shape parameters-based change detection method using one pair of Landsat images from two years and two yearly NDVI time series was developed to calculate the change magnitude with high contrast under the confusions of conflicting differences and improve change detection accuracy. For the NDVI temporal curve, shape parameters were mainly formed by reconstructions of seasonal parameters and the application of zero-crossing rate. The spectral correlation operator was adopted as a special spectral shape parameter. The change magnitudes were derived from differences of all shape parameters and integrated as one final change magnitude. Then, change discrimination was applied based on an auto-threshold method. It mainly consisted of three parts, which were change magnitude calculation, change magnitude integration, and change region discrimination, as described in Section 2. Section 3 presents the experimental results and reports on the accuracy of the proposed method with a comparison regarding the several bi-temporal images or NDVI time series-based methods mentioned above. The discussion and conclusions are provided in Sections 4 and 5, respectively.

Methodology
The main idea of the multiple shape parameter-based (MSP) change detection method is the utilization of shape features from both the spectrum and the NDVI temporal curve to calculate a high-contrast change magnitude that can mostly discriminate the change region from various change magnitude types. It includes three major components: change magnitude calculation (see Section 2.1), change magnitude integration (see Section 2.2), and change region determination (see Section 2.3). Figure 1 gives the framework of the MSP method. In the change magnitude calculation part, five shape parameters (correlation operator for spectral curves, phase angle, baseline cumulant, relative cumulant, and zero-crossing rate for NDVI temporal curves) are adopted to describe the shape features of the spectrum and the NDVI temporal curve and perform measurements from these different parameters. The differences of bi-temporal multispectral image pairs or NDVI time-series pairs are taken to calculate the change magnitudes independently using these shape parameters. In the change magnitude integration part, the units of different change magnitudes are assimilated by standard normalization. These unified change magnitudes are integrated with a weighting function. Finally, the changed and unchanged regions are acquired by comparing the change magnitude with a threshold.  Inter-annual trends in the NDVI time series derived from remote sensing data can be used to distinguish between natural land cover variability and land cover change [20]. Phenological land cover information remains in the shape features (e.g., amplitude and fluctuation) of the NDVI temporal curve at the inter-annual time scale. Throughout the year, seasonal phenologies affect these features during certain periods. Therefore, these features can be united to describe the annual characteristic of the NDVI temporal curve. The features are expressed as shape parameters, which are utilized as indicators to discriminate land cover change between two years. They are the phase angle cumulant, baseline cumulant, relative cumulation rate, and zero-crossing rate.
Phase angle cumulant (PAC). From the beginning of the year, the vegetation growth phases progress from vegetative growth to reproductive growth and withering [30]. The angle from one growth phase to another in the NDVI temporal curve was defined as the phase angle. It can be calculated with the NDVI changes and the time interval. As illustrated in Figure 2a, θ 1 is the first phase angle and θ 2 is the second phase angle of that NDVI temporal curve. Given n + 1 growth phases, n phase angles can be calculated from Equation (1), where V (i,e) and V (i,b) denote the ith NDVI values of the end phase and beginning phase, t (i,e) and t (i,b) denote the end and beginning time.
Baseline cumulant (BC). From the beginning of the growing season to its end, this period of NDVI time series acts as a measurement for most of the photosynthetic biomass accumulation. The straight line linking the beginning to the end of the growing season in the NDVI temporal curve was defined as the baseline. The baseline cumulant is defined by Equation (2), where ∆V i denotes the ith difference between the NDVI value and the baseline NDVI value, and B denotes the set of NDVI values above the baseline. The baseline is illustrated in Figure 2b.
Relative cumulation rate (RCR). Relative to the beginning of the growing season, the rates of the NDVI increase in the following time points represent a distinct absorption of near-infrared and infrared radiation for land cover types. The relative cumulation rate at a single time point is defined in Equation (3), where V 0 denotes the value at the beginning of the growing season and V i denotes the ith NDVI in the temporal curve composed of n NDVI values. The relative cumulation rate of the year is defined in Equation (4) as all rates of NDVI increase during one year: Zero-crossing rate (ZCR). The zero-crossing rate is an indicator that reflects the fluctuations of a curve in a given time interval, and the properties of the curves can be estimated based on the short-time average zero-crossing rate [31]. This indicator plays an important role in automatic recognition systems which process signals that pass through the level of zero, especially in surface electromyogram spectral analysis, voiced classification, and other fields [31][32][33]. Because of natural or artificial disturbance, some land cover statuses change rapidly during a time interval. Therefore, the NDVI values show turbulent amplitude, resulting in a high zero-crossing count around a baseline. Figure 2d illustrates the zero-crossings around the baselines in two periods. The zero-crossing rate can be calculated with Equations (5) and (6), where T is the total time interval, V i is the ith value of the NDVI temporal curve, and V m denotes the mean value of the NDVI baseline during that period: Remote Sens. 2018, 10, x FOR PEER REVIEW 5 of 20 zero-crossing rate can be calculated with Equations (5) and (6), where T is the total time interval, Vi is the ith value of the NDVI temporal curve, and Vm denotes the mean value of the NDVI baseline during that period:

Spectral and Temporal Difference Calculation
For the bi-temporal multispectral images, Pearson's correlation coefficient is adopted as a correlation operator to measure the similarity of two spectra [34]. The spectral change magnitude is calculated by Equation (7), where V is the reflectance of a pixel in band i, n is the band number, m denotes the mean value, and a and b denote the two spectra: The differences of the shape parameters are calculated for the NDVI temporal curves from two years (the time series has an interval of 16 days in each year). The temporal change magnitudes are derived from Equation (8), where Δf is the difference of the same shape parameter in two years, p is the order, and n is the number of vectors in that shape parameter. After calculation, four change magnitude images of the NDVI shape parameters are generated.

Spectral and Temporal Difference Calculation
For the bi-temporal multispectral images, Pearson's correlation coefficient is adopted as a correlation operator to measure the similarity of two spectra [34]. The spectral change magnitude is calculated by Equation (7), where V is the reflectance of a pixel in band i, n is the band number, m denotes the mean value, and a and b denote the two spectra: The differences of the shape parameters are calculated for the NDVI temporal curves from two years (the time series has an interval of 16 days in each year). The temporal change magnitudes are derived from Equation (8), where ∆f is the difference of the same shape parameter in two years, p is the order, and n is the number of vectors in that shape parameter. After calculation, four change magnitude images of the NDVI shape parameters are generated.

Magnitude Normalization
The differences between bi-temporal multispectral images are measured with the spectral correlation operator and the spectral change magnitude is calculated by subtracting the correlation coefficient. Therefore, the unit of the spectral change magnitude is 1. For temporal change magnitudes derived from four independent shape parameters, the units are degree, value, or even ratio. The values of these parameters range widely and differ greatly each other. In consideration of integration, the units of these spectral and temporal change magnitudes should be unified into one form. Standard normalization is done to assimilate these change magnitudes. As shown in Equation (9), where M is the change magnitude, M min /M max is the minimum/maximum of the change magnitude among the whole gray image. The spectral and temporal change magnitudes are transformed into relative change magnitudes with a range between 0 and 1. Remarkably, this linear transformation formula converts the range into 0-1, which remains as the relative change level among the different land cover change types.

Integrating Spectral-Temporal Change Magnitude
The unified change magnitudes are integrated with a weighting function (as shown in Equation (10)), where w is a weighting factor, I is the unified change magnitude, and n is the number of all change magnitudes. In this study, the spectral correlation change magnitude and four temporal change magnitudes are adapted to calculate one spectral-temporal change magnitude with this weighting function.

Change Region Discrimination by Automated Thresholding
After the spectral-temporal change magnitude integration, a single change magnitude grey image is created. Many applications require a binary image which can be used to find change regions. To address this issue, a threshold selection method known as "the histogram concavity analysis" is adapted to provide a good candidate threshold for detecting the changed regions and unchanged regions in the grey image. By analysing the histogram's concavity structure, this threshold selection method defines the optimal thresholds at the shoulders of the grey image histogram peaks [35]. Specifically, this method is an automated way to select a possible threshold. A more precise threshold may still need to be adjusted on the basis of this possible threshold. Then, comparing the grey image with the threshold, all regions that are greater than the threshold are detected as change regions.

Study Area, Data Preprocessing, and Algorithm Implementation
To test the performance of the proposed method, a case study was carried out in the northeast region (117 • 1 -120 • 34 E, 36 • 8 -38 • 1 N) of Shandong Province, China, as shown in Figure 3a. The average temperature of the study region is about 13 • C. This heat condition can satisfy the needs of the crop rotations. In addition, the average annual precipitation is about 555 mm, and the seasonal distribution of precipitation is uneven. Due to influences from biology, climate, and terrains, the soil is diversified in the region. Therefore, this study region has a variety of land cover types (main types are cropland, built-up, and water) as shown in Figure 3b, and the growths of different vegetation covers are varied. These natural factors and human disturbances such as wheat-corn rotation, urban expansion, and forest degradation cause uncertainty in the spectrum and NDVI temporal curves. All of these factors demand robust change detection methods.
Two Landsat images were downloaded (https://earthexplorer.usgs.gov) and transformed into the apparent reflectance (Figure 3c,d, the TM image consists of bands 1, 2, 3, 4, 5, and 7, the OLI image consists of bands 2, 3, 4, 5, 6, and 7). For the years 2010 and 2015, the NDVI time series drawn from the MODIS 16-day composite vegetation indices product (MOD13Q1) were obtained, totaling twenty-three NDVI images for each year (https://lpdaac.usgs.gov/data_access/data_pool). All of these NDVI images were projected from the sinusoidal projection to the WGS84 datum with the help of the MODIS reprojection tool (MRT). Then, the 30-m spatial resolution NDVI time series of 23 images in each year were reconstructed by the unmixing of these coarse NDVI time series [36][37][38][39]. The correlation coefficients of the generated fine NDVI time series and coarse NDVI time series were greater than 0.9 for each year. To evaluate the experimental results, the ground truth points (as shown in Figure 3c) were manually defined using the region of interest (ROI) tool from ENVI with the reference of high-resolution Google Earth images. Among these points, 2100 points had land cover conversions and 2523 points had no land cover conversions from the year 2010 to the year 2015. In the following, the change vector analysis is abbreviated as CVA, the image ratioing is abbreviated as IR; the spectral gradient-based change detection is abbreviated as SGD; the distance of Canberra is abbreviated as Dcanb, the NDVI gradient difference is abbreviated as NDVI-GD and the multiple shape parameter-based method is abbreviated as MSP.The methods of CVA, SGD, and IR used two Landsat images as inputs, NDVI-GD and Dcanb used two yearly fine NDVI time series as inputs. The MSP method used both two Landsat images and two yearly fine NDVI time series as inputs.
The MSP method was implemented with the IDL language and ran on a Windows 7 platform (Intel i3-3110M, 2 GB DDR3). The workflow was as follows (see Figure 1): (1) The spectral change magnitude image of the two Landsat images was calculated using Equation (6).  (6) of the first baseline was set to the means of the 1st to 13th NDVI values in the time series. The V m of the second baseline was set to the means of the 13th to 23rd NDVI values in the time series. The T in Equation (5) was set to 23. Based on these settings, each shape parameter of the NDVI temporal curves was calculated. Four temporal change magnitude images were calculated in Equation (8), with p being 2 for RCR and 1 for other parameters. (3) All of the change magnitude images were integrated into one single change magnitude image, making every w equal to 1 in Equation (10). (4) The histogram of the change magnitude image was analyzed to select a threshold at the shoulders of the histogram peaks, and the change regions were discriminated by comparing them with this threshold.

Change Magnitudes in the Ground Truth Points
In general, the change magnitudes in the real changed regions should have higher responses than in unchanged regions. The difference of the change magnitudes between the changed and unchanged regions could reflect the contrast quantitatively. Based on the two ground truth points group (one consisting of 2100 changed ground truth points, another consisting of 2523 unchanged ground truth points), the mean and median of the change magnitudes in each group and the standard deviation in all 4623 points were calculated. The differences between the changed and unchanged groups are listed in Table 1 with the units of percentage, and the standard deviation given separately. The difference from MSP was greater than the other compared methods. For MSP, the change magnitudes in the changed points were higher than in the unchanged points. This indicates that the MSP could generate high-contrast change magnitudes, which makes the discrimination of changes and no-changes easier.

Change Magnitudes in the Ground Truth Points
In general, the change magnitudes in the real changed regions should have higher responses than in unchanged regions. The difference of the change magnitudes between the changed and unchanged regions could reflect the contrast quantitatively. Based on the two ground truth points group (one consisting of 2100 changed ground truth points, another consisting of 2523 unchanged ground truth points), the mean and median of the change magnitudes in each group and the standard deviation in all 4623 points were calculated. The differences between the changed and unchanged groups are listed in Table 1 with the units of percentage, and the standard deviation given separately. The difference from MSP was greater than the other compared methods. For MSP, the change Remote Sens. 2018, 10, 1251 9 of 20 magnitudes in the changed points were higher than in the unchanged points. This indicates that the MSP could generate high-contrast change magnitudes, which makes the discrimination of changes and no-changes easier. Table 1. The means and medians of the change magnitudes in 2100 changed (C) and 2523 unchanged (UC) ground truth points and the standard deviations (STDs) in these 4623 ground truth points. All change magnitudes have been normalized and multiplied by the factor F = 10,000; CVA: change vector analysis; IR: image ratioing; SGD: spectral gradient-based change detection; Dcanb: distance of Canberra; NDVI-GD: NDVI gradient difference; MSP: multiple shape parameter-based method.   (3) from water to vegetation; (4) from muddy water to bare land; and (5) from crop to water. All these change magnitudes were normalized linearly to facilitate visual comparison and brightness is used to indicate change possibility.
In the Landsat image, the crop type shows a strong difference from the growing season (red) to the harvest season (white). Results of CVA and SGD showed high change magnitudes, while results from IR, NDVI-GD, and Dcanb showed low change magnitudes. The results of the MSP showed moderate but relatively low change magnitudes compared to the others. Because of urban expansion, some parts of the crop areas in 2010 were converted to built-up in 2015, while other parts remained as crop type, but after the harvesting. MSP results showed strong change magnitudes in the changed areas (crop to built-up) and light change magnitudes in the no-change areas (crop in growth to crop in harvest). This gave a great contrast between the change and no-change scenarios in the change magnitude image. The results of IR, CVA, and SGD were indistinguishable, while NDVI-GD and Dcanb were able to distinguish between changes and no-changes. However, they were not distinct enough. When water changed to vegetation, the MSP results showed a clear boundary between the changes and no-changes as shown in the Row 3 of Figure 4. When the yellow river shifted its path over five years, a part of the riverbed on one side was bare and a part of the other side was occupied by crops. The results of MSP gave a strong response to this change (from muddy water to bare land) and a light response to the no-change regions (from crop in growth to crop in harvest). The results of other methods were either indistinguishable or lacked the texture detail necessary to distinguish them accurately. In the case of crop to water, the advantage of keeping the texture detail of the real changes and abandoning no-changes was obvious for the MSP method.

Change Regions
The change magnitude images calculated from all test change detection methods were divided into changes and no-changes expressed with binary images, where the change is symbolized as white and the no-change is symbolized as black. All thresholds were automatically derived from the histogram concavity analysis. The ground truth image of the change regions was extracted interactively using the ROI tool from ENVI and Google Earth images. Partial change detection results of the study area are shown in Figure 5. It shows that the change result of MSP was closer to the ground truth image than other methods. Bi-temporal image-based change detection methods such as IR, CVA, and SGD could not provide accurate change results because of the phenological difference in multispectral Landsat images. Based on the 30-m NDVI time series, the pseudo-changes in the former were reduced (Figure 5g,h). However, some parts of the real changes were missing in the final change results of Dcanb and NDVI-GD. Due to the change magnitudes, these two methods had a low-contrast difference ( Figure 4) and it is difficult for an automated threshold selection method to accurately discriminate the real changes from no-changes. In contrast, the MSP results were more accurate in these conditions.

Change Regions
The change magnitude images calculated from all test change detection methods were divided into changes and no-changes expressed with binary images, where the change is symbolized as white and the no-change is symbolized as black. All thresholds were automatically derived from the histogram concavity analysis. The ground truth image of the change regions was extracted interactively using the ROI tool from ENVI and Google Earth images. Partial change detection results of the study area are shown in Figure 5. It shows that the change result of MSP was closer to the ground truth image than other methods. Bi-temporal image-based change detection methods such as IR, CVA, and SGD could not provide accurate change results because of the phenological difference in multispectral Landsat images. Based on the 30-m NDVI time series, the pseudo-changes in the former were reduced (Figure 5g,h). However, some parts of the real changes were missing in the final change results of Dcanb and NDVI-GD. Due to the change magnitudes, these two methods had a low-contrast difference ( Figure 4) and it is difficult for an automated threshold selection method to accurately discriminate the real changes from no-changes. In contrast, the MSP results were more accurate in these conditions.

Accuracy Assessment
The overall accuracy of the change detection result of MSP was 88.427% and the kappa coefficient was 0.764, as shown in Table 2. A comparison of the overall accuracy and kappa coefficient was also made for the other five change detection methods associated with the automated threshold, respectively, as shown in Table 3. The overall accuracy and kappa coefficient of MSP tended to be higher than those of the other five methods in this study area, indicating that the MSP change result derived from the automated threshold was more reliable and appropriate than the other methods. This advantage is especially important in automated change detection projects.

Accuracy Assessment
The overall accuracy of the change detection result of MSP was 88.427% and the kappa coefficient was 0.764, as shown in Table 2. A comparison of the overall accuracy and kappa coefficient was also made for the other five change detection methods associated with the automated threshold, respectively, as shown in Table 3. The overall accuracy and kappa coefficient of MSP tended to be higher than those of the other five methods in this study area, indicating that the MSP change result derived from the automated threshold was more reliable and appropriate than the other methods. This advantage is especially important in automated change detection projects. In addition, accuracy assessment curves were made to compare the performance of all six methods at different threshold levels. Based on the accuracy assessment curves, the relationship between the threshold levels and accuracies can be presented graphically [40]. In Figure 6, a series of thresholds which were converted into the deviation times of the standard deviation from the mean of the change magnitudes with a step of 0.25 are arranged on the x-axis. From the chart of the overall accuracy in Figure 6a and the chart of the kappa coefficient in Figure 6b, MSP achieved the highest levels in a wide range. With a threshold at seven deviations, MSP reached a peak (overall accuracy was about 88% and the kappa coefficient was about 0.76). Then, the accuracies declined gradually with the increase of the threshold levels. The other five methods showed a poor performance over the whole scene in the study area. CVA, SGD, and NDVI-GD reached a peak at four deviations for both the overall accuracy and the kappa coefficient. IR and Dcanb increased gradually along with the threshold levels and reached a best overall accuracy of 0.65. The kappa coefficients of these two methods were low throughout the entire range. This comparison indicates that MSP was effective in a wide range of thresholds and could simultaneously achieve the best overall accuracy and kappa coefficient.   In addition, accuracy assessment curves were made to compare the performance of all six methods at different threshold levels. Based on the accuracy assessment curves, the relationship between the threshold levels and accuracies can be presented graphically [40]. In Figure 6, a series of thresholds which were converted into the deviation times of the standard deviation from the mean of the change magnitudes with a step of 0.25 are arranged on the x-axis. From the chart of the overall accuracy in Figure 6a and the chart of the kappa coefficient in Figure 6b, MSP achieved the highest levels in a wide range. With a threshold at seven deviations, MSP reached a peak (overall accuracy was about 88% and the kappa coefficient was about 0.76). Then, the accuracies declined gradually with the increase of the threshold levels. The other five methods showed a poor performance over the whole scene in the study area. CVA, SGD, and NDVI-GD reached a peak at four deviations for both the overall accuracy and the kappa coefficient. IR and Dcanb increased gradually along with the threshold levels and reached a best overall accuracy of 0.65. The kappa coefficients of these two methods were low throughout the entire range. This comparison indicates that MSP was effective in a wide range of thresholds and could simultaneously achieve the best overall accuracy and kappa coefficient.

The Influence of Different Order p on Change Magnitudes and Accuracy
The change magnitudes of the NDVI time series are calculated based on PAC, BC, RCA, and ZCR in MSP. The order p in Equation (8) determines the differences of all these temporal change magnitudes (assuming p 1 , p 2 , p 3 , and p 4 were of the order p when calculating the differences of PAC, BC, RCA, and ZCR, respectively). Therefore, the contrast between the changed and unchanged pixels (abbreviated as CCU), the difference between the changed and unchanged pixels (abbreviated as DCU) in the change magnitude of MSP, and the overall accuracy will be affected by the various combinations of the order p. In this section, the values 1, 2, and 3 were assigned to p 1 , p 2 , p 3 , and p 4 , and a total of 81 (3 4 ) combinations were selected to test their influence on the contrast, the difference between changed and unchanged pixels, and the overall accuracy (every weighting factor w was set to 1). As shown in Figure 7, the DCU of combinations ((1, 1, 1, 1), (1, 1, 2, 1), and (1, 1, 3, 1)) were high and the CCU of combinations ((3, 2, 3, 1), (2, 3, 3, 1), and (3, 3, 3, 1)) were also high. Although the higher DCU and CCU could show a high-contrast change magnitude, it is clear that the highest DCU and CCU did not bring the highest accuracy. From the scattered points in Figure 7, we can see that the combinations (2, 2, 3, 1), (3, 2, 3, 1), and (2, 3, 3, 1) are superior for DCU, CCU, and accuracy.

The Influence of Different Order p on Change Magnitudes and Accuracy
The change magnitudes of the NDVI time series are calculated based on PAC, BC, RCA, and ZCR in MSP. The order p in Equation (8) determines the differences of all these temporal change magnitudes (assuming p1, p2, p3, and p4 were of the order p when calculating the differences of PAC, BC, RCA, and ZCR, respectively). Therefore, the contrast between the changed and unchanged pixels (abbreviated as CCU), the difference between the changed and unchanged pixels (abbreviated as DCU) in the change magnitude of MSP, and the overall accuracy will be affected by the various combinations of the order p. In this section, the values 1, 2, and 3 were assigned to p1, p2, p3, and p4, and a total of 81 (3 4 ) combinations were selected to test their influence on the contrast, the difference between changed and unchanged pixels, and the overall accuracy (every weighting factor w was set to 1). As shown in Figure 7, the DCU of combinations ((1, 1, 1, 1), (1, 1, 2, 1), and (1,  1, 3, 1)) were high and the CCU of combinations ((3, 2, 3, 1), (2, 3, 3, 1), and (3, 3, 3, 1)) were also high. Although the higher DCU and CCU could show a high-contrast change magnitude, it is clear that the highest DCU and CCU did not bring the highest accuracy. From the scattered points in Figure 7, we can see that the combinations (2, 2, 3, 1), (3,2,3,1), and (2, 3, 3, 1) are superior for DCU, CCU, and accuracy.

The Functions of the Weighting Factor w in the Integration
The change magnitudes of the proposed shape parameters describe the differences between the land cover conversions in several dimensions. The weighting factors of these change magnitudes determine the proportions of particular change magnitudes in the final MSP change magnitude, and thus, they affect the final change detection accuracy. Values from 1 to 5 were

The Functions of the Weighting Factor w in the Integration
The change magnitudes of the proposed shape parameters describe the differences between the land cover conversions in several dimensions. The weighting factors of these change magnitudes determine the proportions of particular change magnitudes in the final MSP change magnitude, and thus, they affect the final change detection accuracy. Values from 1 to 5 were assigned to each weighting factor w, and the overall accuracies of MSP change detection under 3125 (5 5 ) combinations of all five factors were calculated (every p was set to 1). Figure 8 gives the value distributions of the w factor in the different accuracy groups, where w0, w1, w2, w3, and w4 are the weighting factors of the change magnitudes of SC (spectral correlation), PAC, BC, RCR, and ZCR, respectively. The figure shows that w0 gradually moved from 1 to 5 as accuracy increased. In the groups in which the accuracy was greater than 0.7, most w0 values were 5. The value distributions of w1 and w2 were similar when the accuracy was greater than 0.7. The values 1, 2, and 3 are advisable for w1 and w2 to achieve a higher accuracy. Based on the value distributions of w3, the value shifted sharply from 5 to 1 as the increasing of accuracy. The preferred values of w3 were 1 and 2. When accuracy was greater than 0.8, it is clear that most of w4 gathered around the weighting factor value of 5, and the proportions of w4 in values 4 and 5 were close to each other. Therefore, weighting factor values of 4 and 5 are recommended for w4. It seems that the combinations of (w0, w1, w2, w3, w4) set as (5, 1/2/3, 1/2/3, 1/2, 4/5) could achieve a higher overall accuracy. The accuracy of all these combinations was greater than 0.78, and 49% of them were greater than 0.8.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 20 assigned to each weighting factor w, and the overall accuracies of MSP change detection under 3125 (5 5 ) combinations of all five factors were calculated (every p was set to 1). Figure 8 gives the value distributions of the w factor in the different accuracy groups, where w0, w1, w2, w3, and w4 are the weighting factors of the change magnitudes of SC (spectral correlation), PAC, BC, RCR, and ZCR, respectively. The figure shows that w0 gradually moved from 1 to 5 as accuracy increased. In the groups in which the accuracy was greater than 0.7, most w0 values were 5. The value distributions of w1 and w2 were similar when the accuracy was greater than 0.7. The values 1, 2, and 3 are advisable for w1 and w2 to achieve a higher accuracy. Based on the value distributions of w3, the value shifted sharply from 5 to 1 as the increasing of accuracy. The preferred values of w3 were 1 and 2. When accuracy was greater than 0.8, it is clear that most of w4 gathered around the weighting factor value of 5, and the proportions of w4 in values 4 and 5 were close to each other. Therefore, weighting factor values of 4 and 5 are recommended for w4. It seems that the combinations of (w0, w1, w2, w3, w4) set as (5, 1/2/3, 1/2/3, 1/2, 4/5) could achieve a higher overall accuracy. The accuracy of all these combinations was greater than 0.78, and 49% of them were greater than 0.8.

Change Detection with a Single Shape Parameter
The proposed shape parameters can be used independently to perform change detection. Figure 9 shows the accuracy assessment curves of these indicators with MSP as a comparison. Based on the spectral characteristic of the Landsat image, spectral correlation performed better than NDVI temporal curve-based shape parameters. The baseline cumulant achieved an overall accuracy of 0.75, which was the highest level using the NDVI time series separately. However, the overall

Change Detection with a Single Shape Parameter
The proposed shape parameters can be used independently to perform change detection. Figure 9 shows the accuracy assessment curves of these indicators with MSP as a comparison. Based on the spectral characteristic of the Landsat image, spectral correlation performed better than NDVI temporal curve-based shape parameters. The baseline cumulant achieved an overall accuracy of 0.75, which was the highest level using the NDVI time series separately. However, the overall accuracy and kappa coefficient of these separate indicators were lower than MSP in a wide range, from 3 to 13. Although the accuracy of each shape parameter was low for most thresholds, the overall accuracy of MSP was greater than 0.8. This indicates that the integration of these separate shape parameters improves the accuracy of the final change detection results. Each of them contributes to the accuracy of integration to varying degrees.
Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 20 accuracy and kappa coefficient of these separate indicators were lower than MSP in a wide range, from 3 to 13. Although the accuracy of each shape parameter was low for most thresholds, the overall accuracy of MSP was greater than 0.8. This indicates that the integration of these separate shape parameters improves the accuracy of the final change detection results. Each of them contributes to the accuracy of integration to varying degrees.

The Contributions of Shape Parameter to Land Cover Conversions
As mentioned in Section 3.1, land cover types such as crop, built-up, and water were dominant in the study area. Therefore, the functions of each shape parameter in the conversions between these land cover types should be major concerns. The spectrum and 30-m NDVI time series of four land cover types (crop, built-up, muddy water, and water) were sampled with a pixel point in years 2010 and 2015, as shown in Figure 10. The change magnitudes of these land cover conversions were calculated using each shape parameter and normalized linearly to range between 0 and 1. The land cover types crop, built-up, muddy water, and water are abbreviated as C, B, M, and W, respectively, in Figure 11. The conversions from 2010 to 2015 are expressed as "type-type" in the following. In the case of conversions from crop to other types, as shown in Figure 11a, the baseline cumulation contributed a high change magnitude to changes and low change magnitude to

The Contributions of Shape Parameter to Land Cover Conversions
As mentioned in Section 3.1, land cover types such as crop, built-up, and water were dominant in the study area. Therefore, the functions of each shape parameter in the conversions between these land cover types should be major concerns. The spectrum and 30-m NDVI time series of four land cover types (crop, built-up, muddy water, and water) were sampled with a pixel point in years 2010 and 2015, as shown in Figure 10. The change magnitudes of these land cover conversions were calculated using each shape parameter and normalized linearly to range between 0 and 1. accuracy and kappa coefficient of these separate indicators were lower than MSP in a wide range, from 3 to 13. Although the accuracy of each shape parameter was low for most thresholds, the overall accuracy of MSP was greater than 0.8. This indicates that the integration of these separate shape parameters improves the accuracy of the final change detection results. Each of them contributes to the accuracy of integration to varying degrees.

The Contributions of Shape Parameter to Land Cover Conversions
As mentioned in Section 3.1, land cover types such as crop, built-up, and water were dominant in the study area. Therefore, the functions of each shape parameter in the conversions between these land cover types should be major concerns. The spectrum and 30-m NDVI time series of four land cover types (crop, built-up, muddy water, and water) were sampled with a pixel point in years 2010 and 2015, as shown in Figure 10. The change magnitudes of these land cover conversions were calculated using each shape parameter and normalized linearly to range between 0 and 1. The land cover types crop, built-up, muddy water, and water are abbreviated as C, B, M, and W, respectively, in Figure 11. The conversions from 2010 to 2015 are expressed as "type-type" in the following. In the case of conversions from crop to other types, as shown in Figure 11a, the baseline cumulation contributed a high change magnitude to changes and low change magnitude to The land cover types crop, built-up, muddy water, and water are abbreviated as C, B, M, and W, respectively, in Figure 11. The conversions from 2010 to 2015 are expressed as "type-type" in the following. In the case of conversions from crop to other types, as shown in Figure 11a, the baseline cumulation contributed a high change magnitude to changes and low change magnitude to no-changes. The following is the relative cumulation rate, which performed effectively for no-changes. The phase angle gave a high contrast between changes and no-changes. The spectral correlation gave a low contrast between C-C and C-B. The zero-crossing rate performed well in conversions with water or muddy water. As shown in Figure 11b, in the case of built-up to others, baseline cumulation, phase angle, and relative cumulation rate worked together to give a high change magnitude in B-C. The zero-crossing rate and spectral correlation enhanced the change magnitudes of B-M and B-W, which was weak in B-M. In the case of changes from muddy water to other types, as shown in Figure 11c, the spectral correlation performed well, but not strongly enough to highlight the changes among M-M and M-W. The phase angle gave the lowest change magnitude. The rest were complementary to improve the magnitude of changes. In the case of changes from water to other types, as shown in Figure 11d, the zero-crossing rate gave a high contrast between changes and no-changes. Baseline cumulation performed poorly between W-M and W-W, as did the phase angle between W-B and W-M. Similarly, the relative cumulation rate did not work in discriminating W-M and W-W. Despite these shortcomings, the final change magnitude was enhanced by the integration of different change magnitudes.
Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 20 no-changes. The following is the relative cumulation rate, which performed effectively for no-changes. The phase angle gave a high contrast between changes and no-changes. The spectral correlation gave a low contrast between C-C and C-B. The zero-crossing rate performed well in conversions with water or muddy water. As shown in Figure 11b, in the case of built-up to others, baseline cumulation, phase angle, and relative cumulation rate worked together to give a high change magnitude in B-C. The zero-crossing rate and spectral correlation enhanced the change magnitudes of B-M and B-W, which was weak in B-M. In the case of changes from muddy water to other types, as shown in Figure 11c, the spectral correlation performed well, but not strongly enough to highlight the changes among M-M and M-W. The phase angle gave the lowest change magnitude. The rest were complementary to improve the magnitude of changes. In the case of changes from water to other types, as shown in Figure 11d, the zero-crossing rate gave a high contrast between changes and no-changes. Baseline cumulation performed poorly between W-M and W-W, as did the phase angle between W-B and W-M. Similarly, the relative cumulation rate did not work in discriminating W-M and W-W. Despite these shortcomings, the final change magnitude was enhanced by the integration of different change magnitudes.

The Applicability of Data Quality for Landsat Images and NDVI Time Series
Previous research has demonstrated that multi-spectral images should be pre-processed by radiometric and atmospheric correction for quantitative remote sensing applications [41][42][43]. For the purpose of land cover change detection, some researchers have argued that this type of pre-processing is unnecessary in cases where the images are acquired under the same imaging conditions. Noise reduction techniques for NDVI time series have been widely used to improve the

The Applicability of Data Quality for Landsat Images and NDVI Time Series
Previous research has demonstrated that multi-spectral images should be pre-processed by radiometric and atmospheric correction for quantitative remote sensing applications [41][42][43]. For the purpose of land cover change detection, some researchers have argued that this type of pre-processing is unnecessary in cases where the images are acquired under the same imaging conditions. Noise reduction techniques for NDVI time series have been widely used to improve the monitoring ability of vegetation dynamics. Therefore, the pre-processing changes the data quality of Landsat images or NDVI time series by reducing the difference of radiation or phenology. The applicability of the proposed method for various inputs is important for its wide range of application.
Two kinds of Landsat images (under radiometric correction or both radiometric and atmospheric correction) and two kinds of NDVI time series (with the Savitzky-Golay filter applied or without it) were taken as inputs for MSP. Figure 12 shows that the accuracies of UU and US (UU denotes data with no atmosphere correction and no Savitzky-Golay filter; US denotes data with no atmosphere correction and with the Savitzky-Golay filter) were lower than CS and CU (CS denotes data with atmosphere correction and the Savitzky-Golay filter; CU denotes the data with atmosphere correction and without the Savitzky-Golay filter). This indicates that atmosphere correction is vital for MSP. When atmosphere correction was performed, the accuracy reached a higher level, as seen by the blue and pink curves in Figure 12. When the thresholds were lower than nine times the deviation, the mean of the change magnitudes was selected to detect the changes and the Savitzky-Golay filtering was not necessary. By increasing the thresholds, the Savitzky-Golay filtering becomes important to obtain a high accuracy. This result indicates that using Landsat images after atmosphere correction is essential for MSP and the filtering of the NDVI time series is not indispensable. Landsat images or NDVI time series by reducing the difference of radiation or phenology. The applicability of the proposed method for various inputs is important for its wide range of application. Two kinds of Landsat images (under radiometric correction or both radiometric and atmospheric correction) and two kinds of NDVI time series (with the Savitzky-Golay filter applied or without it) were taken as inputs for MSP. Figure 12 shows that the accuracies of UU and US (UU denotes data with no atmosphere correction and no Savitzky-Golay filter; US denotes data with no atmosphere correction and with the Savitzky-Golay filter) were lower than CS and CU (CS denotes data with atmosphere correction and the Savitzky-Golay filter; CU denotes the data with atmosphere correction and without the Savitzky-Golay filter). This indicates that atmosphere correction is vital for MSP. When atmosphere correction was performed, the accuracy reached a higher level, as seen by the blue and pink curves in Figure 12. When the thresholds were lower than nine times the deviation, the mean of the change magnitudes was selected to detect the changes and the Savitzky-Golay filtering was not necessary. By increasing the thresholds, the Savitzky-Golay filtering becomes important to obtain a high accuracy. This result indicates that using Landsat images after atmosphere correction is essential for MSP and the filtering of the NDVI time series is not indispensable.

Conclusions
In order to calculate a high-contrast change magnitude, a novel change detection method based on multiple shape parameters is proposed in this study. As demonstrated in the experiments, the shape parameters of the NDVI temporal curve and spectrum could be used to detect the inter-annual land cover conversions. These shape parameters were from reconstructions of phenological parameters, existing indicators, and mathematical indexes. Some of them were shown to be effective in particular transformations of land cover types. For example, the zero-crossing rate was valid for detecting the fluctuation of the NDVI temporal curve of water. The phase angle cumulant showed an obvious effect in transformations having crop types. However, none of them was more powerful than the others obviously in detecting various land cover conversions. After the integration of these shape parameters, a high-contrast change magnitude image could be calculated. MSP's change magnitude of the NDVI temporal curve is derived from the annual change

Conclusions
In order to calculate a high-contrast change magnitude, a novel change detection method based on multiple shape parameters is proposed in this study. As demonstrated in the experiments, the shape parameters of the NDVI temporal curve and spectrum could be used to detect the inter-annual land cover conversions. These shape parameters were from reconstructions of phenological parameters, existing indicators, and mathematical indexes. Some of them were shown to be effective in particular transformations of land cover types. For example, the zero-crossing rate was valid for detecting the fluctuation of the NDVI temporal curve of water. The phase angle cumulant showed an obvious effect in transformations having crop types. However, none of them was more powerful than the others obviously in detecting various land cover conversions. After the integration of these shape parameters, a high-contrast change magnitude image could be calculated. MSP's change magnitude of the NDVI temporal curve is derived from the annual change characteristics during the whole year rather than the accumulation of the NDVI differences in each single time point of the year. Accordingly, the offset between the different conflicting differences of single time points is avoided and the changes show a higher response compared with the no-change scenario. Therefore, pseudo-changes are reduced further and the real changes and non-changes are more easily discriminated.
The result of the case study showed that an overall accuracy of about 88% was achieved by MSP, and its kappa coefficient was 0.76. The change regions in the change magnitudes of dominant land covers were distinguishable, which could provide a broad range to the set thresholds. This property is important to detect changes accurately among the various change types. The change magnitude images were distinct from the main land cover conversions in the study area, which ensured that the change regions and no-change regions had a high-contrast. The pseudo-changes caused by the phenological differences of the Landsat images could be eliminated further. This means that we might be able to conclude that MSP is feasible for land cover change detection in regions suffering from pseudo-change errors. The MSP maintained high accuracies in a wide range of threshold levels and depends less on filtering. This might demonstrate that MSP is robust for either automated or optional change detection projects.
The MSP method also has disadvantages. Firstly, this method needs both bi-temporal Landsat images and 30-m NDVI time series as inputs, which is a restriction for some land cover change detection programs that lack these datasets. Another limitation is the low efficiency compared to traditional direct comparison-based change detection methods. Several change magnitudes derived from multiple shape parameters should be calculated separately and integrated into one in MSP. Third, MSP only calculates a binary change image, cannot provide a "from-to" image, and lacks the ability to detect changes within one year. In the future, fewer but more effective features of NDVI temporal curves might be a possible way to improve its efficiency. Due to the many parameters that are involved in MSP, more research needs to be done to identify their detailed function. It would also be valuable to perform the MSP method based on image objects.