OIC-MCE: A Practical Land Cover Mapping Approach for Limited Samples Based on Multiple Classifier Ensemble and Iterative Classification

Land cover samples are usually the foundation for supervised classification. Unfortunately, for land cover mapping in large areas, only limited samples can be used due to the time-consuming and labor-intensive sample collection. A novel and practical Object-oriented Iterative Classification method based on Multiple Classifiers Ensemble (OIC-MCE) was proposed in this paper. It systematically integrated object-oriented segmentation, Multiple Classifier Ensemble (MCE), and Iterative Classification (IC). In this method, the initial training samples were updated self-adaptively during the iterative processes. Based on these updated training samples, the inconsistent regions (ICR) in the classification results of the MCE method were reclassified to reduce their uncertainty. Three typical case studies in the China-Pakistan Economic Corridor (CPEC) indicate that the overall accuracy of the OIC-MCE method is significantly higher than that of the single classifier. After five iterations, the overall accuracy of the OIC-MCE approach increased by 5.58%–8.38% compared to the accuracy of the traditional MCE method. The spatial distribution of newly added training samples generated by the OIC-MCE approach was relatively uniform. It was confirmed by ten repeated experiments that the OIC-MCE approach has good stability. More importantly, even if the initial sample size reduced by 65%, the quality of the final classification result based on the proposed OIC-MCE approach would not be greatly affected. Therefore, the proposed OIC-MCE approach provides a new solution for land cover mapping with limited samples. Certainly, it is also well suited for land cover mapping with abundant samples.


Introduction
Spatially explicit land cover data are important for characterizing anthropogenic activity and biogeographical and eco-climatic diversity [1,2]. Remote sensing is the main approach to produce such land cover data, because it can obtain information on the Earth's surface at a relatively low cost, in a short time and with high efficiency [3,4]. Following the continuous development of satellite sensor technology, computer technology, artificial intelligence technology, and other related technologies, In 2013, Chinese President Xi Jinping put forward the "Belt and Road" (B&R) initiative [30]. The Economic Corridor (EC) has become the basis and entry point for the implementation of this initiative [31]. The planning and construction of ECs could not be separated from the scientific cognitions for the spatial patterns and change processes of the ecological environment along these corridors [32,33]. The long time-series land cover datasets just provided this information from a macro perspective. Except for the FROM-GLC [34] and GlobeLand30 [35], the spatial resolution of existing global-scale land cover products is relatively coarse (≥ 300 m), which made it difficult to describe the spatial pattern of ecosystems in complex areas, such as mountainous area, geographic transition zone [36,37]. For the FROM-GLC and GlobeLand30, their time span was relatively short (≤ 10 years), which could not show the change processes of ecosystem over a long time. Therefore, it is necessary to develop novel and practical land cover mapping and change detection methods to support the generation of land cover datasets. Currently, Landsat images (30 m) can be available over large areas and with a 16-day revisiting period as well as over a forty-year time span. For generating time-series land cover datasets in EC under the B&R initiative, Landsat image is the optimal data source. Most of the ECs span many countries and regions. Limited samples would be the typical characteristic of land cover mapping in these regions. Therefore, it is even more urgent to propose new land cover mapping methods with limited samples.
In conclusion, this paper attempts to propose a novel and practical land cover mapping method with limited samples based on the MCE and IC method. It would realize the self-adaptive updating of samples and gradually improve the classification accuracy during the iterative process. The rest of this paper was organized as follows. Section 2 introduces the study area and the used data. Section 3 describes the core steps of the proposed method in detail. Section 4 shows the newly added samples, classification results, and their accuracy. Section 5 further discusses the effects of the main factors on the conduction of the proposed method. Section 6 draws the conclusion of this paper.

Study Area
The China-Pakistan Economic Corridor (CPEC) is one of the six economic corridors under the B&R initiative and is also the flagship project [38]. Three Landsat scenes (Path/Row: 150/35, 151/40, and 152/43) in the World Reference System-2 (WRS-2), located in the north, middle, and south section of the CPEC, were selected as the typical experiment regions to verify the performances of the proposed algorithm in this study, as shown in Figure 1. The spatial distribution patterns of land cover in three experiment regions are different. Region 1 (150/35, Figure 1b) is located at the intersection area of the Himalayas, the Karakoram, and the Hindu Kush mountains. The mountain landscape is its most prominent feature. Region 2 (151/40, Figure 1c) is located in the transitional area between Indus Plain and Balochistan Plateau. The belt area from the northeast to the southwest of this region is Indus Plain. Balochistan Plateau lies to the northwest, while the Thar Desert occupies the southeast. The most typical plain agricultural landscape and desert landscape in Pakistan were included in Region 2. The mouth of Indus River and the Karachi (the largest city of Pakistan) included in Region 3 (152/43, Figure 1d), which was selected as a representation of the urban landscape and estuarine landscape in Pakistan. In general, three experiment regions covered almost all landscape types in CPEC, which would help to test the performances of the proposed algorithm in various landscapes.

Data
For each Landsat scene, four cloud-free Landsat-8 OLI satellite images were used in this study. Multi-temporal satellite images could describe the difference in phenological characteristics of different vegetation types at different times, which would contribute to improving classification accuracy [39]. The surface reflectance data of all selected Landsat images (Level-2 Science Products) were downloaded from the United States Geological Survey (USGS) EarthExplorer platform. They were calibrated radiometrically, atmospherically, and geometrically. Therefore, pre-processing was not conducted.
In addition to the satellite images, topographical data, including elevation, slope, and aspect, were used as auxiliary data in land cover mapping. The 30 m-resolution SRTM DEM data were adopted. Slope and aspect were calculated from the SRTM DEM data using ERDAS software.

Methodology
A novel and practical Object-oriented Iterative Classification method based on the Multiple Classifiers Ensemble (OIC-MCE) was proposed in this paper. This method systematically integrated object-oriented segmentation, MCE, and IC. Under the condition of limited samples, the training samples were updated self-adaptively during the iterative processes. Based on these updated samples, the ICRs in the classification result of the MCE method were reclassified to reduce their uncertainty. Finally, the accuracy of land cover mapping was improved gradually by the IC method. The flowchart of the proposed OIC-MCE method is shown in Figure 2. The MCE, IC and self-adaptive updating of training samples are the core steps of this method. The following sections introduce this in detail. detailed description of this classification system is shown in Table S1. For each experiment region, all 168 these 23 land cover classes are less likely to exist simultaneously. For example, permanent ice/snow 169 cannot exist in Region 2 and Region 3. Therefore, there may be significant differences in the land 170 cover classes in each experiment region.
collect the initial training samples and validation samples [41]. The salient advantage of using Google

173
Earth is that high-resolution historical images are available for free, while also helping to identify 174 samples better [42]. Usually, the spatial resolution of panchromatic images in Google Earth are 175 typically below 1m, and some of them are less than 0.5m. Considering the following two criteria: (1) 176 homogeneous spatial distribution, (2) the representative of samples, each selected experiment region 177 was divided into approximately 400 grids with 10 km × 10 km, and then 1-2 samples with the 178 prominent land cover classes were collected in each grid [43]. For grids with the heterogeneous landscape, at least two samples were collected, while for grids with the homogeneous landscape,

Initial Training Samples Collection
The classification system of China land cover for carbon budget was adopted in this study [40]. Considering the local characteristics of land cover in CPEC and the interpretation abilities of satellite images, it was further modified. The final classification system contained 23 land cover classes. The detailed description of this classification system is shown in Table S1. For each experiment region, all these 23 land cover classes are less likely to exist simultaneously. For example, permanent ice/snow cannot exist in Region 2 and Region 3. Therefore, there may be significant differences in the land cover classes in each experiment region.
High spatial-resolution historical images in 2015, from the Google Earth platform, were used to collect the initial training samples and validation samples [41]. The salient advantage of using Google Earth is that high-resolution historical images are available for free, while also helping to identify samples better [42]. Usually, the spatial resolution of panchromatic images in Google Earth are typically below 1m, and some of them are less than 0.5m. Considering the following two criteria: (1) homogeneous spatial distribution, (2) the representative of samples, each selected experiment region was divided into approximately 400 grids with 10 km × 10 km, and then 1-2 samples with the prominent land cover classes were collected in each grid [43]. For grids with the heterogeneous landscape, at least two samples were collected, while for grids with the homogeneous landscape, only one sample was needed. For some land cover classes, such as impervious surface, open water, flood plain, which cannot be dominant in the area, samples were randomly collected in their main distribution regions. The geographical location of the collected samples was located at the geometric center of the concentrated distribution regions of these prominent land cover classes in each grid, preferably not mixed with other land cover classes [41]. More importantly, their geographical location must be in areas where the land cover class did not change in 2015. For each experiment region, approximately 600 training samples were collected. Although the sample size was not small and did not match the condition of limited samples, it could provide a reference for further experiments on limited samples.
Based on the overlaying analysis of training samples and segmented objects (see Section 3.2 for details), sample objects were obtained. The area of the sample object should be between 10 and 20 pixels. Because if its area was too small, it might not represent the prominent land cover class in the grid. If its area was too large, other land cover classes might be mixed into it. Each selected sample object, it would be replaced with a neighboring object of the same class that meets this condition, when its area did not meet this condition.

Object-Oriented Multi-Resolution Segmentation and Classification Features Selection
In order to avoid the "Salt and Pepper Noise" often existing in pixel-based classification results, and to take fully into consideration spectral, shape, and texture features, the object-oriented classification method was adopted [44][45][46]. The multi-resolution segmentation algorithm, adopted by the eCognition Developer 8.7, was adopted in this study. All selected OLI images were incorporated in the segmentation process, and all bands of selected OLI images weighted equally. "Scale", "Shape", and "Compactness" parameter are three core parameters for this algorithm, which would determine directly the final segmentation result [47]. By referring to the similar research of author's team [48], the parameters in this study were set as 25, 0.1, and 0.7, respectively.
In this study, spectral features, spectral indices, topographic features, texture features, and shape features were chosen for land cover mapping, which was shown in Table 1. The spectral indices included NDVI (Normalized Difference Vegetation Index) [49], NBR (Normalized Burn Ratio) [50], NDBI (Normalized Difference Built-up Index) [51] and NDWI (Normalized Difference Water Index) [52], which were sensitive to vegetation, fire disturbance, building, and water, respectively. All of the selected features were calculated in the eCognition 8.7 platform. For spectral features, spectral indices, and topographic features, the mean value of the segmented object was used for classification.

Features Calculation Method References
Spectral Features Surface reflectance of b1, b2, b3, b4, b5, b6, and b7 The mean value of each segmented object Spectral Indices Notes: b1, b2, b3, b4, b5, b6, and b7 denote, respectively, the coastal, blue, green, red, near-infrared and two mid-infrared bands of Landsat-8 OLI images. The GLCM is gray level co-occurrence matrix, which has a good ability to distinguish texture features of natural objects. The b v is the border length of the image object, and the p v is the area of the image object.

Multiple Classifiers Ensemble (MCE)
For the MCE method, five common classifiers were adopted, including Random Forest (RF), Support Vector Machine (SVM), Decision Tree (DT), Naïve Bayes (NB) and K-Nearest Neighbor (KNN). Each classifier required to set several optimal initial parameters to maximize the advantages of each classifier, except NB. In this study, the trial and error method was used to determine the optimal values for these initial parameters. Meanwhile, the parameter values of these classifiers in similar works also provided some references for this study [53,54]. The optimal values of the initial parameter for each classifier are shown in Table 2. After setting the initial parameters of each classifier, training was conducted with the initial sample objects and selected features, and then five classification results were generated by using these trained classifiers. For the integration of classification results of the five classifiers, the voting method was adopted [23]. For each object, the land cover class with the most votes was identified as the integrated land cover class.

Self-Adaptive Updating of Training Samples
In the classification result of MCE, there are two notably different regions: CR and ICR. As mentioned above, the classification results of CR have relatively high reliability, which would be the potential regions for extracting new training samples. Generally, the classification results of ICR needed to further reduce uncertainty. However, there was one exception, that is, the classification results of the majority of classifiers (such as four classifiers in this study) are completely consistent, and different from the classification results of other classifiers (such as one classifier in this study). We call it as "consistent regions in inconsistent regions (CRinICR)" (Figure 3). It can be considered that the classification result of the majority of classifiers in CRinICR was credible. While for other classifiers, the CRinICR was just the weak region for their previous classification processes. If training samples can be extracted from CRinICR, this would help reduce the uncertainty of ICR. Based on the above analysis, we proposed a self-adaptive updating method for training samples.
The proposed OIC-MCE method updated independently the training samples for each classifier by adding new training samples selected from the CRinICR. A random stratification sampling approach was adopted to collect training samples in CRinICR. In order to make the training samples statistically representative, the sample size of each land cover class needs to be consistent with their proportion in the initial samples, and the size of newly added training samples for a single iteration should also be the same as the size of initial samples. The area of each newly added sample object also needed to meet the condition between 10 and 20 pixels. Inevitably, a few error samples existed in the initial samples and the newly added training samples. To prevent the influences of these error samples from propagating to the subsequent iterations, they should be eliminated. Due to the high reliability of CR, only the training samples located in the CR would be used in the later iterations in this study.

257
The iterative training sample, iterative region, and iterative time are key factors to be considered 258 in the IC. There are two schemes for iterative training samples. One is that only newly added samples 259 were used for training during each iteration. The other is that newly added samples and previously 260 used samples were adopted for training in each iteration. Considering that newly added samples 261 were only collected from CRinICR where a certain classifier was difficult to classify correctly in the 262 previous iterative processes, these samples might be difficult to fully reflect the overall characteristics of land cover. Therefore, the latter scheme was adopted. Because the classification result of CR has 264 relatively high credibility, only the ICR was selected as the iterative region. This strategy would 265 contribute to improving the computational efficiency of the proposed method. For the iterative time,

266
if it was set too small, the advantage of IC could not be fully realized. Conversely, if it was set too 267 big, the iterative process would spend more time, and more importantly, the last few iterations would 268 not contribute much to improving the classification accuracy. Finally, the iterative time was set to 10 269 in this study. Although this setting might not be the optimal value, it would help to discover the 270 optimal iterative time for land cover mapping in large areas.

Iterative Classification (IC)
The iterative training sample, iterative region, and iterative time are key factors to be considered in the IC. There are two schemes for iterative training samples. One is that only newly added samples were used for training during each iteration. The other is that newly added samples and previously used samples were adopted for training in each iteration. Considering that newly added samples were only collected from CRinICR where a certain classifier was difficult to classify correctly in the previous iterative processes, these samples might be difficult to fully reflect the overall characteristics of land cover. Therefore, the latter scheme was adopted. Because the classification result of CR has relatively high credibility, only the ICR was selected as the iterative region. This strategy would contribute to improving the computational efficiency of the proposed method. For the iterative time, if it was set too small, the advantage of IC could not be fully realized. Conversely, if it was set too big, the iterative process would spend more time, and more importantly, the last few iterations would not contribute much to improving the classification accuracy. Finally, the iterative time was set to 10 in this study. Although this setting might not be the optimal value, it would help to discover the optimal iterative time for land cover mapping in large areas.

Integration of Iterative Classification Result
For each iteration, the classification result was divided into CR and ICR. All CRs formed during the iterative process would be merged to generate the final classification result in this proposed method. In general, only when the classification results of all classifiers were completely consistent, could they be identified as CR. In fact, if the classification results of the majority of classifiers (such as CRinICR) were completely consistent, their credibility is very high. More importantly, CRinICR was the cradle of new training samples. Therefore, when at least four classification results of the five classifiers were consistent, it could be considered as CR in this study. Inevitably, for the last iteration, some ICRs still existed in the classification results. These ICRs would also be merged into the final classification result and their land cover class would be determined by the voting method described in Section 3.3.

Accuracy Assessment
An error matrix was adopted in this study to verify the effectiveness of the proposed OIC-MCE method. The validation samples and training samples were collected independently. For each experiment region, approximately 330 validation sample objects were collected from high spatial-resolution satellite images by team members and international students familiar with the local situation. The spatial distribution of validation sample objects was shown in Figure 2. The error matrix was obtained by checking validation sample objects one by one, and the overall accuracy and total disagreement (sum of quantity disagreement (QD) and allocation disagreement (AD)) [55] of three experiment regions were calculated synchronously. For the accuracy of each land cover class, precision and recall were adopted in this study [56].

The Self-Adaptively Updated Training Samples
The spatial distribution of the self-adaptively updated training samples for the SVM classifier is shown in Figure 4. Overall, their spatial distribution was relatively uniform for each experiment region. Statistically, the number of newly added training samples showed a downward trend with the increase in iterative time (Figure 4a-c). Generally speaking, the area of the ICR was constantly declined following with the iterative process, because more and more new CRs were extracted from previous ICRs. Correspondingly, the CRinICR used to collect new training samples was also decreased. Inevitably, fewer new training samples would be added following the iterative process. Additionally, the size of the newly added samples for some land cover classes was smaller than the size of the initial samples, or even there was no new sample in some iterations. The sample size of each land cover class was relatively close to their actual area ratio, except for some land cover classes with a small area, or which were difficult to map, such as impervious surface, grassland. Figure 5 shows the CR generated in each iteration for three experiment regions. By increasing the iterative time, the area of CR gradually decreased. The total area ratio of CR in three experiment regions after the first five iterations reached 95.24%, 94.31%, and 97.43%, respectively ( Figure 5). This also shows that the increase in iterative time after five iterations was not helpful to improve the overall accuracy [29] because the areal ratio of the remaining ICR was only about 5%. The homogeneous regions were accurately identified in the initial classification and the first few iterations, while the heterogeneous regions and geographical transition zones were gradually identified with the iterative process ( Figure 6).

The Accuracy of the OIC-MCE and Each Single Classifier in the Iterative Process
The overall accuracy and total disagreement of the OIC-MCE and each classifier for each iteration are shown in Figures 7 and 8, respectively. Overall, following the iterative process, the overall accuracy of all classifiers and OIC-MCE increased (Figure 7), and their total disagreement decreased ( Figure 8). The overall accuracies of RF, SVM, and KNN were better than NB, and the accuracy of DT was the worst in three experiment regions (Tables 3-5). With the increase of iterative time, the growth rate of the overall accuracy of each classifier gradually decreased, which also appeared in total disagreement. After about five iterations, the accuracy of each classifier was stable. For the OIC-MCE method, the overall accuracy of three experiment regions after five iterations increased by 8.38%, 7.99%, and 5.58% compared to the accuracy of the initial MCE method, respectively. Therefore, it can be concluded that five iterations were the optimal iterative time for the proposed OIC-MCE method, which is consistent with Han et al.'s findings [29]. The overall accuracy of the OIC-MCE method is significantly higher than that of the single classifier in this study (Tables 3-5). The accuracy of the OIC-MCE was 0.96~3.65% higher than that of the single classifier. The result also indicate that OIC-MCE could not only give full play to the advantage of each classifier but also reduced the potential influences of its shortcomings to some extent.  identified with the iterative process ( Figure 6).

331
For the OIC-MCE method, the overall accuracy of three experiment regions after five iterations 332 increased by 8.38%, 7.99%, and 5.58% compared to the accuracy of the initial MCE method, proposed OIC-MCE method, which is consistent with Han et al.'s findings [29]. The overall accuracy 335 of the OIC-MCE method is significantly higher than that of the single classifier in this study (Tables 336 [3][4][5]. The accuracy of the OIC-MCE was 0.96~3.65% higher than that of the single classifier. The result 337 also indicate that OIC-MCE could not only give full play to the advantage of each classifier but also 338 reduced the potential influences of its shortcomings to some extent. 339 340 Figure 7. The overall accuracy of each classifier and OIC-MCE algorithm in the iterative process for

The Accuracy of Each Land Cover Class in Iterative Classification
The precision and recall of the OIC-MCE method for each land cover class in experiment region 2 are shown in Figure 9. Except for dry land (Figure 9h) and impervious surface (Figure 9j), the precision and recall of remaining land cover classes showed an upward trend with the increase of iterative time. It is further indicated that the proposed OIC-MCE method was helpful for improving the classification accuracy of most land cover classes. Similarly, when the iteration time reached 5, neither the precision nor the recall changed much. With the increase in iterative time, the precision of dry farmlands decreased, but its recall increased, which indicates that although more dry farmlands were identified by the iterative processes, some non-dry farmlands were also included (Figure 9h). The impervious surface showed an opposite trend, indicating that the iterative processes improved continually the purity of impervious surfaces, but more and more impervious surfaces were misclassified into other land cover classes (Figure 9j). The precision and recall of each land cover class were compared between the OIC-MCE method and the RF classifier ( Figure 9). They have similar change trends in the precision and recall for two classification methods. The OIC-MCE method obtained slightly higher precision than the RF classifier for flood plain, dry farmland, sparse vegetation, and bare soil. The recall of the OIC-MCE method is slightly higher than that of the RF classifier in evergreen shrubland, flood plain, and sparse vegetation.  The convenient collection and efficient utilization of training samples have always been one of 377 the core work for supervised classification [9]. A series of methods and strategies have been proposed 378 for the acquisition and utilization of samples. Geo-Wiki [13] and Global Mapper [34] were their 379 typical representatives. They either obtained training samples in the form of crowdsourcing [13] or 380 provided tools to serve the convenient collection of training samples [34]. However, for land cover 381 mapping in large areas, even if some advanced methods for sample collection were available, only 382 limited training samples could be collected and used due to the time-consuming, labor-intensive 383 difference in the classification system and the lack of quality control for these method [10]. How to 384 make full use of these limited samples to produce higher-quality land cover products needed to be 385 further explored [15]. The proposed OIC-MCE method is just one of the solutions for the above

395
The processing time of the OIC-MCE method for a Landsat scene was 1.5 times the total processing 396 Figure 9. The comparison of the precision and recall of each land cover class between the proposed OIC-MCE method and the RF classifier in Region 2. (a-m) represent the precision and recall for evergreen broadleaf forest, evergreen shrubland, grassland, herbaceous wetland, open water, flood plain, paddy field, dry farmland, orchard, impervious surface, sparse vegetation, salina, and bare soil, respectively.

Adaptively Updating of Training Samples
The convenient collection and efficient utilization of training samples have always been one of the core work for supervised classification [9]. A series of methods and strategies have been proposed for the acquisition and utilization of samples. Geo-Wiki [13] and Global Mapper [34] were their typical representatives. They either obtained training samples in the form of crowdsourcing [13] or provided tools to serve the convenient collection of training samples [34]. However, for land cover mapping in large areas, even if some advanced methods for sample collection were available, only limited training samples could be collected and used due to the time-consuming, labor-intensive difference in the classification system and the lack of quality control for these method [10]. How to make full use of these limited samples to produce higher-quality land cover products needed to be further explored [15]. The proposed OIC-MCE method is just one of the solutions for the above problems. It can self-adaptively discover the appropriate samples for each classifier, entirely relying on the information provided by existed training samples and multiple classifiers, and without any manual operations. More importantly, this method not only saved a large amount of time and labor for sample collection but also got an acceptable classification accuracy by verifying the classification result.

The OIC-MCE
The proposed OIC-MCE method integrates the MCE method and the IC method. On the one hand, it solved the problem that the conventional MCE method cannot further reduce the uncertainty of the ICR. On the other hand, it also found a new way of adding new information to the IC method. The processing time of the OIC-MCE method for a Landsat scene was 1.5 times the total processing time of five stand-alone classifiers, because about half of segmented objects were classified at least twice. If parallel computing technology was considered in the proposed OIC-MCE method, the processing time would be shorter in the future. The initial sample size is a significant parameter for the performance of the proposed method. Therefore, further analysis and identification of its optimal value were required.

Effects of Initial Sample Size on the Performance of the OIC-MCE Method
We always expect to obtain the best mapping results by using the minimum sample size, because sample collection is time-consuming and labor-intensive [15]. Whether the OIC-MCE method can satisfy the demands of land cover mapping with limited samples needed to be further discussed. In other words, if the initial sample size was reduced, would this affect the final classification accuracy of the OIC-MCE method? Therefore, the initial sample size of three experiment regions was gradually reduced by 5% in this study until the sample size was only 5% of the initial sample size. Moreover, all the experiments were repeated ten times. The experiment result is shown in Figure 10. The final classification accuracy was decreased, with the reduction of the training sample size. When the training sample size was reduced to 35% of the initial sample size, the classification accuracy was declined by no more than 5%, which is still 0.2%-3.36% higher than the classification accuracy of the MCE method alone for three experiment regions. This implies that when the training sample size declines by about 65% (approximately 200 training samples), there is no significant impact on the classification accuracy by using the OIC-MCE method. Therefore, this provides a new solution for land cover mapping with limited samples.

402
We always expect to obtain the best mapping results by using the minimum sample size, because 403 sample collection is time-consuming and labor-intensive [15].

The Stability of the OIC-MCE Method
Random stratification sampling was adopted for updating training samples in the proposed OIC-MCE method. There is a potential risk regarding the stability of this method because of the difference in the location of the new samples generated in each iteration by this method. Therefore, the stability of this method was verified by repeated experiments in three experiment regions. After ten repeated experiments, it can be seen that the OIC-MCE method has good stability ( Figure 11). After five iterations, the classification accuracy of this method in three experiment regions was 82.13% ± 0.72%, 82.08% ± 0.50%, and 84.32% ± 0.52%. The range of accuracy fluctuations of a single iteration was located in (0.25%, 0.72%). The fluctuation of overall accuracy does not exceed 1%.
Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 21 ten repeated experiments, it can be seen that the OIC-MCE method has good stability ( Figure 11).

425
After five iterations, the classification accuracy of this method in three experiment regions was 82.13%

Merits and Limitations of the Proposed Method
Land cover mapping in large areas often faces a lack of training samples. This paper proposed a novel and practical land cover mapping method (OIC-MCE) under the condition of limited samples. It can adaptively update training samples by using the MCE method and the IC method. The accuracy of land cover classification was improved significantly compared with the single classifier and traditional MCE method, through classification experiments in the CPEC. More importantly, it significantly reduced the requirement of training samples for land cover mapping in large areas, especially for regions where training samples were scarce. Of course, the proposed OIC-MCE method is also well suited for regions with abundant training samples.
Although the proposed method provided a feasible solution for land cover mapping using Landsat images with limited samples, and also improved significantly the classification accuracy, the following four aspects still need to be considered in future studies: (1) Only five classifiers were considered in this paper. Other classifiers, or how to effectively integrate these classifiers [25], or how to optimize the classifier parameters, exceeded the scope of this paper. However, they could be considered in further research. (2) Following the continuous maturity of deep learning algorithms and their wider application in high-resolution land cover mapping [57,58], whether the introduction of deep learning algorithms into this method would contribute to improving the classification accuracy and reduce the limitation of samples still needs to be further studied. (3) It is inevitable that error samples existed in the initial training sample [28]. Although part of the initial error samples was eliminated in the proposed method, it is necessary to further analyze the influences of the existence of error samples for the conduction of the OIC-MCE method. The robustness of the OIC-MCE method will be verified by adding some noise samples or changing the land cover class of part of the initial samples in the future. (4) The proposed OIC-MCE method was validated as a very good contribution for land cover mapping using Landsat images with limited samples. In the future, it can be applied in other images with different spatial resolutions as well as the spectral and temporal resolutions, including ones with a higher spatial resolution or lower spatial resolution.

Application of the Proposed OIC-MCE Method in the Land Cover Mapping of CPEC
The core object of the proposed OIC-MCE method was to serve land cover mapping in large areas, especially for Economic Corridors of B&R initiative. However, this algorithm was proposed for land cover mapping in a single mapping unit (such as a tile of Landsat-8 OLI image). The acquisition time of satellite images used in adjacent mapping units is often different, which would result in the spatial discontinuity of their classification results [48]. Therefore, when the OIC-MCE method was applied to land cover mapping in large areas, it must consider the problem of spatial discontinuity. Another special classification method has been proposed by our research team to deal with this problem. The core idea of that method was similar to OIC-MCE. The CR and ICR were firstly extracted from the overlapping parts of the adjacent mapping units. New training samples were generated from CR by random stratification sampling. The iterative classification method was adopted to reduce ICR. Currently, the paper about this method is being prepared. The preliminary land cover map for CPEC in 2015 has been completed by using the OIC-MCE method ( Figure 12). After preliminary verification, the overall accuracy of land cover map of CPEC is more than 80%. Further accuracy assessment is still in progress.  result in the spatial discontinuity of their classification results [48]. Therefore, when the OIC-MCE 463 method was applied to land cover mapping in large areas, it must consider the problem of spatial 464 discontinuity. Another special classification method has been proposed by our research team to deal 465 with this problem. The core idea of that method was similar to OIC-MCE. The CR and ICR were 466 firstly extracted from the overlapping parts of the adjacent mapping units. New training samples 467 were generated from CR by random stratification sampling. The iterative classification method was 468 adopted to reduce ICR. Currently, the paper about this method is being prepared. The preliminary 469 land cover map for CPEC in 2015 has been completed by using the OIC-MCE method ( Figure 12).

470
After preliminary verification, the overall accuracy of land cover map of CPEC is more than 80%.

471
Further accuracy assessment is still in progress.

476
For land cover mapping in large areas, remote sensing data were no longer a constraint because

Conclusions
For land cover mapping in large areas, remote sensing data were no longer a constraint because a large number of different types of medium or coarse resolution satellite images are freely available. Samples have gradually become a bottleneck in land cover mapping in large areas, due to the time-consuming and labor-intensive for sample collection. How to improve the accuracy of land cover mapping based on available limited samples without increasing the input of manpower, time and money, has become a common concern in this field. This paper introduced a novel and practical land cover mapping method (OIC-MCE) for limited samples. It can self-adaptively update the training samples, entirely relying on the existed training samples and multiple classifiers, and without any manual operations. Meanwhile, it also solved the problem that the conventional MCE method cannot further reduce the uncertainty of the ICR. At least as important, it finds a new way of adding new information to the IC method. Spatially, the distribution of the newly added training samples was relatively uniform. Statistically, the number of newly added training samples showed a downward trend with the increase in iterative time. The overall accuracy of the OIC-MCE method is significantly higher than that of the single classifier. After five iterations, the overall accuracy of the OIC-MCE method increased by 5.58%-8.38% compared to the accuracy of the initial MCE method in three experiment regions. Ten repeated experiments in three experiment regions have shown that the OIC-MCE method has good stability. There is no significant impact on the classification accuracy by using the OIC-MCE method, when the training sample size declined by about 65%. Certainly, the proposed OIC-MCE method not only provided a new solution for land cover mapping with limited samples but is also well suited for land cover mapping with abundant samples. Three case studies and the application in land cover mapping in CPEC proved the effectiveness and stability of the OIC-MCE method for land cover mapping using Landsat images with limited samples. Whether it can be applied to global land cover mapping, further validation studies in different climatic, topographical, and ecological regions are needed, due to the diversity of global ecosystems and the complexity of climate types.