A highly efficient temporal-spatial probability synthesized model from multi-temporal remote sensing for paddy rice identification

ABSTRACT This article develops a temporal-spatial probability synthesized model (TSPSM), in which a metric describing the characteristic of temporal and spatial information is defined to map paddy rice distribution. The purpose is to reduce the effect of cloud contamination on classification. The error matrix and Kappa were used as accuracy measurement. Results showed that TSPSM obtained higher accuracy with significant difference from error matrices of the other two conventional methods, post comparison classification with post-classification comparison and majority voting. Moreover, smaller window was suitable for the area with higher fragmentation, while the larger was suitable for the area with lower fragmentation. It was concluded that TSPSM could help to improve the potentials of temporal optical image to map crops.


Introduction
Paddy rice is one of the most important staple crops, occupying approximately 15% area of the world's arable land (Xiao et al., 2005). It is grown on flooded soils as a main source to release carbon gas, which contributes to global climate change. Moreover, more foods are expected to meet the demand of continuously growing population (Foley et al., 2011). It is essential for crop management and policy decisionmaking to obtain acreage of paddy rice accurately and timely (Sun et al., 2014).
Remote sensing technology has been widely applied for mapping paddy rice because of its abundant temporal and spatial information to record the paddy growing condition. Time series of remote sensing images are usually preferred to monitor paddy rice by tracking its unique temporal variation (Xiao et al., 2002). Generally, there are three types of researches on multi-temporal data for mapping vegetation. The first type is to integrate all date images for classifying (Bailey & Frohn, 2001). Studies attempted to optimize input data sets by using the prior images transformations algorithm (Panigrahy, Ray, & Panigrahy, 2009), for example, principal component analysis (PCA). Bailey and Frohn (2001) using another transformation method, the minimum noise fraction transformation, mapped the wet vegetation and estimated the area of wild rice. An experiment conducted by Panigrahy et al. (2009) achieved good performance with overall accuracy more than 86.0% by using different data sets derived by PCA. Second, the multi-temporal change detection is used to make full use of temporal and spatial information for mapping. Chen, Chen, Shi and Yamaguchi (2012) combined change detection and classification methods to map land cover, reaching to an overall accuracy of 87.01%.
The last approach involving temporal vegetation indices (VIs) is employed to extract crop (Azar et al., 2016). These phenology-and pixel-based methods were used widely to map paddy rice by using VIs, such as normalized difference vegetation index and enhanced vegetation index (EVI) (Chen, Huang, & Hu, 2011;Wardlow, Egbert, & Kastens, 2007;Xiao et al., 2006Xiao et al., , 2002Xiao et al., , 2005. Temporal analysis of VIs is able to capture specific stage of crop phenology, which can present the different characteristics among land covers. VIs derived from various satellite sources, such as AVHRR, MODIS and Landsat, were efficient indices to reflect the vegetation, and all metrics were frequently used to investigate the crop phenology (Bischof et al., 1992;Chen et al., 2011;Roy et al., 2014;Wardlow et al., 2007). Sensors with low spatial resolution have an ability of providing the highly frequent resolution and broad coverage image, which are often employed to mapping paddy rice at regional or global scale. Xiao et al. (2002) introduced an approach for mapping paddy rice focusing on a subregion in Jiangsu province of China, using serial SPOT satellite vegetation data. Xiao et al. (2005), Xiao et al. (2006)) further used multi-temporal MODIS images to analyze the temporal correspondence between land surface water index (LSWI) and VIs for mapping paddy rice in larger area. Both experiments presented that the linear coefficient of determination (R 2 ) varied from 0.42 to 0.87 at province scale. Sun, Huang, Huete, Peng and Zhang (2009) obtained paddy rice distribution with temporal information of EVI and LSWI based on paddy rice growth calendar. Relative errors of estimated paddy rice related to agricultural statistic varied from 3% to 21%. Unfortunately, this algorithm is unsuitable for monitoring interannual study due to the influence of the cloud contamination. Peng, Huete, Huang, Wang and Sun (2011) improved the efficiency of this methodology by gap-filling cloud contaminated pixels with an interpolation filter to solve the issue of all pixels contaminated by the cloud, which was not considered in the Xiao et al. (2005), Xiao et al. (2006) and Sun et al. (2009). By assessing with the land use data, the algorithm acquired R 2 from 0.34 to 0.58 for single, early and late paddy rice. Bridhikitti and Overcamp (2012) further optimized the algorithm conducted by Xiao et al. (2005), Xiao et al. (2006) that is based on temporal profiles of vegetation and water content. However, these algorithms exploiting the coarse resolution data are feasible to underestimate or overestimate paddy rice of the small-scale field (Jeong et al., 2012).
To improve the low performance of crop identification from coarse resolution images in fragmented field, the finer resolution images, such as Landsat 5/7/ 8, have been applying. Crop identification based on single phase and multi-temporal consists of the main measurements. For previous case, various algorithms, such as neural network models and regression tree models, are applied to the single image using different classifiers or classification models to produce a map of paddy rice (Diuk-Wasser et al., 2007;Fang, 1998;Kuenzer & Knauer, 2013). However, the contemporaneous growing vegetation, like corn and soybean, contributes the confusion to the paddy rice that it is hard to ensure the accuracy of paddy rice identification. Hence, multi-temporal satellite sensor images are applied in mapping paddy rice. Similar to the literatures of mapping paddy rice with coarse resolution images, various VIs were also integrated to extract paddy rice from temporal pixel-and phenology-based approach (Dong et al., 2015(Dong et al., , 2016Van Niel, McVicar, Fang, & Liang, 2003). Dong et al. (2015) used time series Landsat data and phenology-based algorithm to develop a Landsat-based paddy rice mapping system achieving satisfied performance from 84% to 95% of overall accuracy and 0.6-0.9 of Kappa, respectively. Dong et al. (2016) extended their algorithm (Dong et al., 2015) to map paddy rice at a larger scale. This phenology-and pixel-based algorithm for mapping paddy rice obtained a good performance with overall accuracy more than 95% for all countries. Besides these algorithms, Yan and Roy (2014) recently proposed an object-based method named as variable region-based geometric active contour method to map crop. This automated computational methodology employed multi-temporal Landsat data to extract the crop field in California, which reached to overall accuracy of 90.1% for field class.
In general, most of the existing approaches used finer resolution images to exploit "cloud-free" images to identify paddy rice. However, applicability of those approaches has limit for large-scale operation, because that it is difficult to get sufficient images for crop identification (Cheng, Shen, Zhang, Yuan, & Zeng, 2014;Jin et al., 2013;Ju & Roy, 2008;Zhu, Chen, Gao, Chen, & Masek, 2010). To avoid the effect of cloud contamination, synthetic aperture radar images were used for mapping rice planting area (Oguro, Suga, Takeuchi, Ogawa, & Konishi, 2001;Shao et al., 2001;Zhang, Wang, Wu, Qi, & Salas, 2009), which could solve the issue of the cloud contamination in optical images. However, its application has been limited by data availabilities of few satellites in-orbit operation (Shao et al., 2001;Zhang et al., 2009). Therefore, it is still a big challenge to eliminate the influences of cloud contamination to map crops at large scale.
The purpose of this paper is to develop a simple and feasible model named temporal-spatial probability synthesized model (TSPSM) to map paddy rice, which aims to reduce the effect of cloud contamination on classification. In this model, the temporalspatial probability (TSP) is defined from posterior probability, and the thematic map of TSP is produced. Then, paddy rice is mapped by extracting from the thematic of TSP with an accurate and binary-based threshold. The error matrix and kappa coefficient were used to assess the accuracy of the proposed model.

Temporal-spatial probability definition
The posterior probability quantifies the membership to one specific land cover, which is the basis for construction of TSPSM. This probability is a procedural variable generated during classifier operation, such as support vector machine (SVM) and maximum likelihood classifier (Rannala & Yang, 1996). Figure 1 shows the process of TSPSM for identification of paddy rice. According to the BBCH phenological scale (Meier, 1997), the principal growth of paddy rice mainly experiences 10 stages: germination, leaf development, tillering, stem elongation, boosting, heading, flowering and anthesis, development of fruit, ripening and senescence. At these stages, the posterior probability maps will be generated by SVM classifier.
Here, we mark the posterior probability, P t ij , of each pixel which belonging to paddy rice as the "spatial probability" of one pixel. P ij t , is calculated by averaging P t ij of the surrounding pixels within a predefined window (m × m, e.g. 3 × 3, see Equation (1)). The P ij t aims to reduce the spectral variability for paddy rice identification. Thus, "spatial probability" maps (i.e. P ij 1 , P ij 2 ,. . ., P ij t , 0 ≤ t ≤ 9) describe the likelihood of pixels belonging to paddy rice integrating spectral and spatial information at the principal growth season. The TSP, P ij , is defined as the average of "spatial probability" of all dates to describe the likelihood of a pixel belonging to paddy rice (see Equation (2)).
where P ij t denotes the posterior probability of paddy rice at i th column and j th row in the image at t th growing stage, m is the size of predefined window and P ij t is the average posterior of P ij t belonging to paddy rice within a predefined window (m × m).
where n is the number of remote sensing images used for identification of paddy rice. Therefore, P ij is the metric considering spatial and temporal information to demonstrate quantitatively the pixel which belongs to paddy rice. It is an inevitable issue that pixels are contaminated by clouds and shadows in remote sensing images. These pixels will be marked as contamination pixels. Then these will be defined as "NULL" and ignored during the process of TSP calculation.

Study area and materials
The study area is located in the south of Liaoning province of China (40°44′-41°30′N and 121°36′-122°24′E). It occupies approximately 2882 km 2 ( Figure 2). The topography of this area is flat, and the agricultural landscape is relatively regular. Paddy rice is the dominant crop in the study area, which is cultivated one time a year. The growing season of paddy rice covers from May to October.
Six Landsat 8 Operational Land Imager (8 OLI) images at 30 m and Thermal Infrared Sensor (TIRS) images at 100 m resolution were obtained from USGS Earth Explorer data center (http://earthex plorer.usgs.gov/), which covered the critical paddy rice growing season. These images were used to produce paddy rice maps and related posterior probability maps. Moreover, the TIRS images were the necessary auxiliary data source for cloud detection and removal.  The detailed description for images was listed in Table 1. BBCH phenological scale (Meier, 1997) is used to describe the phenological growth stages of paddy rice, which can precisely define the phenological growth stages for plant in two-digit code (Ramírez, Fischer, Davenport, Pinzón, & Ulrichs, 2013). According to Meier (1997), the entire cycle of paddy rice subdivided into 10 stages. Based on the BBCH, the images at six key stages that the vegetation characters of paddy rice were easily distinguished, which support the sufficient information for identifying paddy rice.
The OLI and TIRS images were systemically and geometrically corrected to UTM-WGS84 following the professional direction (http://landsat. usgs.gov/documents/LDCM-DFCB-004.pdf). FLAASH tool in ENVI version 5.1 was used to carry out the radiometric correction and atmospheric correction for the images (Cheng et al., 2014).

EUROPEAN JOURNAL OF REMOTE SENSING
The Function of mask algorithm (Fmask) as one of the latest effectively tools is able to detect clouds and related shadows in images, which is based on spectral information related to cloud and shadow physical properties (Zhu, Wang, & Woodcock, 2015;Zhu & Woodcock, 2012). This algorithm has three steps to locate the cloud and related shadow: potential cloud layer-pass one, potential cloud layer-pass two, potential shadow layer. Pass one and two are used to detect potential cloud pixels, which are based on cloud physical properties. Then cloud shadows are detected through its darkening effects in near-infrared (NIR) bands by using flood-fill transformation algorithm. Finally, the object-based cloud and shadow are matched to obtain the distribution of cloud and shadow.
Sample is a basis to guarantee the reasonable and quantitative assessment on crop identification (Stehman & Wickham, 2011). In this paper, two sources of high-resolution images, unmanned aerial vehicle (UAV) and GaoFen-1 satellite (GF-1, meaning "high resolution" in Chinese) (Li, Chen, Tian, Huang, & Feng, 2015) were acquired as reference data to validate the performance of TSPSM. UAV images at 10 cm were acquired from 10-12 August. Nineteen spots from multispectral GF-1 images at 8 m resolution on 10 August were downloaded from China Centre for Resources Satellite Data and Application as complementary materials (http://www.cresda.com/n16/index.html).
These two types of images had a sufficient quality including spectral and spatial features to identify the paddy rice accurately. About 738 samples in the extent of reference data were investigated and recorded the types of crop. Assisted by these ground survey data collected synchronously with UVA mission, the paddy rice was visually interpreted and digitized in vector format. According to the rule of the maximum area within one pixel (Lobell & Asner, 2004), the digitized paddy rice was converted to raster format at 30 m resolution in accordance with the resolution of Landsat OLI images. Figure 3 describes the process of our experiment, which mainly follows four steps: data processing, SVM classification, and TSPSM classification and accuracy assessment.

OLI images classification
The SVM as an excellent supervised classification method can fit an optimal hyperplane among individual land covers only depending on a small training sample set (Foody & Mathur, 2004). The training samples were selected directly from the remote sensing image with the help of field survey data. The 2-7 bands of OLI images as inputs were used to acquire the land cover classification maps and posterior probability maps. In Table 2, the classification schema at individual date was described. Those, for paddy rice, were the periods of irrigation on 23 May and 8 June, when the paddy rice fields were covered by water. Thus, the special classification schema at these two stages included the category of water. In the schema, the terminology of "others" referred to the type of land cover that was easy to be distinguished and we did not concern about, such as, built-ups. The initial results for individual dates from SVM were assessed from error matrix using reference data. We utilized the classification results and associated posterior probability maps for the further operations when the overall accuracy of classification results attained above 80%.  Classification result of SVM Reference data Land cover 5/23 6/8 7/26 8/11 9/12 9/28 Paddy rice

Paddy rice identification
The TSP, P ij , was calculated following Equations (1) and (2). Initially, the window size was set as 3 × 3. Double-window flexible pace searching (DFPS) developed by Chen, Gong, He, Pu and Shi (2003) was introduced to calculate the threshold. It is a feasible and efficient algorithm for binarize the grayscale image to extract the target class. According to the DFPS, typical samples are selected to calculate the "test success rate" that is defined to assess the performance of each potential threshold for identifying the change and no change pixels. In our study, 9457 paddy rice samples and the surrounding 3576 nonpaddy rice samples were selected from TSP map to calculate the threshold for extracting paddy rice. When the "test success rate" is over 90.0%, the corresponding threshold would be accepted to distinguish the paddy rice and non-paddy rice. The final paddy rice map was binarized, 1 for paddy rice, 0 for the background information.

Performance analysis
Two conventional methodologies, post classification comparison (PCC) and majority voting (MV), were also used for identifying paddy rice (Chen et al., 2012;Louisa & Suen, 1997;Lu, Mausel, Brondízio, & Moran, 2004). All three methods were analyzed. PCC is based on the land cover transformation of classified results at two dates in the pixel to pixel manner, and target land covers are determined following the principal of changes from date1 to date2 (Lu et al., 2004). According to phenology, paddy rice showed spectral characteristics as water on 23 May and vigorous characteristics of vegetation on 12 September in the study area. We applied the two date images as input features into SVM classification to identify the paddy rice. The interpreted paddy manner is that the pixels identified as water and vegetation for the respective phenology stages, germination and development of fruit, were labeled as paddy rice.
MV is capable of using the multi-dates images, which is simple to realize pattern recognition as effective as complex combination strategies, such as Bayesian, logistic regression, fuzzy integral and neural network (Chen et al., 2012). This method takes the rule of a majority vote to define the target land cover (Louisa & Suen, 1997). According to the phenology knowledge, paddy rice showed spectral characteristics as water in 23 May, 8 June and characteristics of vegetation from 26 July to 28 November. If pixels were labeled as same as the paddy rice, they would be marked as potential paddy rice. Then if four out of six class labels of pixels were labeled as potential paddy rice, the class label of a pixel on the classification map would be considered as paddy rice depending on the majority principle.

Accuracy assessments
The error matrix is a quantitative measurement to assess the accuracy of classification derived from remote sensing images. In this measurement, four quantitative metrics, including overall accuracy, producer's and user's accuracy, and Z-test with Kappa coefficients (Congalton & Green, 2008), were calculated to assess the accuracy. Overall accuracy represents the proportion of the pixel mapped correctly, which is calculated by Equation (3). Kappa is also a metric of how well of the classification result agrees with actual data, which is calculated by Equation (4). Moreover, Kappa value also is introduced for Z calculation (see Equation (5)). The other two assessment metrics, user's accuracy and producer's accuracy, can explain the characteristics of commission and omission errors, respectively, for the specific land cover type.

Overall accuracy
where n is the number of rows in the matrix, x ii is the number of observation in row i and column i, x iþ and x þi are the marginal totals row i and column i, respectively, and N is the total number of the observation. Congalton and Green (2008) proposed a Z-test method to test significant difference between two different error matrices. This test is possible to compare two algorithms statistically (referred to as pairwise test) to explore which algorithm could produce higher accuracy. Moreover, if the Z ≥ 1.96 (at the 95.0% confidence interval level), it means that two error matrices (i.e. error matrix #1, error matrix #2) derived from two algorithms for mapping land cover are significantly different. It should draw a conclusion that the classification algorithm with higher accuracy is better than the other one (Congalton & Green, 2008;Congalton, Oderwald, & Mead, 1983). The Z-test metric is defined as follows: where K 1 , K 2 denote the estimates of Kappa statistic for error matrix #1 and error matrix #2, respectively. δ 2 1 , δ 2 2 are the corresponding estimates of the variance of K 1 , K 2 , respectively.
In this paper, accuracy of mapping paddy rice was evaluated by using error matrices measurement of TSPSM, PCC and MV.  May and 12 September, images RGB channels were SWIR1, NIR and RED; (c)-(h) were the classification results of TPSFM corresponding to 1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11 window size; (i) was the result of PCC; (j) was the result of MV.

Results and discussion
Accuracy assessment Figure 4(a) is the thematic map of TSP for paddy rice at 3 × 3 window size. The lower TSP in lighter gray color appears at the boundary of the field, gradually growing into the center of the field, demonstrating the changes of TSP trend. The final distribution of paddy rice map is extracted by the threshold of 0.58 for TSP of paddy rice at "test success rate" of 92.04%, derived from DFPS (Figure 4(b). Table 3 presents the quantitative accuracy assessment results. For 3 × 3 window size, according to the Z-test at 95% confidence interval, the result that Z of 100.51 (>1.96) between TSPSM and PCC, and Z of 55.11 (>1.96) between TSPSM and MV, respectively, demonstrates that TSPSM is significantly different from the conventional methods. Also, TSPSM gained Kappa of 0.81 and overall accuracy of 90.47%, which are higher than PCC and MV, respectively. These results demonstrate that TSPSM perform better than the other two methods.
As expected, TSPSM successfully reduces the effect of cloud contamination on mapping paddy rice. To show the advantage of the proposed model, an area contaminated seriously by cloud was chosen ( Figure 5), which was marked in a red color rectangle in Figure 4. As shown from the results in Figure 5(b), the image on 12 September is contaminated by cloud seriously. Figure 5(i) shows paddy rice is misclassified as non-paddy rice because only one date image can offer valid data to identify paddy rice. Figure 5 (d), (j) is the identification result of TSPSM and MV, which presents better performance against PCC. Figure 5(d) intuitively illustrates that the proposed model could effectively map paddy rice from eliminating the influence from cloud contamination because the cloud data can be made up by other date images.
For pixel-based classification, the "salt and pepper" is a common issue caused from the spectral heterogeneity (Oguro et al., 2001). In TSPSM, the TSP, P ij , possesses the spatial traits to quantify the relationship between the center pixel and surrounding pixels, which can reduce the spectral heterogeneity with a similar function of a spatial filter (Ji, Ma, Twibell, & Underhill, 2006;Switzer, 1980;Zhu et al., 2010). Even though the spatial filter can be applied for PCC and MV as a postprocess operation, the process is implemented on the thematic result which considers the occurrence probability of target pixels, without paying attention to the spectrum of the center pixel itself. The TSP considers both spectral and spatial information and realizes the balance to identify the paddy rice.

Accuracy analysis at different window sizes
In the previous analysis, 3 × 3 was initial window size to test the applicability of TSPSM. Also, various windows were conducted for analyzing the performance of TSPSM. The window size varied from 1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9 to 11 × 11. Among these, when window size was set as 1 × 1, there was only one pixel in the window. In this situation, we applied TSPSM without considering the surrounding pixels.
From Table 3, the results of Z-test demonstrate that there is a significant difference (Z > 1.96, at 95% confidence intervals) between TSPSM and PCC because PCC's results are affected by cloud contamination seriously. From individual scales, the Z-test between TSPSM and MV also reports the significant difference, ranging from 9.78 to 55.11. The results of coefficient of Kappa varied from 0.75 to 0.81 when window size ranged from 1 × 1 to 11 × 11, which are still higher than PCC and MV. The overall accuracy, as well as Kappa coefficient, shows a similar trend reaching the highest value at window size 3 × 3. This means window with 3 × 3 size is the best for this experiment. In Figure 5 (c), the scatter pixels of paddy rice are misclassified as the non-paddy rice because spectral heterogeneity occurs in the paddy rice field. This issue is gradually resolved with the window size growth which incorporating the surrounding probability to improve the confidence of paddy rice. Therefore, the user's and producer's accuracy reach the highest simultaneously when window size is 3 × 3, which obtains the balance of commission and omission error paddy rice related to non-paddy rice.

Effects of landscapes fragmentation on paddy rice identification
The degree of landscape fragmentation is a sensitive factor related to the performance of identification for land cover (Chen et al., 2012). Here, fragmentation is the breaking up of a habit or cover type into smaller, disconnected parcels (Turner, Gardner, & O'Neill, 2001). A higher degree of fragmentation implies much more fragmented and smaller parcels. Further analysis was carried out in different fragmented landscape scenarios to test the applicability of the proposed method. There were two subregions selected by visual analysis. This selection was executed referring to reference data and the corresponding high-resolution images. Quantitative accuracy assessment was conducted for these two subregions to analyze the performance of TSPSM.
The subregion with lower fragmentation (Lregion) ( Figure 6) and the other with relatively higher fragmentation (H-region) (Figure 7) are obviously different. L-region has large and compact fields of   paddy rice, whereas H-region is more fragmented with finer small non-paddy rice parcels. In Table 4, the trend of the overall accuracy for L-region initially goes up and down, peaking at 92.79% when window size is 7 × 7. From Figure 6 (b)-(c), ahead of the point of division, the spectral heterogeneity occurring in the paddy rice field dominates the results, and the paddy rice is feasible to be misclassified as non-paddy rice. This issue can be solved to an extent by integrating more spatial information with the window size increasing. However, more commission errors appear with larger window.
The result of overall accuracy shows a similar trend in H-region, but peaking at 87.42% at 3 × 3 window size scale, which is affected more easily by changing the size of window.
The overall accuracies in H-region are higher than in L-region when window size is 1 × 1( Table 4). The reason is that the spectrum of paddy rice is more heterogeneous in L-region, mainly from the discrepancy of crop growth. It can clearly find spectral heterogeneity of paddy rice at L-region. At the initially small window, the traits of spectral mainly contribute to the results, exceeding to the spatial information. For other window from 3 × 3 to 11 × 11, the overall accuracies in H-region show opposite trend, lower than in L-region. These results caused from the spatial scales have made more contribution to the paddy rice results to lessen the influence of spectral heterogeneity, generally as "salt-andpepper" phenomenon.
In summary, window size is a core factor for implementing the proposed method. However, to determine suitable window size for mapping paddy rice, more efforts are needed to focus on the discrepancy of paddy rice landscape heterogeneity.

Conclusions
This paper presented an easy and feasible method, TSPSM, to identify the paddy rice. An index, TSP, as a core element for TSPSM was defined by incorporating spatial and temporal information. Compared with two conventional methods, PCC and MV, TSPSM presented better performance, which demonstrated that the proposed model can avoid the impacts of cloud contamination easily. Several conclusions are drawn.
First, the algorithm of TSP in TSPSM is a sensitive indicator to improve accuracy of mapping paddy rice because that the model synthesizes spatial and temporal information. According to the result of pairwise Z-test between TSPSM and PCC, TSPSM and MV, respectively, Z values were greater than 1.96. These results demonstrate that at 95% confidence interval level, the accuracies have a significant difference between the proposed model and PCC, the proposed model and MV, respectively. The proposed model at individual window size (1 × 1, 3 × 3, 5 × 5, 7 × 7, 9 × 9, 11 × 11) obtained the map of paddy rice with overall accuracy ranging from 87.40% to 90.47%, and Kappa value varying from 0.75 to 0.81. These accuracies for mapping were higher than the other two methods. These results support the idea that the proposed model can solve the issues of cloud contamination and "salt and pepper". Moreover, analysis of the contribution of cloud contaminated images demonstrates that clear pixels in images can provide valuable information for mapping.
Secondly, window size is a sensitive factor for implementing TSPSM. The overall accuracy, Kappa value and user's accuracy decreased with increasing window size, while producer's accuracy increased. These results are mainly caused by spectral heterogeneity. It also was found that 3 × 3 was the optimal choice for this experiment. The reason is that using this scale acquired the highest overall accuracy and kappa coefficient, while balanced the commission and omission error for mapping. Moreover, the proposed method can reduce effects of "salt and pepper" at some degree due to the incorporation of spatial information of probabilities of surrounding pixels.
Furthermore, character of the landscape (i.e. fragmentation, heterogeneity) is another sensitive factor to affect the performance of TSPSM. For instance, the results in areas with different fragmentation showed that optimal window size should be defined separately. In an area with higher fragmentation, smaller window should be set for TSPSM, while larger window was better in the area with lower fragmentation.

Future work
Although TSPSM can reduce the constraints of the cloud contamination and "salt and pepper", several issues still need to be concerned. First of all, the result reveals that smaller window is suitable for areas with higher fragmentation, while larger window is appropriate for areas with lower fragmentation. It is necessary to define a quantitative index representing the fragmented degree, and analyze the relationship between the index and paddy rice identification for future work.
Another limitation is the proposed method involves more complexity if applying it for multitarget crops identification with different or similar crop calendar. To identify multi-target crops, individual type of crop at a time is extracted in TSPSM, and all crops are extracted iteratively. Eventually, the principal is that pixel with the highest probability belonging to certain crop class is assigned to that one, comparing the values of P ij among all the crops.
With increasing acquisition of different source of satellite sensors, other sensors of the satellite, such as SPOT 5, Sentinel-2 and so on, will provide more choices to calculate TSP by integrating images of one or more sensors. Currently, Google Earth Engine has become a powerful platform, which allows users to exploit the image cube to create a cloud-free image to monitor land surface timely and accurately (Patel et al., 2015). TSPSM can be realized in this platform for large-scale paddy rice identification. In addition, classification and regression tree is another tool for crop mapping by using multi-temporal remote sensing images (Lawrence & Wright, 2001;Shao & Lunetta, 2012). Further study will focus on how to integrate the posterior probability, multi-temporal spectral and phenology information to eliminate the effect of cloud contamination accurately.