Year-to-year cyclical variations on global long time series of clumping index product
-
摘要:
植被聚集指数（clumping index，CI）是表征植被冠层叶片空间聚集程度的结构参数。基于作者所在研究团队研发生产的全球长时序（2001—2019年）、逐月的植被聚集指数(CI)遥感产品，使用傅里叶分解的方法逐像元地对CI的年际周期变化规律进行探究．步骤如下：1）预处理去除数据缺失较多的像元并填补部分缺失较少的像元，以产生相对完整的年际时间序列，预处理筛选出有效研究区占全球植被区面积的75.58%；2）假设CI年际变化可以分解为由多组余弦波信号与随机噪声构成的序列，使用离散傅里叶变换提取时间序列中振幅最大的余弦波（主波），发展了表征CI年际周期（1 a 1周期）变化参数指标体系；3）基于模拟数据和部分高质量数据检验了该方法的抗噪性，并将提取结果与MODIS物候产品（MCD12Q2）的1 a物候周期数（NumCycles，NC）及峰值时间（Peak_1）进行对比验证．结果表明：主波周期12（月）像元占比显著大于其他周期像元，占研究区的76.22%，表明CI最显著、最普遍的年周期性变化特征是周期长度12个月的年际变化；在研究区中主波周期等于12（月）区域与物候产品的1 a物候周期数（NC = 1）高度重合，精度达到96%；对主波周期12月的像元，主波低峰值月与物候产品peak值较接近，全球平均差异为1.37月，且CI低峰值月（植被最聚集月份）普遍提前于peak值，说明年周期植被叶片的最聚集状态（CI的季节变化低峰值）要普遍早于该周期植被叶片的绿度峰值（物候产品的peak最大值）．该研究为理解植被聚集效应的年际周期变化提供了证据．
-
关键词:
- 聚集指数 /
- 核驱动模型 /
- 二向性反射分布函数 /
- NDHD /
- MCD12Q2时间序列分析
Abstract:Vegetation clumping index (CI) is a structural parameter characterizing spatial clustering degree of vegetation canopy leaves, and playing an important role in extracting other vegetation parameters and in modeling land surface． Currently, remote sensing CI products are limited to a few multiangular satellite-borne sensors, compared to other remote sensing products (e.g., leaf area index, LAI). Therefore, study on long time series of CI products is needed: many models do not consider the cyclical variations of CI product in many applications． The global long-time series (2001-2019) and monthly CI products developed and produced by the author’s research team are therefore used here, to explore interannual periodic variations of 19-year CI product on a global scale, with Fourier decomposition technique． The CI product is first preprocessed by screening out and smoothening outliers, and filling in gaps, to generate a relatively complete year-to-year time series of CI product, with effective pixel proportion accounting for 75.58% of vegetation area on the global scale． Assuming that interannual variations of time-series CI product can be decomposed into a sequence composed of multiple sets of cosine wave signals with possible additive noise of random distribution, discrete Fourier transform is used to extract the cosine wave with the largest amplitude (main wave) in the time series． This is further used as the main index to characterize interannual periodic variation frequency． The period and amplitude of main waves are analyzed and further cross-compared with annual phenological NumCycles (NC) and peak time (Peak_1) of MODIS phenological products (MCD12Q2). Anti-noise performance of this method is tested by simulation data and high-quality CI pixels． The pixel proportions with = 12 are found significantly larger than in other periods, accounting for 76.22% of the study area． This indicates that interannual variations with 12-month period length is most significant among various periodic variations of CI time frequencies． Validation results show that the period with = 12 in the study area highly coincided with annual phenological period of MODIS phenological products with = 1, with validation accuracy reaching to 96%． For the pixel with = 12, low-peak month of the main wave is close to value of MODIS phenological product and average difference between them is about 1.37 months． In general, low-peak month of CI product (i.e., month with most clumping vegetation foliage) is earlier than peak value of MODIS phenological product, probably related to greenness of vegetation foliage． This implies that maximum clumping status is probably earlier to maximum greenness status for vegetation foliage． The present study provides evidence for improved understanding of interannual periodic variations for vegetation foliage clumping effect．
-
0 引言
植被聚集指数（clumping index，CI）是定量表征地表植被叶片聚集程度的参数，对提高叶面积指数（leaf area index, LAI）反演精度有重要作用[1-2]．Nilson于1971年通过马尔可夫链分析，首次对Beer定律的间隙率模型进行校正，定义了一个反映植被集聚程度的参数Ω，其表达式为[1]
$$ P(\theta ) = {{\mathrm{e}}}^{-G\left(\theta \right)I\varOmega /{\mathrm{cos}}\;\theta } , $$ (1) 式中：P(θ)为间隙率，即天顶角为θ的光线穿过冠层的概率；G(θ)是叶片投影函数，定义为单位叶面积在观测方向法平面上的投影面积，当叶倾角为球状分布时，其值为0.5；I表示LAI，其定义为单位地表上所有叶片的单面面积总和；Ω为聚集指数，LAI与Ω的乘积称为有效叶面积指数(effective LAI, LAIe)[2]，即，ILAe = Ω×I．自然生长植被冠层叶片常常呈聚集性分布，CI分布范围一般<1，且叶片分布越随机，CI越接近1．
现实场景中CI无法直接测量．地面获取CI一般采用异速生长法和间接光学测量得到．异速生长法基于目标参数和其他易测量的植被结构参数之间的经验关系来进行参数估算；间接光学测量则根据光学仪器测量冠层间隙率，并通过Beer-Lambert定律反演CI[3]．地基测量的数据相对精确，但难以大范围采集，一般仅用于验证；CI的遥感估算方法主要分为被动光学遥感估算法和激光雷达估算法[3]．大范围的遥感CI产品一般采用被动光学遥感方法估算[4-6]．多角度光学遥感提取CI的计算原理通常是基于CI与归一化热点和暗点差异指数（normalized difference between hotspot and darkspot，NDHD）的线性关系[4]．NDHD是Chen等[6]基于多角度观测数据提供的二向反射特征提出的一种角度植被指数，其表达式为
$$ I_{\mathrm{NDHD}} = \frac{{R}_{{\mathrm{Hot}}}-{R}_{{\mathrm{Dark}}}}{{R}_{{\mathrm{Hot}}}+{R}_{{\mathrm{Dark}}}} , $$ (2) 式中RHot 和RDark分别为热点和暗点方向的反射率．热点是当光照射方向与观测方向重合时观测到的反射率的一个峰值，此时可观测到的阴影面积最小[7]；暗点通常定义为前向散射方向处反射率的一个极小值，在该方向处观测到的阴影面积最大，暗点信号提供了聚集效应的大部分信息．Chen等[4]、He等[8]的研究表明，NDHD和聚集指数间存在着较好的线性关系．
CI与LAI等遥感参数密切相关，应用前景广泛[9−10]，但目前地表建模应用一般仅在每类地物类型中取单一CI值[11−13]．这主要是因为CI在各地物类之间的差异较为显著，相对其他植被参数，如LAI，CI随时间变化规律更为复杂[14]；目前，很少有针对全球CI随时间变化规律的研究，这对CI的实际应用产生了一定的限制，例如，一些研究在应用过程中常常不考虑CI随时间变化的特征[10]．
CI作为植被的一种结构参数，反演过程中除了随机噪声的影响，多角度NDHD方法中的热点反演的不确定性，使CI在反演过程中尤其在时间尺度存在较大不确定性，这些潜在的问题为CI随时间变化的研究带来很大挑战．目前，对CI时相变化特征的分析仍停留在一些简单的方法上，例如Wei等[14]提取全球 CI产品时涉及季节变化分析，把季节变化量化为1 a内月尺度CI数据的标准差除以整年均值，该方法较难以直接量化CI周期性变化特征，尤其是年际变化特征，虽然，文中地物类型的长时序CI散点图暗示CI可能存在周期性的季节变化，但并没有直接给出相关的定量分析和表征方法．
本文基于作者所在研究团队研发的长时序CI产品，假设CI长时序的变化是由一系列周期不同的正（余）弦波信号与随机噪声叠加而成，使用傅里叶分解分离出这一系列周期信号并提取振幅最大的波动作为主波，代表CI随时间变化的主要特征．主波的振幅大小代表CI随时间变化的显著性，且该振幅与CI值具有相同量纲，二者可直接加和运算；主波的周期代表CI主要变化的年际变化周期．当主波振幅显著且周期＝1 a时，说明CI变化主要有季节性年际变化特征；显著>1 a时，说明CI值可能存在某种长期变化；<1 a时，可能代表某些周期<1 a的物候现象，如一年两熟的作物等．
本文主要试图在全球范围定量地回答以下问题：
1）全球范围19 a的CI产品（2001—2019）是否存在主要变化周期为1 a的年际周期性变化，尤其在中纬度季相特征明显的区域；2）CI可能的年际周期性变化的振幅、相位等规律特征；3）通过与物候产品对比验证，研究CI的周期性变化特征与植被物候规律的可能关系．
1 数据和方法
1.1 数据
本文使用CI数据为作者所在研究团队研发生产的、基于MODIS/BRDF数据生成的、全球500 m分辨率、 2001—2019年逐月数据的CI卫星遥感产品[15]．该产品利用Jiao等人的MODIS热点反射率的校正模型，且对含雪像元进行异常值处理[4−5,16−17]，包括CI值与质量指标QA值（表1）；分辨率为500 m×500 m，全球不考虑海洋，分为315块区域，每块区域像元个数为
2400 ×2400 ，每块区域19 a逐月共288景数据．CI值范围为[0.33,1]，对有CI值的像元，QA分为0~3级，共4级，其中0、1表示CI是利用NDHD与CI线性关系的主算法估算的，CI的质量为0代表原始MODIS/BRDF产品的高质量（MODIS_BRDF_QA）为0或1；CI的质量为1代表原始MODIS/BRDF产品质量值为2或3；CI的质量为2、3是基于备用算法估算的[15]（表3）．特殊值共3类，即无植被覆盖区（32 765）、水域（32 766）和填充值（32 767）；CI读取值范围为3 300~10 000，尺度因子为0.000 1．主算法是基于热点改进版本的RossThick-LiSparseReciprocal（RTCLSR）模型，备用算法是当主算法计算值超出规定范围[0.33, 1]，使用基于各向异性平整指数AFX（anisotropic flat index）通过NDHD与CI建立的查找表进行估算[15−16]．表 1 聚集指数数据产品波段 取值范围 描述 CI [0.331]，特殊值 CI QA 0,1,2,3，特殊值 QA 表 3 聚集指数产品质量保证（QA）值描述QA值 描述 0 基于主算法计算，MODIS_BRDF_QA = 0~1 1 基于主算法计算，MODIS_BRDF_QA = 1~3 2 基于备用算法计算（无积雪覆盖），MODIS_BRDF_QA = 0~3 3 基于备用算法计算（有积雪覆盖），MODIS_BRDF_QA = 0~3 >30 000 无有效CI数据（填充值） 地表物候数据来源于美国地质勘探局，由波士顿大学研发的MODIS地表物候（MCD12Q2）产品，可通过https://lpdaac.usgs.gov/products/mcd12q2v061/ 网站获取 [18]，空间分辨率500 m，每年12幅影像．产品基于矫正的MODIS观测数据，利用每年8 d增强植被指数EVI（enhanced vegetation index）生成，该产品被广泛用来研究全球地面的季节性物候特征和年际变化[19]．
MCD12Q2产品主要用于验证和分析CI年际周期变化，其第1层数据（波段1）NumCycles（NC），代表周期数（Nc）．第2层数据（波段2）Peak_1，代表1 a内EVI曲线达到第1个周期峰值的日期（表2），时间为2020年．1 a内周期数指的是处理后的EVI曲线1 a内达到峰值的次数，Peak_1波段的数值是距1970年1月1日以来的时间，将其转化为月份后得到peak．有研究表明，全球植被的生长开始日期（SOS）与结束日期（EOS）的年际变化幅度较小，幅度约0.1~0.2 d·a−1[20]，而EVI峰值介于二者之间，这种差别对于本文所采用的月尺度CI研究数据而言其变化可以忽略，因此，我们在研究中只采用2020年的物候数据．
表 2 MCD12Q2数据产品波段 取值范围 描述 NumCycles（NC） 0,1,2,3,$ \cdots$,7 1 a内周期数 Peak_1（peak）* [ 11138 ,32766 ],整数1 a内达到第1个周期
峰值的时间注：* Peak_1单位为d， peak为将其计算为月份后的值（范围0~12）. 1.2 数据预处理
MODIS数据计算的CI产品中存在缺失值及可能的随机噪声．预处理包括：1）将19 a数据及质量指标升尺度到1 500 m；2)逐像元计算数据缺失率及异质数据率并筛选出有效像元；3)对有效像元进行插值处理使数据连贯．
对数据进行升尺度可以一定程度上平滑随机噪声，同时降低运算量；升尺度过程中以9个邻域像元的标准差（standard deviation，std）值为阈值筛除异质性较高的像元，减少像元异质性的干扰．阈值根据直方图法选择．以植被占比较多的h10v04区为例，图1展示了2019年不同月份MODIS h10v04区的std直方图，各月份直方图均相似．大部分像元std为0~0.3，极少量>0.5．最终阈值选择为0.3，以剔除这些异常值（图1）．
数据缺失值处理是预处理的关键部分：升尺度后含有缺失值的粗像元（
1500 m×1500 m）标记为缺失值，每个像元时间序列上缺失值>10%则被删除．保留下来的像元使用常用的3次样条插值补全时间序列（图2）．![]()
图 2 以一个3×3邻域像元为例显示数据预处理前后的变化a．升尺度前折线；b．升尺度后折线图及插值点；c．插值曲线及均值滤波曲线．该像元位于h10v04区，中心像元经纬度坐标44.3°N、109.7°W．升尺度过程中，在空间维度上使用取平均数并向上取整的方法，时间维度上对QA均采用取众数的方法得到升尺度后像元的QA．
对MCD12Q2产品的2个波段同样进行升尺度处理，并将Peak_1数值转化为月份，分别得到NC与peak（表2）．
1.3 方法原理
本研究主要采用傅里叶变换的方法逐像元处理时间序列，提取振幅最大余弦波（主波），记录主波的参数(振幅、周期、相位、均值)，每个参数都对应1个数据矩阵，根据这些参数，推导CI主波低峰值月，计算主波的振幅占比（表4），基于这些参数和指标，进一步分析CI年际周期变化属性．其中，主波表达式为
表 4 主波参数及衍生指标参数 范围 描述 主波参数 最大振幅$ A_{\mathrm{m}} $* (0,1] 主波振幅 (详见1.3) 最大振幅周期Tm（月） [2, 76] 主波周期（详见1.3） 初相位p [−inf, inf] 主波初相位(详见1.3节) 均值Icmean* [0.33,1] 对预处理后的CI序列取均值(详见1.3节) 指标 振幅占比Ap (0,1] 2倍$ A_{\mathrm{m}} $除以Icmean（式（1）~式（3））** CI低峰值月Mmin（月） 0,1,2,$ \cdots$,12 根据p及式（1）、（2）计算的主波达到低峰值的月份 注：*Am与 CI的量纲含义一致，将其直接与CI值进行加减运算是有意义的．**主波最大最小值之差为2倍$ {A}_{\rm{m}} $．该参数直接衡量了主波振幅相对均值带来的变化幅度． $$ {I}_{{\mathrm{cmain}}}\left(t\right) = A_{\mathrm{m}}\times \mathrm{cos}\left(\frac{2{\text{π}} }{T_{\mathrm{m}}}(t-p)\right)+{I_{{\mathrm{cmean}}}} \text{，} $$ (3) 式中：$ {I}_{{\mathrm{cmain}}} $为CI主波；$ A_{\mathrm{m}} $为最大振幅周期（与CI具有相同量纲）；t为时间（月）；$ p $为初相位（月），$ T_{\mathrm{m}} $为函数周期，当$ T_{\rm m} = 12 $时，可以确定CI峰值所在月份．$ I_{{\mathrm{cmean}}} $为CI均值．$ A_{\mathrm{m}}、p、T_ {\mathrm{m}}$可以通过FFT方法以及后续处理得到，$I_{{\mathrm{cmean}}} $则通过直接对预处理的时间序列产品取均值获得．
通过初相$ p $可计算出主波低峰值月．当Tm = 12时，以t为自变量，对于所有t满足t−p = 6+12i时，该式取得最小值．计算式为
$$ {M}_{{\mathrm{min}}} = {\mathrm{mod}}(6+p,12) \text{，} $$ (4) 式中：$ {M}_{{\mathrm{min}}} $为CI主波低峰值月；$ {\mathrm{mod}} $为取余函数；$ p $为式中的初相位．CI的低值指示了年周期变化中的植被叶片的最聚集状态，可能与植被生长期叶片的茂密程度有关．
振幅占比$ A_{\mathrm{p}} $定义为2倍$ A_{\mathrm{m}} $与均值之比，因为余弦波变化范围为2倍振幅，代表主波相对均值上下偏移的比例．
$$ A_{\mathrm{p}} = \frac{2A_{\mathrm{m}}(T_{\rm m})}{I_{{\mathrm{cmean}}}} \text{，} $$ (5) 式中：$ A_{\mathrm{m}}\left(T_{\mathrm{m}}\right) $表示最大振幅周期为Tm的主波振幅；Icmean使用的是序列插值后计算的均值；$ A_{\rm p} $可以用来衡量主波在CI均值附近造成的变化幅度．但需要注意，CI时间序列的极差值一般显著大于2倍$ A_{\mathrm{m}} $（图3、4），这是由于时间序列中可能叠加了噪声及振幅小于主波的周期性波动．
![]()
图 3 傅里叶变换分解提取主波注：模拟信号由2个三角函数组分与正态分布噪声及常数组成（如分图a、b、c、d），三角函数振幅分别为1.0与0.6，周期分别为12与6．噪声为服从正态分布的随机噪声（均值0，标准差1.2）；分图f中展示了提取出的主波；傅里叶分解展示了较好的抗噪性能，成功提取出了信号组成函数的周期（12）．噪声对振幅（Am）、相位、均值的提取影响较小．![]()
图 4 方法抗噪性检验注：左图为对典型像元添加随机噪声前后提取的主波差异展示；右图为向典型像元添加不同倍率的预期值为0的正态分布随机噪声，各参数相对没有噪声时的变化情况（添加噪声/无噪声）．噪声倍率为随机噪声标准差与原始Am的比值．噪声倍率较小时，各参数波动幅度均在0.9~1.1，受到影响较小．添加随机噪声后，2倍Am将明显小于极差．数据处理具体采用的方法为离散傅里叶分析变换(discrete Fourier transform，DFT)，通过DFT方法可以将时间序列分解为若干振幅、周期不同的正弦波组分的组合．编程应用中可以采用快速傅里叶变换 (fast Fourier transform，FFT) [21]实现DFT．当输入信号长度远大于可能的波动周期时，DFT方法可很好地判别主波周期，因此19 a长周期时间序列的数据有利于更精确地判断CI年际周期性变化特征．
1.4 检验方法
1.4.1 方法抗噪性检验
分别使用模拟信号及高质量数据对傅里叶变换方法进行了抗噪性测试．模拟信号由一系列振幅不同、周期为12的主波组分，1个振幅较小、周期为6的次波组分，以及1个均值为0、标准差略大于主波振幅的正态分布随机噪声相加．高质量数据假设为若干波形周期性较为明显的像元，以44.3°N、109.7°W像元为例，该典型像元波形呈现明显的周期性变化（图4），添加标准差为其主波振幅0~1.4倍的正态分布噪声（主波振幅为0.17），对比不同噪声程度下主波各参数相对不添加噪声时的变化百分比（添加噪声与无噪声之比）．其余高质量像元处理与该像元相同，最终结果较为相似．
1.4.2 计算结果对比验证
Tm = 12代表主波周期为1 a，与NC = 1的含义一致．最终本文在基于CI产品筛选出的有效范围内与存在MCD12Q2产品值的区域内对二者进行比较．流程如图5所示，比较Tm = 12像元与NC = 1的像元之间的重合率，较高的重合率指示CI存在1 a周期的季节变化，并在此基础上进一步对比全球Mmin与peak的差异．
Mmin与peak为取值在0~12中循环的参数，其值本身不代表大小，在长度12个月的周期循环中，二者可能的差异最大为提前或滞后6月．采用取模后取较小值方式计算这种差异：
$$ [{\mathrm{Dis,od}}] = \mathrm{m}\mathrm{i}\mathrm{n}\left({\mathrm{mod}}\right(P-M,12),{\mathrm{mod}}(M-P,12\left) \right) ，$$ (6) 式中：M为$ M_{\mathrm{min}} $；P为$ {\mathrm{peak}} $；Dis为计算出的P与M的差异； od为较小值序号，例如，$ {\mathrm{od}} = 1 $代表$ {\mathrm{mod}}(P-M,12) $较小，即$ {\mathrm{Dis}} = {\mathrm{mod}}(P-M,12) $，此时P较为滞后，od = 2则相反．
2 计算结果与分析
2.1 抗噪性检验分析
使用模拟信号对提取主波方法进行了测试．图3以主波组分振幅为1.0，次波组分振幅为0.6，正态分布噪声标准差1.2为例，显示傅里叶分解很好地提取出该序列的主波周期．示例中，主波各参数受到随机噪声的影响均较小．对噪声最不敏感的变量为周期．
图4显示了对坐标44.3°N、109.7°W的像元添加标准差为主波0~1.4倍（下称噪声倍率）的正态分布随机噪声后，主波参数的变化．当噪声倍率<0.7时，各参数变化幅度很小（90%~110%）．随着噪声倍率增加，受到随机噪声影响最小的是$ I_{{\mathrm{cmean}}} $，这是由于添加随机噪声的期望值为0．其次是$ T_{\mathrm{m}} $．噪声倍率较大时，$ M_{\mathrm{min}} $随机波动幅度最大，$ T_{\mathrm{m}} $则呈现递减趋势．其余高质量像元的检验结果较为相似．值得注意的是，主波变化范围（2倍$ A_{\mathrm{m}} $）一般明显小于原数据极差．综上所述，本方法对于假设CI长时序变化可能出现的随机噪声展示出了较好的抗噪性．
2.2 计算结果分析
经过对数据的预处理，以2019年6月CI产品标识的植被区为准，筛选出了占全球植被区面积75.58%的有效像元（图6、表5）．
表 5 各类像元面积占植被区面积比% 纬度 像元占比 实验区 Tm = 12 $ {A}_{{\mathrm{p}}} $>0.10 90~60°N 18.18 11.25 3.02 <60~30°N 90.02 76.61 35.69 <30~0°N 47.36 34.26 13.10 0~30°S 97.78 62.08 17.68 >30~60°S 94.08 79.25 44.03 总计 75.58 57.59 23.84 注：60~90°S 陆地为南极洲，没有植被覆盖，未参与统计．Tm是主波周期，Ap为振幅占比（表4）．0.1为所有 Tm = 12像元的Ap均值．使用的参考植被区选取为CI产品2019年6月全球非填充值区（表3）． 图7-a为计算出的全球Tm分布图，Tm = 12的像元占研究区面积76.22%，远大于其余Tm类别，表明研究区植被CI主要存在以12个月为周期的年际变化．图8展示了NC = 1与Tm = 12区域重合情况．表6列举了NC = 1与Tm = 12的重合率：研究区中70.90%的数据存在对应的MCD12Q2产品值，NC = 1占研究区面积67.24%；Tm = 12且有NC值的区域占研究区面积55.86%，其中Tm = 12与NC = 1的重合区占研究区面积54.04%，相对Tm = 12的重合率为(54.04/55.86) = 96.74%，相对NC = 1的重合率为(54.04/67.24) = 80.37%，证明植被1 a为1周期的物候规律同样在植被CI中显著呈现．
表 6 周期12月像元验证占比% NC Tm 总计 12 其他值 占实验区 相对* 1 54.04 13.20 67.24 80.37 其他值 1.82 1.84 3.66 49.37 占实验区 55.86 15.04 70.90 相对* 96.74 85.71 注：总计栏为MCD12Q2产品在实验区有值像元占实验区的比．相对栏（列）为第1栏（列）除以总计栏（列）的值（参数解释见表4）．在相对*栏，注意到NC = 1与Tm = 12重合区的相对占比分别为96.74%与80.37%． 图7-d为Tm = 12的像元的主波相位推算CI低峰值月Mmin．Mmin分布呈现一定的纬度地带性，Mmin在北半球多为5—7月，南半球主要为1、11—12月，均为当地夏季（基本为植被生长旺季），表明CI表征的叶片最聚集状态总体上出现在植被生长最旺盛的季节．图9对Tm = 12与NC = 1的重合区计算了Mmin与物候产品peak的差异，全球二者平均时差约1.37月，二者均有值区域76%为CI比物候产品的peak值（即EVI所代表的绿度极大值）提前达到低峰值（图9）．
图7-c所示为全球$ A_{\mathrm{p}} $值分布，美国东西部、西欧、南非等地$ A_{\mathrm{p}} $值较高；图10-b为$ A_{\mathrm{p}} $直方图，为单峰右偏分布，代表其均值小于概率分布曲线的峰值，且大部分数据大于峰值．计算不同Tm数据的$ A_{\mathrm{p}} $直方图峰值并画折线(图10-c)发现，当Tm = 12时$ A_{\mathrm{p}} $直方图峰值最大（7.18%），说明CI变化幅度最大的周期性变化是以1年为周期的年际变化．相比之下不同Tm = 12的CI均值直方图有多峰且分布较为相似（图10-a）．表5根据纬度带计算了各类像元占比，发现温带区域Tm = 12且$ A_{\mathrm{p}}$大于均值的像元比例最多．
3 结论与讨论
3.1 结论
本文使用作者所在团队研发的MODIS全球19 a（2001—2019年）长时序植被CI产品，对全球CI时间序列进行傅里叶变换分析，通过提取主波的主要参数及其衍生指标，结合MCD12Q2产品对结果进行验证，主要结论如下．
1）本文首先对像元尺度的CI数据19 a时间序列产品进行预处理，筛选出研究区域面积占全球植被区为75.58%．通过对傅里叶方法进行抗噪声检验，表明傅里叶变换方法对随机噪声有较强的抗噪性，对周期的识别受干扰小，可以很好地识别主波周期，主波各参数相对无噪声时变化幅度也很小．因此，使用傅里叶变换提取CI主波是简便有效的方法．
2）结合CI生物物理参数，为CI长时序的傅里叶变换参数赋予物理意义，并发展了相应的衍生指标．分析表明，年际周期变化是CI最显著的周期性变化特征．所有Tm类型中，Tm = 12（月）的像元占76%，主要分布于温带，且该类型$ A_{\mathrm{p}} $值显著大于其他Tm类型；Tm = 12的像元与NC = 1的像元高度重合，相对Tm = 12区域，重合率96.74%；相对NC = 1区域，重合率80.37%；表明1 a周期的变化同时体现在植被叶片的聚集状态变化（CI）和植被叶片的绿度状态变化（EVI）中．
3）对Tm = 12的像元，其CI主波达到低峰值的月份分布呈现出纬度地带性，时间以当地夏季为主，北半球CI低值出现月份多在5—7月，南半球则多在11—次年1月．对比MCD12Q2产品EVI波峰曲线时间，计算出CI达到低峰值与EVI达到峰值时间全球平均差异约1.37月（全球平均），主要为CI低峰值提前．这与相关研究中LAI达到峰值时间比冠层叶绿素含量达到峰值时间提前[22]的结论吻合．可能的原因是植被生长期间，植被叶片率先生长达到较聚集状态，而叶片中代表绿度的叶绿素浓度则存在逐渐增大的过程[22−23]．同时该结果也表明，当植被达到生长旺季，CI值会降低，这符合CI的物理意义，与前人研究结果吻合[14]．
3.2 讨论
在研究过程中发现，仍存在一些问题需要后续研究进一步开展讨论．
1）在不同Tm值的$ A_{\mathrm{p}} $峰值图中，Tm = 6的像元$ A_{\mathrm{p}} $值同样存在一个小高峰（可能暗指一年两熟的农作物），但对于Tm = 6与NC = 2的区域，并没有类似Tm = 12与NC = 1那样的高重合率，故本文未展示其结果．其原因可能为：① 本文使用数据时间分辨率为月，当提取的目标周期较小时，一个完整周期内有效数值较少，因此，易受到数据质量（如缺失值）的干扰；②预处理过程中对时间序列进行邻域均值滤波及插值处理，该处理虽然有利于CI年际变化的提取，但对年内短周期变化可能有一定平滑作用．
2）本文提出的参数Am可以与CI均值直接加减，具备一定的实际意义，可以很好地回答CI周期变化是否显著的问题，但该指标是基于长时间序列计算的一个总体趋势，默认所有年份振幅都是一个值，因此并不能体现特殊年份中较大或较小的变化幅度．它仅是衡量CI总体周期变化幅度的一个总体估计，指示大多数年份下，CI的周期变化幅度约为2倍Am．
3）本研究虽然表明CI有显著的年际周期性变化，但如何在应用中考虑CI值，例如，是需要代入每月的值，还是代入季度平均值，以及考虑时间变化后不同地物类的CI值如何选取等，需要进一步研究，因为代入每月值同时也意味代入CI噪声的不确定性，不同地物类的变化周期及变化幅度也可能存在显著差异．研究检测出了1 a周期的CI振幅占比一般为5%~15%，这仅是对长时序CI变化范围的总体估计．特殊年份CI值可能围绕均值产生比该范围更大的差异．不同应用情景下这种差异造成的影响仍需进行评估．
-
![]()
图 2 以一个3×3邻域像元为例显示数据预处理前后的变化
a．升尺度前折线；b．升尺度后折线图及插值点；c．插值曲线及均值滤波曲线．该像元位于h10v04区，中心像元经纬度坐标44.3°N、109.7°W．
![]()
图 3 傅里叶变换分解提取主波
注：模拟信号由2个三角函数组分与正态分布噪声及常数组成（如分图a、b、c、d），三角函数振幅分别为1.0与0.6，周期分别为12与6．噪声为服从正态分布的随机噪声（均值0，标准差1.2）；分图f中展示了提取出的主波；傅里叶分解展示了较好的抗噪性能，成功提取出了信号组成函数的周期（12）．噪声对振幅（Am）、相位、均值的提取影响较小．
![]()
图 4 方法抗噪性检验
注：左图为对典型像元添加随机噪声前后提取的主波差异展示；右图为向典型像元添加不同倍率的预期值为0的正态分布随机噪声，各参数相对没有噪声时的变化情况（添加噪声/无噪声）．噪声倍率为随机噪声标准差与原始Am的比值．噪声倍率较小时，各参数波动幅度均在0.9~1.1，受到影响较小．添加随机噪声后，2倍Am将明显小于极差．
表 3 聚集指数产品质量保证（QA）值描述
QA值 描述 0 基于主算法计算，MODIS_BRDF_QA = 0~1 1 基于主算法计算，MODIS_BRDF_QA = 1~3 2 基于备用算法计算（无积雪覆盖），MODIS_BRDF_QA = 0~3 3 基于备用算法计算（有积雪覆盖），MODIS_BRDF_QA = 0~3 >30 000 无有效CI数据（填充值） 
下载: 导出CSV
表 2 MCD12Q2数据产品
波段 取值范围 描述 NumCycles（NC） 0,1,2,3,$ \cdots$,7 1 a内周期数 Peak_1（peak）* [ 11138 ,32766 ],整数1 a内达到第1个周期
峰值的时间注：* Peak_1单位为d， peak为将其计算为月份后的值（范围0~12）.
下载: 导出CSV
表 4 主波参数及衍生指标
参数 范围 描述 主波参数 最大振幅$ A_{\mathrm{m}} $* (0,1] 主波振幅 (详见1.3) 最大振幅周期Tm（月） [2, 76] 主波周期（详见1.3） 初相位p [−inf, inf] 主波初相位(详见1.3节) 均值Icmean* [0.33,1] 对预处理后的CI序列取均值(详见1.3节) 指标 振幅占比Ap (0,1] 2倍$ A_{\mathrm{m}} $除以Icmean（式（1）~式（3））** CI低峰值月Mmin（月） 0,1,2,$ \cdots$,12 根据p及式（1）、（2）计算的主波达到低峰值的月份 注：*Am与 CI的量纲含义一致，将其直接与CI值进行加减运算是有意义的．**主波最大最小值之差为2倍$ {A}_{\rm{m}} $．该参数直接衡量了主波振幅相对均值带来的变化幅度．
下载: 导出CSV
表 5 各类像元面积占植被区面积比
% 纬度 像元占比 实验区 Tm = 12 $ {A}_{{\mathrm{p}}} $>0.10 90~60°N 18.18 11.25 3.02 <60~30°N 90.02 76.61 35.69 <30~0°N 47.36 34.26 13.10 0~30°S 97.78 62.08 17.68 >30~60°S 94.08 79.25 44.03 总计 75.58 57.59 23.84 注：60~90°S 陆地为南极洲，没有植被覆盖，未参与统计．Tm是主波周期，Ap为振幅占比（表4）．0.1为所有 Tm = 12像元的Ap均值．使用的参考植被区选取为CI产品2019年6月全球非填充值区（表3）．
下载: 导出CSV
-
[1] NILSON T. A theoretical analysis of the frequency of gaps in plant stands[J]. Agricultural Meteorology,1971,8:25 doi: 10.1016/0002-1571(71)90092-6
[2] 阎广建,胡容海,罗京辉,等. 叶面积指数间接测量方法[J]. 遥感学报,2016,20(5):958 [3] 方红亮. 真实和有效叶面积指数及聚集指数的尺度效应[J]. 地球信息科学学报,2021,23(7):1155 doi: 10.12082/dqxxkx.2021.200609 [4] CHEN J M,MENGES C H,LEBLANC S G. Global mapping of foliage clumping index using multi-angular satellite data[J]. Remote Sensing of Environment,2005,97(4):447 doi: 10.1016/j.rse.2005.05.003
[5] DONG Y D,JIAO Z T,YIN S Y,et al. Influence of snow on the magnitude and seasonal variation of the clumping index retrieved from MODIS BRDF products[J]. Remote Sensing,2018,10(8):1194 doi: 10.3390/rs10081194
[6] CHEN J M,LIU J,LEBLANC S G,et al. Multi-angular optical remote sensing for assessing vegetation structure and carbon absorption[J]. Remote Sensing of Environment,2003,84(4):516 doi: 10.1016/S0034-4257(02)00150-5
[7] HAPKE B,DIMUCCI D,NELSON R,et al. The cause of the hot spot in vegetation canopies and soils:shadow-hiding versus coherent backscatter[J]. Remote Sensing of Environment,1996,58(1):63 doi: 10.1016/0034-4257(95)00257-X
[8] HE L M,CHEN J M,PISEK J,et al. Global clumping index map derived from the MODIS BRDF product[J]. Remote Sensing of Environment,2012,119:118 doi: 10.1016/j.rse.2011.12.008
[9] FANG H L,LI S J,ZHANG Y H,et al. The relationship between canopy clumping index (CI),fractional vegetation cover (FVC),and LEAF area index (LAI):an analysis of global satellite products[C]//IGARSS 2020 - 2020 IEEE International Geoscience and Remote Sensing Symposium. Waikoloa,HI,USA:IEEE,2020:4120
[10] FANG H L. Canopy clumping index (CI):a review of methods,characteristics,and applications[J]. Agricultural and Forest Meteorology,2021,303:108374 doi: 10.1016/j.agrformet.2021.108374
[11] HE M Z,JU W M,ZHOU Y L,et al. Development of a two-leaf light use efficiency model for improving the calculation of terrestrial gross primary productivity[J]. Agricultural and Forest Meteorology,2013,173:28 doi: 10.1016/j.agrformet.2013.01.003
[12] ZHOU Y L,WU X C,JU W M,et al. Global parameterization and validation of a two-leaf light use efficiency model for predicting gross primary production across FLUXNET sites[J]. Journal of Geophysical Research:Biogeosciences,2016,121(4):1045
[13] BI W J,HE W,ZHOU Y L,et al. A global 0.05° dataset for gross primary production of sunlit and shaded vegetation canopies from 1992 to 2020[J]. Scientific Data,2022,9(1):213 doi: 10.1038/s41597-022-01309-2
[14] WEI S S,FANG H L,SCHAAF C B,et al. Global 500 m clumping index product derived from MODIS BRDF data (2001–2017)[J]. Remote Sensing of Environment,2019,232:111296 doi: 10.1016/j.rse.2019.111296
[15] JIAO Z T,DONG Y D,SCHAAF C B,et al. An algorithm for the retrieval of the clumping index (CI) from the MODIS BRDF product using an adjusted version of the kernel-driven BRDF model[J]. Remote Sensing of Environment,2018,209:594 doi: 10.1016/j.rse.2018.02.041
[16] JIAO Z T,HILL M J,SCHAAF C B,et al. An anisotropic flat index (AFX) to derive BRDF archetypes from MODIS[J]. Remote Sensing of Environment,2014,141:168 doi: 10.1016/j.rse.2013.10.017
[17] JIAO Z T,SCHAAF C B,DONG Y D,et al. A method for improving hotspot directional signatures in BRDF models used for MODIS[J]. Remote Sensing of Environment,2016,186:135 doi: 10.1016/j.rse.2016.08.007
[18] FRIEDL M,GRAY J,Sulla-Menashe D. MODIS/Terra+ Aqua Land Cover Dynamics Yearly L3 Global 500 m SIN Grid V061[Z]. Missoula,MT,USA:NASA EOSDIS Land Processes DAAC,2022
[19] XIAO W W,SUN Z G,WANG Q X,et al. Evaluating MODIS phenology product for rotating croplands through ground observations[J]. Journal of Applied Remote Sensing,2013,7(1):073562 doi: 10.1117/1.JRS.7.073562
[20] WU W,XIN Q C. A global annual vegetation phenology dataset derived from GIMMS LAI 3G time series for 1982–2015[C]//IGARSS 2022 - 2022 IEEE International Geoscience and Remote Sensing Symposium. Kuala Lumpur,Malaysia:IEEE,2022:6130
[21] 刘威,吕洁. 傅里叶变换与离散傅里叶变换相结合的信号谱分析教学新方法探讨[J]. 黑龙江科学,2022,13(3):24 doi: 10.3969/j.issn.1674-8646.2022.03.010 [22] LI J,LU X H,JU W M,et al. Seasonal changes of leaf chlorophyll content as a proxy of photosynthetic capacity in winter wheat and paddy rice[J]. Ecological Indicators,2022,140:109018 doi: 10.1016/j.ecolind.2022.109018
[23] CROFT H,CHEN J M,ZHANG Y. Temporal disparity in leaf chlorophyll content and leaf area index across a growing season in a temperate deciduous forest[J]. International Journal of Applied Earth Observation and Geoinformation,2014,33:312 doi: 10.1016/j.jag.2014.06.005
计量
- 文章访问数: 184
- HTML全文浏览量: 35
- PDF下载量: 9
