2020 RSE《Mapping cropping intensity in China using time series Landsat & Sentinel-2 images and GEE》

Introduction

        阐述 cropland 和 croping intensity 的定义、croping intensity 数据的重要性。

Liu, J., Liu, M., Tian, H., Zhuang, D., Zhang, Z., Zhang, W., Tang, X., Deng, X.J.R., 2005. Spatial and Temporal Patterns of China’s Cropland during 1990–2000: An Analysis Based on Landsat TM Data. 98. pp. 442–456.

       已有研究：使用 MODIS（中等分辨率、高时间分辨率）识别作物密度。利用物候，通过一年中的植被指数峰值数量，对作物密度进行量化。

        MODIS 存在 sub-pixel heterogeneity。需要使用更高空间分辨率的遥感影像。阈值法要求每个生长季都至少有一张高质量影像，但 Landsat 时间分辨率较低，可用性较低。

        Sentinel-2 的两颗卫星的合并重放频率较高（5天），且其光谱信息和空间分辨率与 Landsat 相似。融合 Landsat 和 Sentinel 数据的意义：可以有效解决 Landsat 可用性较低 和 MODIS 空间分辨率较低的问题。

Materials and methods

2.1. Study area

2.2. Landsat and Sentinel-2 images and pre-processing

2.2.1. Image selection

        为了绘制2017年的作物种植密度图，考虑到作物密度包括冬季作物（冬小麦、冬油菜），选择了 2016年8月 ~ 2018年6月 的Landsat和Sentinel-2影像。使用 GEE提供的 TOA 数据（Sentinel MSI、Landsat TM、Landsat ETM+、Landsat OLI）。

TOA：top-of-atmosphere reflectance

ETM+：Enhanced Thematic Mapper Plus

OLI：Operational Land Imager

2.2.2. Image quality assessment

        低质量数据包括 cloud、cloud shadow、cirrus、snow、Landsat 7 ETM+ SLC-off gaps，大概占据一张影像 22% 的像元。使用 FMask 算法识别 cloud、cloud shadow、cirrus、snow。使用 ETM+ metadata 识别 SLC-off gaps。

SLC：scan line corrector

        统计了研究区内的每个像元在研究时期内的高质量观测数据。Sentinel MSI 的可用数据比Landsat ETM+ 和 OLI 之和 还多 2 倍。

        在整个研究期，研究区 A - G 中，拥有>200个高质量观测的像元比例分别为 87.11、89.97、61.53、44.10、44.66、 47.52、73.34 %。>100个高质量观测的像元比例分别为 87.82、 98.84、70.71、56.03、18.11、48.37、13.24 %。中国北部比南部的高质量数据多。

2.2.3. Imagery harmonization

       由于 ETM+、OLI、MSI 各波段的波长不同，必须进行统一化以获取可比较的结果，本实验以 OLI 数据为标准，对 ETM+ 和 MSI 进行调整(harmonize)。使用最小二乘回归系数，转换 ETM+ 的 波段 3 (红)、4 (近红外)、 5 (短波红外)（Roy等, 2016）和 MSI 的 波段 4 (红)、8A (近红外)、(短波红外)（Zhang 等, 2018）。

Roy, D.P., Kovalskyy, V., Zhang, H.K., Vermote, E.F., Yan, L., Kumar, S.S., Egorov, A., 2016. Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity. Remote Sens. Environ. 185, 57–70.

Zhang, H.K.K., Roy, D.P., Yan, L., Li, Z.B., Huang, H.Y., Vermote, E., Skakun, S., Roger, J.C., 2018. Characterization of sentinel-2A and Landsat-8 top of atmosphere, surface, and nadir BRDF adjusted reflectance and NDVI differences. Remote Sens. Environ. 215, 482–494.

2.2.4. Vegetation indices

        NDVI=(NIR=R)/(NIR+R)   (Tucker, 1979) 

        LSWI=(NIR-SWIR)/(NIR+SWIR)   (Xiao et al., 2005)

        注：red (630–680 nm)、near infrared (845–885 nm)、shortwave infrared (1560–1660 nm)

2.2.5. Preparation of 10-day composite images

        由于各个传感器的数据具有高度的 overlap 和 sidelap，许多地区可能会在同一天同时存在不同时间的数据，在同一天不同时间的植被指数会有微小区别。在一些天内，不同传感器可能会扫描同一个地区多次，从而造成植被指数的差异。

        本研究通过计算十天内所有观测数据的 NDVI 最大值，获取 10天 合成 NDVI；通过计算十天内所有观测数据的 LSWI 平均值，获取 10天合成 LSWI。(Running 等., 1995)

Running, S.W., Loveland, T.R., Pierce, L.L., Nemani, R., Hunt, E.R., 1995. A remote-sensing based vegetation classification logic for global land-cover analysis. Remote Sens. Environ. 51, 39–48.

2.2.6. Gap filling and smoothing

        由于云、雪等影响，有些区域有时可能没有高质量观测数据。可以基于 该 time step 之前或之后的高质量观测数据，进行时序线性插值，来填充这些 gap。(Kandasamy et al., 2013)

Kandasamy, S., Baret, F., Verger, A., Neveux, P., Weiss, M., 2013. A comparison of methods for smoothing and gap filling time series of remote sensing observations–application to MODIS LAI products. Biogeosciences 10, 4055–4071.

        为了去除残留噪声，需要进行平滑（smooth）。本研究使用 Savitzky–Golay filter 来重建 NDVI 数据集。具体参数：a moving window of 9 observations、a filter order of 2  (Fischer et al., 2002)。由于 LSWI 在干湿情况下变化较大，LSWI 的平滑没有必要。

Fischer, G., Van Velthuizen, H., Shah, M., Nachtergaele, F., 2002. Global Agro-Ecological Assessment for Agriculture in the 21st Century: Methodology and Results. International Institute for Applied Systems Analysis.

2.3. Algorithms for Identifying cropping intensity in a year for individual pixels

2.3.1. Phenology and signature analyses of cropping intensity

        一般特征：在 green-up 时期 NDVI 逐渐上升，成熟(mature)时 NDVI 达到顶峰，随后，从衰老期(senescence)一直到收获(harvest)，NDVI 逐渐下降。因此，可以从 NDVI 时序中提取作物生长周期的开始和结束时间，从而确定 cropping intensity。

        噪声：① 如 Fig. 6b 所示，冬小麦的NDVI在生长期并不是持续上升，因为冬小麦的NDVI在春化期(vernalization)会下降。② 如 Fig. 6d 所示，有些蔬菜的 NDVI 时序也存在抖动。③ 如 Fig. 6c 所示，有些水稻种植区，早稻收获和晚稻移栽会在2个礼拜之内进行。晚稻移栽后，LSWI 上升很快。如果这两个礼拜没有合适的观测数据，可能会观测不到裸土，从而无法利用裸土来判断作物生长周期的起始时间。④ 在作物种植之前或收获之后，田地里可能会生长杂草。

        对于①和②，利用裸土(bare soil)来判断生长周期的开始和结束时间。裸土的 LSWI 比绿色植被小很多。对于③，采用一个额外算法来识别这种 short rice harvesting–transplanting period。对于④，杂草的NDVI峰值一般小于作物，可以通过 NDVI 最值来判断作物的生长周期。

2.3.2. Algorithm for crop intensity mapping

Step 1：查找峰值（peak-finding）

        通过移动窗口，查找 NDVI 时序中的所有峰值。如果某个时间点的 NDVI 值大于该时间之前和之后的 NDVI，就被认为是一个峰值( peak )。同理，若某个时间点的 NDVI 值小于该时间之前和之后的 NDVI，就被认为是一个谷值( trough )。

Step 2：识别裸土，从而确定生长周期的起始时间

        计算整个低谷时期(trough period) LSWI 最大和最小值，并通过一个动态阈值方法识别裸土。动态阈值由下式定义，LSWI小于阈值()的时期，被认为是裸土。

       上式中，表示潜在的 LSWI 阈值，为最终阈值。和分别是两年内 LSWI 的最大值和最小值。

        研究表明，裸土的识别阈值：中国北部 LSWI<0 (Biradar 和 Xiao, 2011; Dong 等, 2015)，中国南部 LSWI<0 (Chen et al., 2018)。

Biradar, C.M., Xiao, X.M., 2011. Quantifying the area and spatial distribution of double- and triple-cropping croplands in India with multi-temporal MODIS imagery in 2005. Int. J. Remote Sens. 32, 367–386.

Dong, J., Xiao, X., Kou, W., Qin, Y., Zhang, G., Li, L., Jin, C., Zhou, Y., Wang, J., Biradar, C., 2015. Tracking the dynamics of paddy rice planting area in 1986–2010 through time series Landsat images and phenology-based algorithms. Remote Sens. Environ. 160, 99–113.

Chen, B., Xiao, X., Ye, H., Ma, J., Doughty, R., Li, X., Zhao, B., Wu, Z., Sun, R., Dong, J., 2018. Mapping forest and their spatial–temporal changes from 2007 to 2015 in tropical Hainan island by integrating ALOS/ALOS-2 L-band SAR and landsat optical images. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 11, 852–867.

Step 3：识别水稻

        若 0.2 < NDVI < 0.5，则该像元为裸土和植被混合体。若 NDVI > 0.5，则认为该像元全部是植被（Sobrino 等, 2001）。

Sobrino, J., Raissouni, N., Li, Z.-L.J.R.S., 2001. A Comparative Study of Land Surface Emissivity Retrieval from NOAA Data. 75. pp. 256–266.

       在 6 ~ 9 月之间，两个波的峰值 NDVI 都 > 0.5，且两个波之间的低谷 < 0.5，则这两个波都被认为是作物生长周期。

Step 4：若 NDVI 波峰 < 0.5，则该波不是作物生长周期。

Step 5：确定作物生长周期的起止时期（SOS 和 EOS）

       NDVIratio = (NDVI - NDVImin) / (NDVImax - NDVImin) 

        上式中，NDVImax 是作物生长周期内的NDVI最大值，NDVImin 是两年内的NDVI最小值。SOS 和 EOS 的阈值分别为 0.1 和 0.19。

Step 6：确定 MCI（multiple cropping index）

       如果某个生长季的 SOS 或 EOS 在2017，则该生长季的 MCI = 0.5。如果 SOS 和 EOS 都在 2017 年，则 MCI = 1。2017年的 MCI 为所有 MCI 之和。

2.3.3. Regional implementation of algorithms

        单季：1 ≤ MCI < 2。双季：2 ≤ MCI < 3。三季：MCI ≥ 3。

2.4. Accuracy assessment of cropping intensity maps

Results

3.1. Accuracy assessment of cropping intensity maps

3.2. Comparison of cropping intensity maps with national statistical data

3.3. Maps of cropping intensity in 2017

Discussion

4.1. Integration of times series Landsat and Sentinel-2 imagery

4.2. Algorithm development and critical conditions

4.3. Potential sources of uncertainty

4.4. Advantages of the Google Earth Engine

4.5. Implications and future work

Conclusion