Mapping rapeseed in China during 2017-2021 using Sentinel data: an automated approach integrating rule-based sample generation and a one-class classifier (RSG-OC)

ABSTRACT

Rapeseed mapping is important for national food security and government regulation of land use. Various methods, including empirical index-based and machine learning-based methods, have been developed to identify rapeseed using remote sensing. Empirical index-based methods commonly employed empirically designed indices to enhance rapeseed’s bright yellow spectral feature during the flowering period, which is straightforward to implement. Unfortunately, the heavy cloud cover in the flowering period of China would lead to serious omission errors; and the required flowering period varies spatially and yearly, which often cannot be acquired accurately. Machine learning-based methods mitigate the reliance on clear observations during the flowering period by inputting all-season imagery to adaptively learn features. However, it is difficult to collect sufficient samples across all of China considering the large intraclass variation in both land cover types of rapeseed and non-rapeseed. This study proposed an automated rapeseed mapping approach integrating rule-based sample generation and a one-class classifier (RSG-OC) to overcome the shortcomings of the two types of methods. First, a set of sample selection rules based on empirical indices of rapeseed were developed to automatically generate samples in cloud-free pixels during the predicted flowering period throughout China. Second, all available features composited based on the rapeseed phenological calendar were used for classification to eliminate the phenology differences in different regions. Third, a specific sample augmentation that removes the observation in the flowering period was employed to improve the generalization to the pixels without cloud-free observation in the flowering period. Finally, to avoid the need for diverse samples of nonrapeseed classes, a typical one-class classifier, positive unlabeled learning implemented by random forest (PUL-RF) trained by the generated samples, was applied to map rapeseed. With the proposed method, China rapeseed was mapped at 20 m resolution during 2017–2021 based on the Google Earth Engine (GEE). Validation on six typical rapeseed planting areas demonstrates that RSG-OC achieves an average accuracy of 94.90%. In comparison, the average accuracy of the other methods ranged from 83.33% to 88.25%, all of which were poorer than the proposed method. Additional experiments show that the performance of RSG-OC was not sensitive to cloud contamination, inaccurate predicted flowering time and the threshold of sample selection rule. These results indicate that the rapeseed maps produced in China are overall reliable and that the proposed method is an effective and robust method for annual rapeseed mapping across China.

KEYWORDS:

• Rapeseed mapping
• rule-based sample generation
• sample augmentation
• positive and unlabeled classification
• random forest

1. Introduction

Rapeseed (Brassica napus, Brassica campestris L) is one of the most widely planted oil crops, serving as an important source of edible oil, animal feed, and green biofuel (Raboanatahiry et al. 2021; Firrisa, van Duren, and Voinov 2013; van Duren et al. 2015; Zhang and He 2013). Global rapeseed production has undergone sustained growth over the past 20 years, becoming the second-most produced oilseed behind soybeans (FAO 2021, 2022; Carré and Pouzet 2014). And China is one of the largest producers and consumers of rapeseed in the world. (FAO 2022; USDA 2022; Hu et al. 2017; Tian et al. 2018). With the development of China’s livestock industry and renewable energy sector, the demand for rapeseed has increased in recent years (Tian et al. 2021; Yin and Wang 2012). However, rapeseed production in China is constrained by limited cultivated land resources (Bonjean, Dequidt, and Sang 2016; Fu et al. 2016) Therefore, monitoring the extent of rapeseed in China is important for macro policy regulation of land use and national food security (Tao et al. 2019; Zang et al. 2020).

Remote sensing is an effective tool for large-extent crop mapping (Dong et al. 2016; Weiss, Jacob, and Duveiller 2020; Gumma et al. 2018; Huang et al. 2022; Mansaray et al. 2019; Zhou et al. 2019). Several methods have been proposed to identify rapeseed using remote sensing, which can be grouped into empirical index-based and machine learning-based methods (Ashourloo et al. 2019; Han et al. 2021; Tao et al. 2019; Zang et al. 2020; Sulik and Long 2020; Zhang, Liu, and Zhang 2022). Empirical index-based methods commonly utilize the spectral features of flowering, i.e. the distinctive bright yellow color of rapeseed flowers (Fang et al. 2016; Sulik and Long 2015). Thus, several yellowness indices, which enhance the spectral feature of increasing reflectance of red and green bands for the rapeseed flower, were developed to identify rapeseed (Ashourloo et al. 2019; Sulik and Long 2016, 2020; Zang et al. 2020). For example, Ashourloo utilized a Canola Index time-series, computed as (Green+Red)×Nir, to classify rapeseed in 4 Iran and America test areas (Ashourloo et al. 2019). Han et al. developed an empirical index-based method to map rapeseed in 33 countries by combining the high normalized difference yellow index (NDYI) values on Sentinel-2 images during the flowering period and the high VH values on Sentinel-1 images during the podding period (Han et al. 2021). Sulik & Long combined NDYI and visible atmospherically resistant index (VARI) to identify flowering transition date of rapeseed (Sulik and Long 2020). Zang et al. developed an enhanced area yellow index (EAYI) by combining the features of the different yellow index (DYI) and NDVI during flowering period and mapped rapeseed flower coverage in five typical rapeseed planting areas in China (Zang et al. 2020). However, these methods heavily rely on cloud-free images acquired during the flowering period. Unfortunately, cloud-free observations during the flowering period in China are often lacking because rapeseed often flowers in the rainy season (Zang et al. 2020). This issue would lead to serious omission errors for empirical index-based methods. Furthermore, the flowering period is needed to obtain flowering images, which is difficult to be predicted accurately because rapeseed phenology varies spatially and yearly. Machine learning-based methods can address this issue because both flowering and nonflowering features can be learned from the training samples (Meng et al. 2020; Tao et al. 2019). For example, Tao et al. utilized an artificial neuron network to classify rapeseed on the Jianghan Plain in China (Tao et al. 2019). Meng et al. used random forest to classify winter wheat and rapeseed, which showed that postflowering features are also important for classification (Meng et al. 2020). However, collecting sufficient training samples required in machine learning-based methods is laborious and expensive (Ashourloo et al. 2019; Zhong, Gong, and Biging 2014). To overcome the above shortcomings, a hybrid method combining the advantages of empirical index-based and machine learning-based methods named seamless automated rapeseed mapping (SARM) was proposed very recently (Zhang, Liu, and Zhang 2022). First, a temporal difference in the winter rapeseed index (WRI) during the flowering season and before-flowering season, which empirically captures the unique phenological feature of rapeseed, was employed to automatically select rapeseed samples from cloud-free pixels. Then, based on these samples, multiple random forest classifiers were trained to classify each temporal image of the time series, and the classified results were integrated by voting on the classification probabilities of all single temporal images. This method works effectively in mapping rapeseed in the Yangtze River Basin. However, it still cannot meet the needs of large-scale and long-term rapeseed mapping in China for the following reasons. First, rapeseed is widely planted with a large phenological and spectral variation in China (Hu et al. 2017; Zang et al. 2020), and the non-rapeseed classes are more diverse than in the Yangtze River Basin. The WRI difference was designed for distinguishing winter rapeseed and winter wheat, which is inapplicable for other regions with more complicated rapeseed and non-rapeseed land cover types. Second, the SARM classified each temporal image independently to avoid the classifier learning the WRI difference itself as determining features; however, it would neglect the relationship among multitemporal features.

To address the above issue, a method integrating rule-based automatic sample generation and one-class classification (RSG-OC) was proposed to map rapeseed in China. A set of rules was developed based on empirical indices of flowering features to automatically select rapeseed samples from cloud-free pixels during the predicted flowering period. The samples were then augmented to improve the generalization to the pixels without cloud-free flowering features. Finally, a typical one-class classifier, positive unlabeled learning (PUL) implemented by random forest (PUL-RF), which only requires rapeseed samples and unlabeled samples randomly generated from all pixels, was applied to avoid the need for samples of diverse non-rapeseed classes (Chen et al. 2016; Lu and Wang 2021).

With the proposed method, rapeseed was mapped in China at 20 m resolution for 2017–2021 based on the Google Earth Engine (GEE) (Gorelick et al. 2017). Quantitative accuracy evaluation was conducted in 6 typical rapeseed planting areas, and 3 methods of rapeseed mapping were compared with the proposed method.

The rest of the paper is organized as follows. The materials and methodology are described in Sections 2 and 3, respectively. The results and discussion are presented in Section 4 and Section 5. Section 6 concludes the study.

2. Data and study area

2.1. Study area

The study aims to map rapeseed across China (73.5°E-135°E, 4°N-53.5°N) (Figure 1(a)). The main planting area of rapeseed is distributed widely in China, spanning from most southern in Yunnan (24.5°N, 101.5°E) to most northern in Hulun Buir (119.77°E, 49.21°N) and from most eastern in Jiangsu (119.8°N, 32.9°E) to most western in Yili (84.99°E, 44.45°N) (Wang, Guan, and Zhang 2007; Zang et al. 2020), with various climate types, including subtropical monsoon, temperate monsoon, and temperate continental climates. As a result, rapeseed shows a large difference in phenology and cultivar (Qian et al. 2006). Based on the field phenological observation records (Zhang, Liu, and Zhang 2022; Zhang et al. 2022; Li et al. 2021), the rapeseed phenological calendar in four typical rapeseed planting sites, Qujing (site 1), Jingzhou (site 2), Haibei (site 3), and Hulun Buir (site 4), is presented in Figure 1(b). The phenological stages of rapeseed include seeding, emergence, budding, flowering, podding, mature, and harvest. All of these phenological stages are delayed from warm areas to cold areas. For example, the flowering dates ranged from February in the southernmost region (site 1) to August in the northmost region (site 4). Moreover, the length of phenological stages could differ among different cultivars. Especially, the winter rapeseed sowed from October to November in the southern region has a much longer emergence stage than spring rapeseed sowed from April to July in the northern region.

Figure 1. (a) Study area and distribution of rapeseed phenological calendar. (b) Phenological calendar of rapeseed in four typical sites. “E, “M,” and “L” represent the early, middle, and late periods of the month, respectively.

2.2. Data

2.2.1. Satellite data

Based on the growth calendar of rapeseed, Sentinel-1 and Sentinel-2 imagery during October 2016-October 2021 accessed on the GEE platform were used in this study.

Sentinel-1 is a C-band Synthetic Aperture Radar (SAR) satellite with dual polarization (VV and VH). To match Sentinel-2, we used the ground-ranging detection (GRD) product in IW mode, of which the spatial resolution is 10 m and the temporal resolution is 12 days (Torres et al. 2012). Previous studies have shown that the Sentinel-1 data has specific responses to the canopy structure of rapeseed, which helps to distinguish rapeseed from other crops (D’Andrimont et al. 2020; Veloso et al. 2017; Han, Zhang, and Cao 2021). Specifically, the VV scattered signal drops in the flowering season, and the VH scattered signal increases in the podding stage. Therefore, Sentinel-1 data plays an important role in identifying rapeseed in the case of cloudy weather.

Sentinel-2 is a multispectral optical imaging satellite with 10–60 m spatial resolution, 5–10 days temporal resolution, and 12 bands (Drusch et al. 2012). Both the top-of-atmosphere reflectance products at the L1C level (hereafter referred to as L1C data) and ground reflectance products at the L2A level (hereafter referred to as L2A data) were used in this study. The L2A data are used for sample generation with empirically designed rules because the thresholds in the rules could be more stable if the atmospheric effects are corrected (Zang et al. 2020). However, the availability of L2A data on GEE for most areas of China (December 2018-present) was shorter than the availability of L1C data (August 2015-present). Therefore, L1C data-derived features were used for the classification of rapeseed considering that machine-learning techniques are not that sensitive to atmospheric effects (Jin et al. 2019; Pott et al. 2021). Moreover, the Sentinel-2 Cloud Probability product was used to remove clouds (ESA 2020).

2.2.2. Land cover data

We focused our study on cultivated land considering that most rapeseed is planted in cultivated land except for a negligibly small amount of rapeseed planted in urban green spaces. Such exclusion of non-cultivated land helps to save computational resources and avoid the interference of the non-cultivated land in rapeseed mapping (Konduri et al. 2020; Zhang et al. 2017). To avoid the omission error of available cultivated land cover datasets, we integrated three widely used global land cover products, GlobeLand30–2020 (Chen et al. 2015), ESA World Cover V100 (Zanaga et al. 2021), and ESRI 2020 Global Land Use Land Cover (Karra et al. 2021), to take the union of their cropland layers to derive the regions to map.

2.2.3. Reference data

Ground truth data for evaluating the mapping results were obtained in six typical rapeseed planting regions by field survey or interpretation from very high-resolution images. As shown in Figures 2, 3 field surveys were conducted in Qujing city (Figure 2(d)), Sichuan Province (Figure 2(c)), and Qinghai Province (Figure 2(b)) during the flowering season (i.e. February, March, and July, respectively) of 2021. The GPS coordinates of crop parcels were recorded using mobile GPS devices. A total of 576 rapeseed and 494 non-rapeseed parcels were collected. These parcels were converted to binary images with 20 m resolution, and 105,179 rapeseed samples and 28,663 non-rapeseed samples were obtained. Two high-resolution Google Earth images in Taizhou (Figure 2(e)) and Lincang (Figure 2(g)) at a 0.5 m resolution and an unmanned aerial vehicle (UAV) image in Jingzhou (Figure 2(F)) at a 0.07 m resolution were collected during the rapeseed flowering period in 2020. All the high-resolution images were manually interpreted as rapeseed and non-rapeseed, resampled to 20 m resolution, and randomly selected as validation samples. A total of 13,744 rapeseed samples and 14,979 non-rapeseed samples were selected from the images as ground truth data for accuracy evaluation.

Figure 2. The locations of ground truth data. (b) to (d) are the distributions of validated parcels obtained by field campaigns in Qinghai, Sichuan, and Qujing. (e) to (g) are rapeseed maps interpreted from very high-resolution images acquired during the flowering period in Taizhou, Jingzhou, and Lincang.

Figure 3. Flowchart of the proposed method. (a) Data preprocessing, (b) sample generation, (c) feature composition, (d) sample augmentation, and (e) PUL-RF classification.

Census data of rapeseed area of 32 provinces in 2017–2020 were acquired from the China statistical yearbooks released by the National Bureau of Statistics of China (http://www.stats.gov.cn/tjsj./ndsj/) to compare with the mapped rapeseed area at the provincial level.

3. Methodology

The proposed method includes five main stages: data preprocessing, sample generation, feature composition, sample augmentation, and PUL-RF classification (Figure 3).

3.1. Data preprocessing

As mentioned above, the input data for classification consist of Sentinel-1, Sentinel-2, and cropland data. The preprocessing of the Sentinel-1 data includes thermal noise removal, radiometric calibration, and topographic correction, implemented by a GEE toolbox (Thorp and Drajat 2021). The preprocessing of the Sentinel-2 data includes cloud filtering, band selection, index calculation, compositing, and interpolation. First, the Sentinel-2 images with cloud cover below 95% were selected, and the pixels with cloud probability above 50% were masked out using the ESA Sentinel-2 cloud probability product (ESA 2020). Second, the blue (B2), green (B3), red (B4), and near-infrared (B8) bands with a 10-m spatial resolution and the four red edge (B5 to B7, and B8A) and two shortwave infrared (B11 and B12) bands with a 20-m resolution were used. Third, the commonly used NDVI, DYI (Zang et al. 2020), NDYI (Sulik and Long 2016), and LSWI (Chandrasekar et al. 2010) were calculated for each pixel according to Eq. (1)(1) $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{B}8-\mathrm{B}4}{\mathrm{B}8+\mathrm{B}4}$(1) to (4)(4) $\mathrm{L}\mathrm{S}\mathrm{W}\mathrm{I}=\frac{\mathrm{B}8-\mathrm{B}11}{\mathrm{B}8+\mathrm{B}11}$(4) ,

(1) $\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{B}8-\mathrm{B}4}{\mathrm{B}8+\mathrm{B}4}$(1)

(2) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{Y}\mathrm{I}=\mathrm{B}3-\mathrm{B}2\end{array}\end{array}$(2)

(3) $\mathrm{N}\mathrm{D}\mathrm{Y}\mathrm{I}=\frac{\mathrm{B}3-\mathrm{B}2}{\mathrm{B}3+\mathrm{B}2}$(3)

(4) $\mathrm{L}\mathrm{S}\mathrm{W}\mathrm{I}=\frac{\mathrm{B}8-\mathrm{B}11}{\mathrm{B}8+\mathrm{B}11}$(4)

The Sentinel-2 time series of each band and index were then composited by the 10-day median synthesis. Finally, the null values in the time series are linearly interpolated using the valid values from the adjacent times.

3.2. Generation of rapeseed samples

Sufficient samples are the basis for machine learning-based methods to perform large-area mapping (Foody and Mathur 2004; Maxwell, Warner, and Fang 2018). However, sample acquisition usually requires costly manual selection or field investigation (Friedl et al. 2002; Tian et al. 2020; Xiong et al. 2017). Considering that the yellow rapeseed flower is a distinctive feature compared to other crops (Fang et al. 2016; Sulik and Long 2020), a series of rules can be developed to select rapeseed samples automatically from cloud-free pixels, which only requires at least one clear observation during the flowering period (Zhang, Liu, and Zhang 2022).

The flowering period of each area was first determined by geolocation using a multilinear regression model according to (Zang et al. 2020). The model (Equation 5(5) $\begin{array}{c}\begin{array}{c}T=7.07\times \mathrm{L}\mathrm{a}\mathrm{t}+1.508\times \mathrm{L}\mathrm{o}\mathrm{n}+0.03\times \mathrm{A}\mathrm{l}\mathrm{t}-318.11\end{array}\end{array}$(5) ) was established between the multiyear average flowering date (T) recorded by agrometeorological stations and longitude (Lon), latitude (Lat), and altitude (Alt) (R² close to 0.9).

(5) $\begin{array}{c}\begin{array}{c}T=7.07\times \mathrm{L}\mathrm{a}\mathrm{t}+1.508\times \mathrm{L}\mathrm{o}\mathrm{n}+0.03\times \mathrm{A}\mathrm{l}\mathrm{t}-318.11\end{array}\end{array}$(5)

Previous studies have examined several temporal features of rapeseed (Ashourloo et al. 2019; Han et al. 2021; Zang et al. 2020). Rapeseed flowers can increase the DYI and decrease the NDVI (Shen et al. 2014, 2010; Zang et al. 2020). Thus, we analyzed the DYI and NDVI of rapeseed on the Sentinel-2 time series. In addition, we also analyzed the LSWI, which can reflect the water content of rapeseed due to the high water content of rapeseed (Chandrasekar et al. 2010; Mandal et al. 2020). Comparing the time series of DYI, NDVI, and LSWI of rapeseed and other crops in two typical rapeseed growing areas, Jingzhou and Qinghai, we found several temporal features that can be used to develop the rules for sample selection (Figure 4). During the predicted flowering period, the DYI of rapeseed peaked, while the NDVI decreased to a local minimum because of the yellow flowers (Zang et al. 2020). In contrast, the DYI values of other crops remained low, and there was no such NDVI valley for them during the period. The LSWI of rapeseed is high during the flowering period, which will be used to remove interference pixels.

Figure 4. Temporal features of rapeseed and other crops. One parcel of each crop was selected for analysis. (a) Winter rapeseed, (b) winter wheat, and (c) water field in Jingzhou; (d) spring rapeseed, (e) barley, and (f) corn in Qinghai. The red dashed line indicates that the time when the DYI reaches a maximum is equal to the time when the NDVI reaches a minimum; the black dashed line is the threshold of the NDVI and LSWI.

The DYI peak in the Moderate Resolution Imaging Spectroradiometer (MODIS) time series has been demonstrated to be a consistent feature of rapeseed in different regions of China by subtracting the base DYI at nonflowering periods from the peak DYI (Zang et al. 2020). We adapted the idea for the Sentinel-2 time series and derived the differential-DYI (DDYI) with a modified base DYI (Eq. (7)(7) $\begin{array}{c}{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=\frac{{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\left[T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-60,T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-15\right]}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}+{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\left[T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}+15,T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}+60\right]}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}}{2}\end{array}$(7) ) due to the poorer data availability of Sentinel-2 compared with MODIS. DDYI, calculated by Eqs. 6(6) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}={\mathrm{D}\mathrm{Y}\mathrm{I}}_{\left[T-30,T+30\right]}^{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}\end{array}\end{array}$(6) to 8(8) $\begin{array}{c}\begin{array}{c}{T}_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}=argmax\mathrm{D}\mathrm{Y}{\mathrm{I}}_{\left[T-30,T+30\right]}\end{array}\end{array}$(8) , is a global adaptable flowering index across China and is used as a main feature for selecting rapeseed samples.

(6) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}={\mathrm{D}\mathrm{Y}\mathrm{I}}_{\left[T-30,T+30\right]}^{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}\end{array}\end{array}$(6)

(7) $\begin{array}{c}{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=\frac{{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\left[T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-60,T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-15\right]}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}+{\mathrm{D}\mathrm{Y}\mathrm{I}}_{\left[T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}+15,T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}+60\right]}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}}{2}\end{array}$(7)

(8) $\begin{array}{c}\begin{array}{c}{T}_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}=argmax\mathrm{D}\mathrm{Y}{\mathrm{I}}_{\left[T-30,T+30\right]}\end{array}\end{array}$(8)

where $T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}$ is the time when DYI achieves the maximum during the period $\left[T-30,T+30\right]$, and ${\mathrm{D}\mathrm{Y}\mathrm{I}}_{\mathrm{b}\mathrm{a}\mathrm{s}\mathrm{e}}^{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$ is the basic DYI during the nonflowering period $\left[T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-60,T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-15\right]$ and $\left[T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}+15,T_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}+60\right]$. Note that a time window of [T-30, T + 30], which was longer than the flowering period (20–40 days in usual), was applied here to ensure that the flowering image can be found within the period. A pixel with a higher DDYI indicates that the pixel is more likely to be flowering rapeseed. To select rapeseed samples as accurately as possible, an empirically strict threshold is applied to DDYI:

(9) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}>0.04\end{array}\end{array}$(9)

Some non-rapeseed crops could also have a high DDYI if their greening-up stage overlaps with the flowering period of rapeseed. In this case, the reflectance in the green band increases dramatically, leading to a high DYI value. Fortunately, the NDVI valley induced by yellow flowers could be used for distinguishing rapeseed from such crop types (Zang et al. 2020). Therefore, these nonrapeseed crop pixels can be excluded by equation (10)(10) $argmax\mathrm{D}\mathrm{Y}\mathrm{I}{ }_{\left[T-30,T+30\right]}=argmin\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}{ }_{\left[T-30,T+30\right]}$(10) :

(10) $argmax\mathrm{D}\mathrm{Y}\mathrm{I}{ }_{\left[T-30,T+30\right]}=argmin\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}{ }_{\left[T-30,T+30\right]}$(10)

Which means that the observed DYI peak should correspond to an NDVI valley for the rapeseed pixels (Zang et al. 2020). Equation 9(9) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}>0.04\end{array}\end{array}$(9) and are the main rules for selecting rapeseed samples across all of China.

However, some occasional non-rapeseed pixels with undesired high DDYI values remain, including rainbowlike clouds (Frantz et al. 2018; Han et al. 2021) and blue impervious surfaces (Figure S1). Therefore, two additional criteria Equation (11)(11) $\begin{array}{c}\begin{array}{c}\mathrm{L}\mathrm{S}\mathrm{W}{\mathrm{I}}_{T\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}>0.4\end{array}\end{array}$(11) and Equation (12)(12) $\begin{array}{c}\begin{array}{c}\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{T\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}>0.4\end{array}\end{array}$(12) were used to remove such pixels from the samples selected by the main rules of Eq. (9)(9) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}>0.04\end{array}\end{array}$(9) and (10)(10) $argmax\mathrm{D}\mathrm{Y}\mathrm{I}{ }_{\left[T-30,T+30\right]}=argmin\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}{ }_{\left[T-30,T+30\right]}$(10) :

(11) $\begin{array}{c}\begin{array}{c}\mathrm{L}\mathrm{S}\mathrm{W}{\mathrm{I}}_{T\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}>0.4\end{array}\end{array}$(11)

(12) $\begin{array}{c}\begin{array}{c}\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{T\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}>0.4\end{array}\end{array}$(12)

As such types of nonrapeseed pixels have much lower LSWI or NDVI values than rapeseed pixels, the thresholds are relatively easy to determine. A more detailed description of these two additional criteria is included in the supplementary materials.

Finally, with the proposed rules (Equation 9(9) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}>0.04\end{array}\end{array}$(9) -Equation 12(12) $\begin{array}{c}\begin{array}{c}\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{T\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}>0.4\end{array}\end{array}$(12) ), sufficient rapeseed samples from 2019 to 2021 were generated across China based on the Sentinel-2 L2A data (Figure 6). Meanwhile, the same number of unlabeled samples were also randomly generated from all cultivated land pixels, which is also an essential input of the one-class PUL classifier that will be introduced in Section 3.5. Then, the spatial locations of these generated samples across different years were merged to generate a training sample set based on Sentinel-2 L1C data.

3.3. Feature composition based on phenology calendar

Variation in rapeseed phenology in different regions of China and availability of cloud-free satellite data challenges classification. To alleviate their negative effects, the phenology of rapeseed was first divided into three periods according to the flowering date (T): preflowering ([T-45, T-15)), flowering ([T-15, T + 15)) and postflowering ([T + 15, T + 45]). The three periods could roughly correspond to the budding, flowering, and podding stages, which all last approximately 30 days in different geographic locations (Figure 1). Then, the median composition was employed for each time series at each stage (Phan, Kuch, and Lehnert 2020), and a total of 48 composited feature images were generated (three stacks of 16 features from Sentinel-1 and −2) (Figure 5(a)) for the consequent classifier training and classification.

Figure 5. Schematic diagram of (a) feature composition and (b) sample augmentation.

3.4. Sample augmentation

Note that the above generation of samples and the corresponding composited features were mainly based on the flowering signal. These samples thus would fail to represent the rapeseed pixels without adequately clear observation of flowering due to cloud blocking, because the interpolation step in Section 3.1 could produce distorted values in case of continuous cloudy weather. To improve the representativeness of the rapeseed samples, the samples were augmented by purposefully replacing the original cloud-free observations in the flowering period ([T-15, T + 15)) with the linearly interpolated values from the valid observations in the pre and postflowering periods (i.e. [T-45, T-15) and, [T + 15, T + 45]) (Figure 5), thus simulating the poor observation condition. Note that the simulation was employed for all the features except VV and VH because Sentinel-1 is also effective under cloudy weather (D’Andrimont et al. 2020; Mercier et al. 2020). Sample augmentation was applied to both the rapeseed samples and the unlabeled samples. The expanded sample set merged the original and augmented sample sets in equal proportions. With the expanded sample set including adequate clear and unclear samples, machine-learning method is expected to learn and adaptively use both flowering and non-flowering features for rapeseed classification.

3.5. Application of PUL-RF

Considering the diversity of the nonrapeseed classes, the one-class random forest classifier implemented by the PUL strategy (PUL-RF), which only needs the positive samples of the target land cover type and the unlabeled samples randomly generated from all pixels in cultivated land layers (Chen et al. 2016; Elkan and Noto 2008; Li, Guo, and Elkan 2010), was employed to map rapeseed. Here, the unlabeled samples include both the positive and negative samples, thus, were generated automatically without laborious labeling efforts.

PUL is designed to transform the traditional binary classifier into a one-class classifier based on the Bayesian rule and has been proven effective for mapping other land cover types (Lei et al. 2021; Lu and Wang 2021). It includes two steps. First, a binary random forest classifier is trained to calculate the probability that a pixel is labeled. Then, the probability that a pixel belongs to rapeseed is adjusted by a factor (c) that represents the labeling probability of positive pixels

(13) $\begin{array}{c}f\left(y=1|x\right)=\frac{p\left(s=1|x\right)}{c}\end{array}$(13)

where $f\left(y=1|x\right)$ represents the probability that a pixel is rapeseed, x is the feature of the pixel and y∈{−1,1} represents the land cover types (non-rapeseed and rapeseed). $p\left(s=1|x\right)$ is the probability that a pixel is labeled, and s∈{0,1} represents whether it is labeled or not. The adjusting factor (c) was estimated through an independent validation set (V).

(14) $\begin{array}{c}c=\frac{1}{n}\sum _{x\in V}p\left(s=1|x\right)\end{array}$(14)

Here, the independent validation set (V) with n pixels was divided from the expanded rapeseed samples. This validation set (V) is excluded in the training phase but used for estimating the adjusting factor (c) and determining the hyperparameter of the random forest. Finally, a single PUL-RF classifier trained with the expanded sample set was used to map rapeseed during 2017–2021.

3.6. Accuracy assessment

To evaluate the mapping accuracy, the overall accuracy (OA), F1-score, producer accuracy (PA), and user accuracy (UA) were calculated using the ground truth data in the six experimental regions (Section 2.2.3). The mapping accuracy and spatial details of RSG-OC were also compared with those of the NDYI-based method (Han et al. 2021), the CI-based method (Ashourloo et al. 2019), and SARM (Zhang, Liu, and Zhang 2022).

4. Results

4.1. Spatial distribution of rule-based generated samples

Based on the sample generation rules and sample augmentation, a total of 61,464 samples were generated, including 32,090 rapeseed samples and 29,374 unlabeled samples. The number of rapeseed samples for 2019, 2020, and 2021 was 7944, 12164, and 12,000, respectively. Generally, the spatial distribution of the generated rapeseed samples is consistent with the main rapeseed production regions, and it varies slightly from year to year due to the cloud cover variation during the flowering period (Figure 6).

Figure 6. (A) Spatial distribution of generated rapeseed samples; (b) Number of generated rapeseed samples for 2019–2021.

4.2. Map of rapeseed in China

The generated expanded sample set is divided into a training set and a validation set (V) by a ratio of 7:3 to train the PUL-RF. The tree number of the RF was set as 225 based on the highest OA on the validation set. The parameter c in PUL was derived from Eq.14(14) $\begin{array}{c}c=\frac{1}{n}\sum _{x\in V}p\left(s=1|x\right)\end{array}$(14) , equaling to 0.87 in this study. The other parameters of the PUL-RF, such as the number of variables per split, were set to the default values in the GEE function of smileRandomForest. Using the trained PUL-RF, 20 m rapeseed maps of China for 2017–2021 were produced. The data link is provided in the data availability section. Consistent with prior knowledge (Wang, Guan, and Zhang 2007; Bonjean, Dequidt, and Sang 2016), rapeseed was widely cultivated in China in several regions, including the middle-lower Yangtze Plain, Sichuan Basin, Yunnan-Guizhou Plateau, northeastern Qinghai Tibet Plateau, and northeastern and northwestern North arid and semiarid regions (Figure 7(a)). The mapped area of rapeseed in China during 2017–2021 is 2.46, 2.34, 2.35, 2.81, and 2.75 million hectares respectively (Figure 7(b)). Although the mapped area seems to have a similar temporal fluctuation with the census data, there is a large underestimation of the mapped area compared to the census data, indicating a nonnegligible difference between the two data sources.

Figure 7. (a) Chinese rapeseed coverage map at 500 m resolution in 2020 (the original 20 m resolution was resampled to 500 m resolution for a better illustration in the thumbnail figure); (b) China rapeseed mapping area and census area, 2017–2021 (Note: census area in 2021 is unavailable).

To further illustrate the reasonability of rapeseed change in the produced annual maps, Figure 8 compared the temporal change in the rapeseed map over 2017–2021 in detail in four rapeseed planting sites. The Sentinel-2 images at the flowering period are also presented for visual comparison if they are cloud-free. Otherwise, the most cloud-free and temporally close images were used, and fortunately, rapeseed can still be identified on most of these images because the green color of rapeseed is subtly different from that of other vegetation. Visually, the maps match the rapeseed parcels every year and present clear details.

Figure 8. Detailed comparison of multiyear mapping results in 4 typical rapeseed planting areas. (a) Jingzhou, (b) Qujing, (c) Hulun Buir, (d) Qinghai. The first line in each region is the Sentinel-2 true color images, and the second line of each area is the mapped rapeseed. Images with ‘*’ after the annotated date were acquired in the nonflowering period.

4.3. Visual comparison with other methods

Figure 9 compares the Sentinel-2 images during the flowering period, the ground truth (ground truth rapeseed maps interpreted from very high-resolution images in Jingzhou, Taizhou, and Lincang and crop parcels acquired from field surveys in Sichuan, Qujing, and Qinghai from Section 2.2.1), and the mapping results of the NDYI-based method, the CI-based method, SARM, and the proposed RSG-OC at six validation areas. It should be noted that owing to constant cloud cover, the compared Sentinel-2 images in Sichuan and Qujing are not appropriate (and the one is at the nonflowering period and the other is partly covered by clouds) but are already the best we could find. In general, the mapping results of RSG-OC were in good agreement with the spatial distribution at all six validation areas, whereas the other methods produced large classification errors in some areas. Specifically, the NDYI-based and CI-based methods have a large omission error in the Qujing and Sichuan because they rely heavily on cloud-free observations during the flowering period. The results of SARM are poorly consistent with the actual rapeseed distribution at the Qinghai and Sichuan areas because their sample selection rule is inapplicable at these two areas.

Figure 9. Comparison of the mapping details of the four methods in the validation areas (a) Jingzhou, (b) Taizhou, (c) Lincang, (d) Sichuan, (e) Qujing, and (f) Qinghai. The first column shows the RGB images of Sentinel-2 during the flowering periods in the validation areas except for Sichuan area, where the flowering images were totally cloud contaminated and the image acquired after flowering was shown instead. The second column shows the ground truth maps, which were interpreted from very high-resolution images for (a-c) and obtained from field-survey for (d-f). Columns 3 to 6 show the mapping details of CI, NDYI, SARM, and proposed RSG-OC respectively.

4.4. Quantitative evaluation

The quantitative assessment is summarized here using the ground truth dataset (Table 1) and the agricultural census area data (Figure 10). Among all of the compared methods, the proposed RSG-OC has the best performance with an average OA of 94.90% and an average F1-score of 93.50%. RSG-OC outperformed the other three methods in Taizhou, Sichuan, and Qinghai and achieved the second-highest OA in Jingzhou, Lincang, and Qujing. Although RSG-OC achieved stable performance in different areas, the PA in Sichuan was relatively low, indicating a large omission in this region. It is probably induced by very high cloud coverage (>54%) in Sichuan areas. SARM performed more poorly in Sichuan and Qinghai, with OAs of 81.58% and 56.59% and F1-scores of 61.75% and 67.70%, respectively, because the sample generation rules of SARM do not apply to some land cover in these two regions. The CI-based method has large omission errors in Sichuan and Qujing, with PAs of 42.69% and 43.69%, respectively, in these two regions. Similarly, the NDYI-based method has a poorer performance, with PAs of 0.35% in Sichuan and 22.13% in Qujing.

Figure 10. Comparison mapped rapeseed area with census area at the provincial level for (a)2017, (b)2018, (c)2019, (d)2020.

Table 1. Comparison of the accuracy of the four methods.

Region	Metrics	NDYI	CI	SARM	RSG-OC
Jingzhou	OA	91.52%	95.39%	88.76%	92.64%
	UA	96.58%	91.97%	87.19%	87.58%
	PA	86.05%	99.46%	90.85%	99.53%
	F1-score	91.01%	95.57%	88.98%	93.17%
Taizhou	OA	91.17%	93.72%	93.89%	94.43%
	UA	98.51%	94.30%	92.62%	93.37%
	PA	80.69%	90.88%	93.21%	93.71%
	F1-score	88.71%	92.56%	92.91%	93.54%
Lincang	OA	90.87%	86.44%	94.41%	93.47%
	UA	98.09%	97.77%	96.33%	97.86%
	PA	83.36%	74.58%	92.34%	88.89%
	F1-score	90.13%	84.62%	94.29%	93.16%
Sichuan	OA	83.24%	86.47%	81.58%	95.49%
	UA	57.14%	64.64%	47.38%	90.66%
	PA	0.35%	42.69%	88.63%	81.53%
	F1-score	0.70%	51.42%	61.75%	85.85%
Qujing	OA	64.79%	74.53%	98.91%	98.78%
	UA	100.00%	100.00%	100.00%	100.00%
	PA	22.13%	43.69%	97.59%	97.30%
	F1-score	36.24%	60.81%	98.78%	98.63%
Qinghai	OA	78.52%	92.93%	56.59%	94.58%
	UA	97.66%	99.17%	87.44%	99.69%
	PA	75.74%	92.19%	55.23%	93.72%
	F1-score	85.31%	95.55%	67.70%	96.61%
Average OA*		83.35%	88.25%	85.69%	94.90%
Average F1-Score*		65.35%	80.09%	84.07%	93.50%

*Average OA and average F1-score denote the averaged values of OA and F1-score of the six areas.

The RSG-OC mapped rapeseed area in each province was further compared with the official census area for 2017–2020. As shown in Figure 10, the mapped area is significantly correlated with the census area. However, the mapped areas are much lower than the census areas for most provinces; and the coefficients of determination (R²) are generally low, with some outliers like Jiangsu Province. The low consistency between the mapped area and census data might be induced by the uncertainties in both the remotely sensed maps and the census data (Liu et al. 2020), which should be further explored in the future.

5. Discussion

5.1. Dependence on clear observation in the flowering period

As aforementioned, both the index-based (CI and NDYI) and the machine learning-based methods (SARM and RSG-OC) require clear flowering observations but to different extents. The key difference here is the method’s sensitivity to the availability of clear observations in flowering season and should be investigated. Considering that RSG-OC employed sample augmentation and Sentinel-1 SAR data to address the issue of cloud contamination, two incomplete RSG-OC methods without sample augmentation and Sentinel-1 data were also compared to illustrate the importance of these two steps. To simulate the time-series data with weakening flowering signals induced by cloud contamination, we sequentially removed the observations during the flowering period with the top 33%, 66%, and 100% (namely, all removed) DYI values and checked how the performances of the compared methods worsened. As shown in Figure 11, both RSG-OC and SARM were stable against the decreasing flowering observations, whereas the average OA and F1-score of the CI-based and NDYI-based methods dropped largely as the clear flowering observations decreased, and they completely failed when there were no clear flowering observations. The performance of two incomplete RSG-OC methods without samples augmentation and Sentinel-1 data also decreased with decreasing clear observation in the flowering season. Especially, when none of the clear observations was available during the flowering season, two incomplete RSG-OC methods performed much poorer than standard RSG-OC. These results confirmed the contribution of both the Sentinel-1 SAR data and sample augmentation on resisting the issue of cloud contamination.

Figure 11. The sensitivity of the different methods to clear observation in the flowering period.

The role of sample augmentation could be further explained by the differences between the feature importance of RFs in standard RSG-OC and the incomplete RSG-OC without sample augmentation. As Figure 12 shows, the optical flowering features, e.g. LSWI_flower, NDYI_flower, and B11_flower, exhibited great importance for RSG-OC without sample augmentation. These optical features unfortunately would be largely weakened by cloud contamination, and thus is highly dependent on adequate clear observation. For standard RSG-OC, the importance of optical flowering features decreased, whereas SAR features and postflowering optical features became more important. Specifically, the VH_post, VH_flower, and VV_pre ranked 1^st,3^rd, and 5^th respectively, and B11_post, B5_post, and LSWI_post ranked 2^nd, 4^th and 6^th respectively. Therefore, the RSG-OC method is not that dependent on cloud-free observations during the flowering period.

Figure 12. Ranking of the feature importance in (a) RSG-OC without sample augmentation and (b) standard RSG-OC. The subscripts after the features indicate the period of feature composition. “pre” represents the “preflowering” period, “flower” represents the “flowering” period, and “post” represents the “postflowering” period.

5.2. Comparison with SARM

RSG-OC adapted the strategy of SARM that combines empirical index-based sample generation and machine-learning methods to map rapeseed in China. Compared with SARM, there are three main improvements of RSG-OC. First, an empirical rule set with a stricter “yellowness” feature compared with the WRI index was employed to automatically select rapeseed samples from cloud-free pixels during the flowering period. The stricter “yellowness” rule helps to better select correct rapeseed samples across the whole of China; however, the representativeness of the generated samples could be biased. Therefore, the second improvement compared to SARM is a specially designed sample augmentation that removes the flowering features of the generated rapeseed samples. Such sample augmentation helps to avoid the classifier learning the strict yellowness rule as a classification feature and improves the generalization to the pixels without cloud-free flowering features. Finally, a typical one-class classifier, positive unlabeled learning (PUL) implemented by random forest (PUL-RF), which only requires rapeseed samples, was applied to avoid the need for samples of diverse nonrapeseed classes (Chen et al. 2016; Lu and Wang 2021).

Especially, the rules employed in SARM method was designed for rapeseed in Yangtze River Basin (Zhang, Liu, and Zhang 2022), which was not applicable for other areas and lead to poor performance in Qinghai area (Table 1). The SARM selects rapeseed samples with following rules,

(15) $\begin{array}{c}\begin{array}{c}\mathrm{N}\mathrm{D}\mathrm{V}{\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}>0.5\end{array}\end{array}$(15)

(16) $\begin{array}{c}\begin{array}{c}\mathrm{W}\mathrm{R}{\mathrm{I}}_{\mathrm{b}\mathrm{f}}-\mathrm{W}\mathrm{R}{\mathrm{I}}_{\mathrm{f}}>0\end{array}\end{array}$(16)

(17) $\mathrm{W}\mathrm{R}\mathrm{I}=\frac{\mathrm{B}8-\mathrm{B}3}{\mathrm{B}8+\mathrm{B}3}\times \frac{\mathrm{B}2}{\mathrm{B}3+\mathrm{B}4}$(17)

where NDVI_max is the maximum NDVI during the growing season, WRI_bf is the WRI before flowering and WRI_f is the WRI at flowering. In the Yangtze River Basin, the above rules are adequate to distinguish the winter rapeseed from other winter crops. However, such rules failed to accurately select rapeseed samples when they were applied to a larger extent. Some other usual vegetation types outside the Yangtze River Basin were likely selected by the SARM rules. In Qinghai, for example, the barley and natural grasslands also satisfy the rules of Equation 16(16) $\begin{array}{c}\begin{array}{c}\mathrm{W}\mathrm{R}{\mathrm{I}}_{\mathrm{b}\mathrm{f}}-\mathrm{W}\mathrm{R}{\mathrm{I}}_{\mathrm{f}}>0\end{array}\end{array}$(16) (Figure 13(a,b)). In Sichuan, the evergreen forest also satisfies the rules of Equation 16(16) $\begin{array}{c}\begin{array}{c}\mathrm{W}\mathrm{R}{\mathrm{I}}_{\mathrm{b}\mathrm{f}}-\mathrm{W}\mathrm{R}{\mathrm{I}}_{\mathrm{f}}>0\end{array}\end{array}$(16) (Figure 13(c)). In summary, the sample selection rule of SARM is too relaxed to be extended to a larger extent. In contrast, the strict sample selection rules of RSG-OC are more universal to generating accurate rapeseed samples, and these types of vegetation will not be selected by the rules of RSG-OC due to the low DDYI. However, a strict rule of flowering features could generate a biased sample set that cannot represent the rapeseed pixels without clear flowering features. Therefore, another improvement of RSG-OC compared to SARM is the employment of a specifically designed sample augmentation that removes the observations during the flowering period. Such augmentation avoids the classifier from learning the biased selecting rule itself, thus enhancing the generalization performance of the trained classifier.

Figure 13. The WRI and DYI time series of (a) barley, (b) grass in Qinghai and (c) evergreen forest in Sichuan. The time series was selected from one typical parcel for each crop type.

5.3. Sensitivity of sample selection rule

In this study, the sample generation rules can generate accurate rapeseed samples with easy-determining thresholds. To illustrate the robustness to the varying threshold, the number of pixels satisfying the sample selection rules, the total area of mapped rapeseed (denoted as A_rapeseed), and the classification accuracy at different thresholds of the main rules (Eq. (9)(9) $\begin{array}{c}\begin{array}{c}\mathrm{D}\mathrm{D}\mathrm{Y}\mathrm{I}>0.04\end{array}\end{array}$(9) ). Considering that the OA and F1-score do not reflect the error in the mapped area due to the inconsistency between the proportion of rapeseed in the ground truth data and the proportion of the actual area, the proportion correct (PC) was also used as an indicator of accuracy (Pontius and Millones 2011), which was calculated using Equation 18(18) $\begin{array}{c}\begin{array}{c}PC=\sum _{i=1}^J{p}_{ii}\end{array}\end{array}$(18) to Equation 19(19) $\begin{array}{c}{p}_{ij}=\frac{n_{ij}}{{\sum }_{j=1}^J{n}_{ij}}\times \frac{N_i}{{\sum }_{i=1}^J{N}_i}\end{array}$(19) :

(18) $\begin{array}{c}\begin{array}{c}PC=\sum _{i=1}^J{p}_{ii}\end{array}\end{array}$(18)

(19) $\begin{array}{c}{p}_{ij}=\frac{n_{ij}}{{\sum }_{j=1}^J{n}_{ij}}\times \frac{N_i}{{\sum }_{i=1}^J{N}_i}\end{array}$(19)

where J is the number of classes, which is two here; $n_{ij}$ is the number of samples attributed to class i that are classified as class j; and $N_i$ is the total number of pixels classified as category i.

As shown in Figure 14, as the DDYI threshold becomes strict, the number of pixels satisfying the sample selection rules always decreases steadily. However, A_rapeseed decreases at different rates at different thresholds. When the threshold is relaxed, the A_rapeseed drops largely. As the threshold tends to be strict, the A_rapeseed gradually stabilizes. Similarly, as the threshold becomes stricter, the OA, F1-score, and PC generally become stable, although the PC is much lower when the threshold is relaxed. The cause can be explained as follows. When the threshold is relaxed, the selected samples are mixed with many inaccurate samples. Thus, more nonrapeseed was misclassified as rapeseed, leading to large classification error. As the threshold becomes stricter, the flowering features of selecting samples could be more biased for representing ordinary rapeseed pixels. Fortunately, the sample augmentation that removes the flowering features allows the classifier to focus more on nonflowering features. Therefore, bias in flowering features does not lead to large changes in A_rapeseed and accuracy. Therefore, it is easier to select a strict global threshold for large-extent rapeseed mapping.

Figure 14. The variation of (a) number of pixels satisfying the sample selection rules, (b) A_rapeseed, and (c) accuracy with DDYI threshold.

5.4. Cross-year generalization ability of RSG-OC

In this study, the PUL-RF trained by the samples generated from 2019–2021 was also used for mapping rapeseed in 2017–2018 (Figure 8). Unfortunately, there was a lack of ground truth samples in 2017–2018 for evaluation. Thus, to confirm the reliability of such a cross-year mapping strategy, we regenerate two subsample sets by excluding the generated samples in 2020 and 2021 and retrain two PUL-RF classifiers by these two regenerated sample sets (abbreviated as the non2020-model and non2021-model, hereafter). Then, the ground truth samples in 2020 and 2021 were used to evaluate the non2020-model and non2021-model, respectively. The performance was also compared with the standard PUL-RF classifier trained by all samples (abbreviated as all-model hereafter).

As shown in Figure 15, the average OA of the non2020-model and non2021-model only decreased by 0.14% and 0.32% compared to the all-model, respectively. Similarly, the average F1-score of the two models only decreased 0.08% and 0.28% compared to the all-model. Thus, classifiers trained by the generated samples from other years are able to map rapeseed in the years without generated samples, suggesting that the 2017–2018 rapeseed delineated by the classifier for 2019–2021 is also reliable.

Figure 15. Comparison of (a) average OA and (b) average F1-score before and after excluding the 2020 and 2021 training samples.

5.5. Influence of the predicting error of flowering date

In this study, the flowering features used for PUL-RF were extracted from the flowering period of [T-15, T + 15], where the flowering date T is predicted at each location by our previously proposed model (Zang et al. 2020). However, the predicted flowering date may be biased because the actual flowering date varies from year to year due to climate interannual variability. Therefore, we examined how the predicting error of flowering date could influence the classification accuracy. Firstly, we compared the 734 flowering dates recorded at the agrometeorological stations from 2003–2017 and their corresponding predicted flowering dates. The predicted RMSE of flowering date was 13.09 days (Figure 16(a)), and the predicting error of flowering date obeyed a normal distribution with a mean of 0.23 and a standard deviation of 13.09 (Figure 16(b)). It means that the predicting error is less than 13 days with a probability of 68.2% and is less than 26 days with a probability of 95.5%. Thus, we manually added −30 to 30-day offsets (one offset every 5 days) to the predicted flowering dates to simulate the possible errors of the predicted flowering date; and then evaluated the classification performance under each offset (Figure 16(c)). It is found that average OA and average F1-score stabilized above 94% and 93% respectively if the predicting error of flowering date is less than 15 days. If the predicted flowering date deviated 25 days, the average OA and average F1-score could slightly drop to 89% and 87% respectively. Thus, the proposed RSG-OC could tolerate the predicting error of flowering date to a large degree, which helps to relieve the influence of interannual fluctuations in the flowering time of rapeseed. However, if more accurate flowering dates were predicted by some advanced phenological models (Chuine and Régnière 2017) in the future, the mapping accuracy could be further stabilized.

Figure 16. (a) Relationship between predicted flowering date and actual flowering date. (b) Distribution of errors in predicting flowering date. (c) Effect of bias in predicted flowering date on classification accuracy.

5.6. Implications for future research

This study firstly produced the annual rapeseed maps at 20-m resolution for the whole of China, which exhibited detailed spatial patterns and temporal variation of rapeseed cultivation in China. This data could provide an important basis for agricultural subsidies, land use planning, and other policy formulation or evaluation (Zhang, Liu, and Zhang 2022; Zang et al. 2020). For example, a jump in rapeseed areas from 2019 to 2020 could reflect the effect of the policy reform in agriculture. In 2015, the Chinese government abolished the national rapeseed storage policy (Han 2015), resulting in a relatively low rapeseed price and less motivation for rapeseed planting for farmers in 2017–2019. Since 2019, several provinces with leading rapeseed productions started introducing subsidizing policies to promote rapeseed cultivation (e.g. General Office of Zhejiang Provincial Government 2020; Jiangsu Rural Statistics Division 2019; Anhui Provincial Department of Agriculture 2019), stimulating a rapid increase of rapeseed area from 2019 to 2020. The produced rapeseed maps provide evidence independent from census data for evaluating the policy effect. With much more spatial details than census data, the annual rapeseed maps derived from remote sensing could further support fine-grain analysis of the spatial pattern and the driving forces of the rapeseed area change. Moreover, the rapeseed maps are also the basic data layer for rapeseed growth monitoring with remote sensing. It thus provides opportunities for revealing more detailed spatial patterns of the rapeseed phenology, yield, and other agricultural information derived from remote sensing, compared to traditional plot-scale studies (D’Andrimont et al. 2020; Sulik and Long 2020; Mercier et al. 2020).

The proposed RSG-OC employed a promising strategy for effectively combining expert (or domain) knowledge and machine-learning techniques. Recently, deep learning has shown great superiority in improving crop classification accuracy (Xu et al. 2021; Yang et al. 2022; Ge et al. 2021). However, deep learning techniques often rely on a large number of training samples, which is challenging for large-scale crop mapping. On the other hand, expert knowledge for crop classification has been intensively explored for a long time in remote sensing community, especially for the widely-cultivated crop types like wheat, corn, peanut and rice (Dong et al. 2016; Qiu et al. 2017, 2018, 2021; Zhan, Zhu, and Li 2021). Effectively integrating the domain knowledge might help to reduce the relying on the large training samples for deep (Dash et al. 2022; Xie et al. 2021). RSG-OC converts the expert knowledge into automatically generated samples, which provides a straightforward strategy to integrate expert knowledge into machine learning techniques. Such a strategy does not need to modify the available machine learning techniques, and thus should not be difficult to be applied for deep learning and also for the other crop classification.

5.7. Remaining issues

There remain several issues in this study. Firstly, the ground truth data for accuracy evaluation were only obtained in the major rapeseed planting regions, thus is not adequate to represent the whole of China. Especially, the commission error in non-major rapeseed planting regions might be underestimated. For example, we found that some classified rapeseed areas (35.04N, 116.37E) in Shandong province might be actually garlic (Figure 17) based on our recent communication with local farmers. Therefore, more validations in more areas should be further conducted to better evaluate the accuracy of the produced maps. Secondly, the rapeseed area derived from the produced map is much less than the census data. Such inconsistency could be induced by the asymmetric commission/omission classification error, the effect of mixing pixels, and the uncertainty of census data (Boschetti, Flasse, and Brivio 2004; Czaplewski and Catts 1992; Gallego 2004; Waldner and Defourny 2017; Ozdogan and Woodcock 2006; Liu et al. 2020). Even with accurate land cover maps, there is still a possible area for the data with mixed pixels, especially in the case of fragmented landscape (Ozdogan and Woodcock 2006). The exact reasons for the inconsistency between mapped area and census data would be further explored in future research. Finally, the RSG-OC generated samples from the whole of China and trained only one national-wide classifier for rapeseed mapping, which might lead to underfitting for the regions where the generated samples are inadequate. Introducing a zoning strategy on sample generation and classifier training might further improve the performance of RSG-OC.

Figure 17. Example of misclassified rapeseed fields (35.04N, 116.37E) without yellowness in flowering season. (a) classification map, (b) true color image in flowering season.

6. Conclusion

Chinese rapeseed for 2017–2021 was mapped using the RSG-OC method in this paper. With specifically designed selection rules and augmentation methods for generating rapeseed samples, RSG-OC trained a classifier with less dependence on optical flowering features and is thus generalizable to pixels with serious cloud contamination during the flowering period. With the implementation of the PUL one-class classification strategy, the RSG-OC method avoids the difficulty of collecting diverse nonrapeseed samples across all of China. Evaluating the ground truth samples confirms that RSG-OC performs better in rapeseed mapping than the other methods. The rapeseed map produced for 2017–2021 filled the gap in long-term-span and large-extent rapeseed maps in China, which is of importance to national food security and land use regulation.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The Sentinel data (Sentinel-1, Sentinel-2 L2A, Sentinel-2 L1C) used in this paper are available in GEE; The produced rapeseed maps are also shared in GEE (https://code.earthengine.google.com/b457540737690943e0f0ac1d7dd41ed6) and Zenodo (https://doi.org/10.5281/zenodo.7047270). The processing code is available at (https://code.earthengine.google.com/7e7484b3ca4d9d1bd5896956c9c45ddc).

The validation data in this study are available from the corresponding author, X.H.Chen, upon reasonable request.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2022.2163576.

Additional information

Funding

This study is supported by the National Key Research and Development Program of China (No. 2022YFD2001101) and National Natural Science Foundation of China (No. 41871224).