2001—2015年海南岛橡胶林物候时空变化特征分析

胡盈盈， 戴声佩， 罗红霞， 李海亮，李茂芬， 郑倩， 禹萱， 李宁

(1.中国热带农业科学院科技信息研究所/海南省热带作物信息技术应用研究重点实验室，海口 571101；2.农业农村部农业遥感重点实验室，北京 100081； 3.中国热带农业科学院环境与植物保护研究所/海南儋州热带农业生态系统国家野外科学观测研究站，海口 571101)

0 引言

植被物候是植物长期适应季节性变化而形成的生长发育节律，是目前指示植物生长发育的季节特征和生理活动最常用的一个指标[1]。森林植被作为陆地生态系统的重要组成部分，在全球物质和能量循环中起着重要作用。开展森林植被物候研究不仅有助增进森林植被对气候变化响应的理解，对评估植被生产力和碳平衡也具有重要意义[2-3]。海南岛是重要农产品天然橡胶保护区划定的重点区域，天然橡胶是一种必不可缺的战略物资，监测海南岛橡胶林物候变化特征对于重要农产品保护以及生态保护具有重要意义。

长时间序列物候监测和评估成为目前植被物候研究的主要方向，而遥感作为一门对地观测综合学科，由于其范围广、时序长、可获取性强，最近几十年来被广泛应用在植被物候监测领域[4-6]。MODIS遥感影像由于其描述植被信息受到的干扰较小，时间序列较长和分辨率较高，许多学者利用其归一化植被指数(normalized difference vegetation index，NDVI)对植被物候变化进行研究[7-8]。目前，已有学者对橡胶物候进行了研究，Razak等[9]使用Landsat NDVI数据提取了马来西亚橡胶林物候参数，实现了橡胶林制图； Kou等[10]和Zhai等[11]分别利用Landsat数据和物候观测资料研究了西双版纳橡胶林物候特征； 陈汇林等[12]和陈小敏等[13]使用MODIS时间序列年内均值变化及邻域变化量来提取海南岛橡胶林物候参数； 田光辉等[14]基于2001—2011年海南岛橡胶林物候特征提取了橡胶林分布信息。上述研究或局限于橡胶物候的应用，或研究物候参数不全面，对橡胶林生长季长度鲜有涉及。

海南岛作为一个独立的生态系统，在全球气候变化背景下，橡胶林物候会出现什么样的规律特征，需要进一步探讨。因此，本研究基于2001—2015年MODIS NDVI数据，采用Savitzky-Golay(S-G)滤波法、动态阈值法提取海南岛橡胶林展叶期开始日、落叶期开始日和生长季长度，对15 a间橡胶林的物候变化进行动态监测，并基于橡胶林像元分析物候变化趋势和显著性，以期摸清海南岛橡胶林物候时空变化规律，为其他研究提供基础数据和参考。

1 研究区概况及数据源

1.1 研究区概况

海南岛(N18°10′～20°10′，E108°37′～111°30′)地处中国最南端，是海南省的主体，是我国第二大岛，全岛面积约3.54万km2。海南岛是我国橡胶林的主要种植基地之一，中西部和中部地区多为橡胶林种植区，如儋州市、澄迈县、白沙黎族自治县、琼中黎族苗族自治县等[15-16]。海南岛属热带季风气候，年均气温为22～26 ℃，年均降水量在1 600 mm以上，地势中间高四周低。综合地理位置、气候、胶林面积、胶龄、经营主体等多种因素，在海南岛西部、南部和中部分别选择3个典型样区。海南岛典型样区及数字高程模型(digital elevation model,DEM)如图1所示。

1.2 数据来源与预处理

遥感数据为2001—2015年MOD13Q1 NDVI数据，空间分辨率为250 m，时间分辨率为16 d，每年有23期影像数据，数据格式为EOS-HDF。使用MODIS Reprojection Tool工具对海南岛原始影像进行镶嵌、裁剪以及投影转换； 然后对NDVI数据去除无效值、归一化处理，将值控制在-1～1； 物候参数的提取在Timesat软件中完成，需将原始数据HDF格式转换成Timesat有效的DAT格式，归一化处理和格式转换均在ArcGIS软件中完成。

海南岛是天然橡胶保护区，橡胶为多年生热带经济作物，经济寿命长，种植至少在30 a以上[17-18]，受《森林法》保护，不可大面积砍伐。Dong等[19]研究表明基于Landsat/PALSAR橡胶林分类精度比MODIS/PLASAR高，橡胶林划定以Chen等[18]基于Landsat /PALSAR(25 m) 绘制的2010年海南岛橡胶林分布图为准，利用海南岛橡胶林分布图对NDVI时间序列数据进行掩模处理，以便精准提取橡胶林。

2 研究方法

2.1 物候提取方法

为了减少噪声的影响，利用S-G滤波对NDVI时间序列进行去噪平滑处理。S-G滤波[20]是利用滑动窗口对曲线进行分段的多项回归，并加权平均，拟合公式为：

(1)

式中：Y*为拟合后值；Yi+j为第i+j个像元的原始NDVI值；j为滑动窗口内纵向原始数据下标；Ci为第i个像元的滤波系数，由最小二乘法确定；n为窗口宽度(n=2m+1)，其中m为半个平滑窗口的宽度。采用S-G滤波对时间序列数据拟合时，移动窗口的宽度n和迭代次数d是应用中的关键，n不仅决定了平滑度，也影响了生长季节开始时跟踪曲线快速变化的能力[21]；d越高可以去除更多异常值，但是也会使结果噪声更多。本研究依据经验[17]和样区尝试，设定n和d都为3。

考虑到研究区橡胶林NDVI的变化特征，本研究采用动态阈值法提取橡胶林的物候参数。动态阈值法的提出者Jönsson等[21]建议植被生长季开始日期阈值为20%左右，但是对森林植被生长季开始与结束期提取时的阈值设定各有不同。根据已有经验并在多次试验的基础上，将生长季开始和生长季结束开始时间阈值分别设定为30%和60%。

2.2 物候趋势分析

为研究海南岛橡胶林每个像元点15 a间的物候变化趋势，基于最小二乘法的线性回归方法逐像元计算橡胶林物候变化率，并利用F检验分析其显著性，具体算法详见文献[22]，其公式为：

(2)

式中：Gslope为物候变化趋势;Gi为第i年物候开始(结束)日期，N为年变化量(N=15)。当Gslope>0，物候推迟； 当Gslope<0，物候提前。使用F检验进行显著性分析，当P<0.05时，趋势显著。

3 结果与分析

3.1 典型样区NDVI值拟合曲线特征

海南岛3个典型橡胶林样区Google Earth影像纹理如图2(a)—(c)所示。本研究在典型样区NDVI曲线重构时，前后各添加1 a的虚拟时间，图2(d)—(f)拟合曲线时序有69期，只选取中间时序23—46期进行分析。儋州市两院试验农场橡胶林展叶期开始日和落叶期开始日阈值分别设置为30%和50%。图2(d)可以看出，儋州市橡胶林在12月下旬—次年2月初明显落叶，2月初NDVI值最低，为落叶盛期； 在2月下旬开始展叶，6月下旬进入稳定期； 时序38—41间NDVI值出现下降，经分析该时段正处于8—9月雨季，云量较多，故NDVI值受影响下降。同理，保亭黎族苗族自治县热作所阈值设定为30%和60%。图2(e)所示保亭黎族苗族自治县橡胶林12月下旬开始落叶，2月为落叶盛期，3月下旬开始展叶，6月下旬—7月上旬NDVI值下降，可能是热作所橡胶林地块琐碎(图2(b))，像元混杂，或天气等因素造成的NDVI值下降。琼中黎族苗族自治县阳江农场阈值设定为30%和60%。图2(f)显示琼中黎族苗族自治县12月初橡胶林开始落叶，2月初—2月下旬达落叶盛期，3月中旬为橡胶林展叶期，之后NDVI值迅速上升，在5月下旬达到稳定期。根据典型样区试验，设定海南岛橡胶林展叶期开始日和落叶期开始日阈值分别为30%和60%。

(a) 儋州两院试验农场四队影像(b) 保亭热作所十一队影像(c) 琼中阳江农场老田队影像

(d) 儋州橡胶林NDVI时间序列曲线(e) 保亭橡胶林NDVI时间序列曲线(f) 琼中橡胶林NDVI时间序列曲线

图2 典型样区影像及拟合曲线

3.2 橡胶林物候结果验证

采用2017年儋州市(19°35′N,109°30′E)、白沙黎族自治县(19°22′N,109°28′E)、琼中黎族苗族自治县(19°14′N,109°48′E)橡胶林物候观测数据对研究结果进行验证，数据来源于海南儋州热带农业生态系统国家野外科学观测研究站。观测记录数据虽然比影像时间晚，但仍可用，因为观测站点橡胶林龄为15～20 a。本研究选取展叶期开始日、落叶期开始日2种观测数据，因生长季长度观测站点未记载，根据观测点记录的第一篷叶抽发期到落叶始期之间的时间计算得出。地面观测站点数据与相同站点15 a的平均物候数据对比，对比结果见表1。3个站点观测记录的展叶期开始日和遥感提取的结果相差不大。观测数据的落叶期开始日均大于遥感手段监测天数，相差6～11 d，原因是每年遥感影像最后一期(第23期)正好对应353 d，小于全年实际天数365 d，但2种结果均表现出橡胶林在年末开始落叶。橡胶生长季长度观测数据均大于遥感手段监测数据，由于人工观测和遥感监测的尺度差异较大，且遥感手段监测到的物候期和地面观测的物候期不具有传统意义上的一致性，区域内最早展叶的植物出现绿色，并不一定表明植物遥感生长季的开始，因为遥感监测的生长季开始日期是区域内多数植物开始展叶的日期，遥感监测的物候期相对地面单一植株的物候期有延迟，因此遥感监测的橡胶林生长季长度短于地面观测记录的生长季长度。2001—2014年海南省儋州市南丰镇那王村(19°24′20″N,109°33′44″E)站点观测的橡胶展叶期物候[13]平均在3月6日，本研究春季物候平均出现在3月4日，与14 a平均观测数据相差不大。本文研究结果与地面观测物候结果一致，证明了本研究结果的可靠性。

表1 不同物候测量方法的结果对比

3.3 橡胶林物候年际变化

从时间上看，海南岛橡胶林展叶期开始日呈提前趋势(图3)，从2001年的第70天(3月上旬)，波动提前到2015年的第42天(2月中旬)，约提前28 d，平均每年提前0.94 d(R2=0.16，P>0.05)。落叶期开始日整体呈延迟趋势，从2001年的第330天(11月下旬)延迟到2015年第356天(12月下旬)，约延迟了26 d，平均每年延迟0.84 d(R2=0.04，P>0.05)。生长季长度总体呈延长趋势，从2001年的第246 天(8个月)波动延长到2015年的第300 天(10个月)，大致延长了54 d，平均每年延长1.79 d(R2=0.10，P>0.05)。橡胶林展叶期开始日的提前趋势(0.94 d/a)，与1950—2015年全国木本植物春季物候变化(0.28 d/a)相比提前较多[23]，与1982—2006年东部南北样带亚热带落叶阔叶林返青期开始日变化(0.58 d/a)相比提前也较多[24]，比2001—2014年海南省儋州市南丰镇那王村观测点橡胶林展叶期变化(1.12 d/a)提前较少[13]，但在变化趋势上具有一致性； 秋季落叶期开始日的延迟趋势(0.84 d/a)，与1950—2015年全国木本植物秋季物候变化(0.19 d/a)相比延迟较多[23]，同1982—2006年东部南北样带亚热带落叶阔叶林休眠期起始日变化(0.80 d/a)相差不大[24]； 生长季长度的延长趋势(1.79 d/a)，较1950—2015年全国木本植被生长季长度变化(0.21 d/a)延长较多[23]，且大于1982—2006年东部南北样带亚热带落叶阔叶林生长季长度延长趋势(1.40 d/a)[24]。本研究春季展叶期开始日提前趋势是秋季落叶期开始日延迟趋势的1.1倍，推测生长季长度延长是由春季展叶期开始日提前引起的。

图3 2001—2015年海南岛橡胶林物候年际变化

3.4 橡胶林物候空间分布

海南岛15 a间展叶期开始日自中向西和东逐渐提前，落叶期开始日则相反(图4(a)—(b))，研究区海拔由中间向四周降低，表明物候的开始、结束日期与海拔有着密切关系。橡胶林平均展叶期开始日在第17—81天，即1月中旬—3月下旬(图4(a))，其中第17—81天开始展叶的橡胶林占全岛的81.38%，第49—81天占58.11%(图5(a))，即全岛一半左右的橡胶林在3月份开始展叶； 展叶期较早的地区集中在研究区西部和东部，均在第1—49天开始展叶。橡胶林平均落叶期开始日在第321—353天，11月中旬—12月下旬(图4(b))，第321—次年第1天开始落叶的橡胶林有71%(图5(b))，即绝大多数橡胶林在11月中旬—12月下旬落叶； 落叶期较早的地区集中在中部山区，落叶期较晚的地区有东方市、琼海市和万宁市等。橡胶林生长周期约持续209～305 d，达7～10个月左右(图4(c))，其中生长季长度225～273 d的橡胶林占全岛的54.04%(图5(c))，即海南岛大多数橡胶林生长季长度约在8～9个月； 琼海市和万宁市生长周期较长约在9～10个月左右。

(a) 展叶期开始日平均物候(b) 落叶期开始日平均物候(c) 生长季长度平均物候

(a) 展叶期开始日 (b) 落叶期开始日 (c) 生长季长度

3.5 橡胶林物候时空变化趋势

从空间上看，15 a来，海南岛橡胶林展叶期开始日总体上呈提前趋势(图6(a))，变化幅度为-2.3～-0.1 d/a，展叶期开始日显著变化地区较少(图6(d))，约占8.7%，变化显著区域中约81%的橡胶林显著提前，主要集中在白沙黎族自治县、屯昌市、琼中黎族苗族自治县、万宁市和琼海市等地。橡胶林落叶期开始日呈延迟趋势(图6(b))，延迟幅度集中在0.5～2.7 d/a，显著延迟地区占33.4%(图6(e))。海南岛橡胶林生长季长度呈增加趋势(图6(c))，增长趋势为0.05～0.9 d/a，生长季长度显著延长区域占13%(图6(f))，主要分布在儋州市和白沙黎族自治县等。对比橡胶林时空变化趋势，落叶期开始日显著延迟是橡胶林物候变化最显著的特征。

(a) 展叶期开始日年变化趋势(b) 落叶期开始日年变化趋势(c) 生长季长度年变化趋势

(d) 展叶期开始日通过(P<0.05)显著性水平检验(e) 落叶期开始日通过(P<0.05)显著性水平检验(f) 生长季长度通过(P<0.05)显著性水平检验

橡胶林春季物候的提前和秋季物候的延迟，这与近几十年来基于地面观测[13]和遥感监测[25-26]的研究结果一致。此外，本研究获取的橡胶林生长季长度是指从橡胶展叶期开始日到橡胶林落叶期开始日之间的时间累积，实际生产过程中，橡胶割胶时间是从篷叶发育完好后开始计算的，故本研究的橡胶林生长季长度较实际产胶时间偏长。就橡胶生长机理来看，该结果提取的生长季长度符合橡胶生长机理，可以准确地反映橡胶林的生长季长度。

橡胶林春季物候提前，秋季物候延迟，生长季长度延长，割胶期有可能提前和延长，产胶量可能会相应提高。但是，橡胶展叶期提前也会带来一些危害的提前发生，例如第一篷叶展叶期的提前橡胶树极大可能受到寒害，而且这个时期橡胶病虫害极易发生，尤其是白粉病，导致防治时间会提前和延长。因此，随着橡胶物候期的变化，还应加强橡胶树对极端气象灾害和病虫害的预防措施。

本研究中使用2010年橡胶林空间分布数据，并以此分析2001—2015年橡胶林物候时空变化特征。结合橡胶种植面积统计数据来看，2001—2010年和2010—2015年期间橡胶种植面积、种植区域均发生了变化。但是根据Chen等[27]研究表明，近几十年来新增种植橡胶林主要来源于老胶园更新，其次为耕地和林地转化，这表明橡胶林种植面积虽然有变化，但老胶园种植区的位置变化不大。从海南省各市县橡胶热带作物面积统计数据来看，2001—2010年海南岛南部、东南部内陆地区橡胶种植面积增长幅度较小，而北部橡胶种植面积增长幅度较大，主要有儋州市、澄迈县和临高县。因此，本研究结果可能由于影像时间和土地利用变化之间的差异存在误差，其中海南岛北部橡胶林物候结果差异较大，东部误差较小。

4 结论

本研究基于MODIS NDVI数据探讨了2001—2015年海南岛橡胶林物候的时空变化特征，基于S-G滤波和动态阈值法提取的物候结果与地面观测结果较一致，说明对海南岛橡胶林物候监测是可靠的，可为橡胶林生长监测和保护提供一定的参考依据。

1)本研究设定展叶期开始日和落叶期开始日阈值分别为30%和60%，可以准确地监测海南岛橡胶林物候特征。典型样区大约在12月开始落叶，2月进入落叶盛期，2月下旬—3月下旬陆续进入展叶期，典型样区物候提取方法可为局部地区物候监测提供有效的技术参考。

2)空间上，2001—2015年海南岛橡胶林从1月中旬—3月下旬开始展叶，11月中旬—12月下旬逐渐开始落叶，生长季长度长达7～10个月。同一时间段内，海南岛橡胶林物候空间分布与海拔空间分布具有较好的一致性。

3)时间上，15 a间海南岛橡胶林物候变化趋势不显著。展叶期开始日提前约28 d，落叶期开始日约推迟26 d，生长季长度延长约54 d，展叶期开始日提前是生长季延长的主要原因。

4)物候在空间尺度上变化趋势存在显著空间差异。15 a间橡胶林展叶期开始日变化趋势普遍提前，落叶期开始日变化趋势整体推迟，生长季长度变化趋势普遍延长，橡胶林落叶期开始日推迟成为橡胶林较显著的物候变化特征之一。