花青素对植物反射特性的影响及遥感估算: 叶片尺度

梁守真,隋学艳,王 猛,王 菲,韩冬锐,王国良,李洪忠,马万栋

1. 山东省农业科学院农业信息与经济研究所,山东 济南 250100 2. 中国科学院深圳先进技术研究院先进计算与数字工程研究所,广东 深圳 518055 3. 生态环境部卫星环境应用中心生态保护红线监管中心,北京 100029

引 言

植被是陆地土壤、 大气的重要“纽带”,植被生化组分直接或间接地参与生物地球化学循环,影响和控制着生态系统生态与物理过程,尤其是植物色素,其动态变化直接影响着碳、 氮和水循环,在生物圈中具有极其重要的作用。 因此,定量估算植被色素含量成为地球科学、 生物科学领域的重要研究内容之一。 传统的植被色素含量测定主要通过在实验室内对植物样本进行化学分析来实现,传统方法能够得到比较准确的色素含量,但获取信息有限(点信息),耗时,对植株有一定的破坏性,难以有效跟踪植株随环境变化的特征。 遥感将传统的点测量信息扩展为连续的面信息,且不会对生态系统造成破坏,省时、 省力,已被广泛用于估算叶片色素,但主要集中于叶绿素[1-2]。 可能来自两方面的原因,一是叶绿素与氮关联性强,关注度高,同时叶片叶绿素含量高,对植被可见光波谱影响大; 二是在早期辐射传输模型中,叶绿素便是一个必要的输入参数,遥感监测机理相对清楚,反演精度高。 在全球气候变化的大环境下,除了叶绿素,需要更多植被色素(如花青素、 类胡萝卜素)的动态来说明与植被相关的各种环境问题,比如环境胁迫、 植物病害等。 花青素是自然界三大主要植物色素之一,具有调节光合作用、 保护光合组织免受强光伤害、 避免环境胁迫(如冷冻、 污染、 植物病原体)的能力[3-4],其含量与植物的生长发育、 生长状态以及生长环境状况密切相关,被认为是理解环境胁迫或季节波动引起的植物生理反应与抗性的关键[5-6],花青素的定量监测可提供植物对环境胁迫的响应与适应信息。 由于相关研究的缺乏,花青素对叶片光学特性的影响特征仍不清晰。 随着包含花青素参数的辐射传输模型PROSPECT-D的出现,从机理上开展花青素的光学影响研究成为可能。 研究以PROSPECT-D为基础,采用敏感性分析方法,系统量化花青素的光学影响,从叶片尺度探讨花青素对叶片光学属性的支配作用,筛选花青素敏感波段,并初步探讨花青素的遥感反演方法,为更广泛地开展花青素遥感研究提供理论支撑。

1 实验部分

1.1 叶片辐射传输模型

POSPECT是Jacquemoud和Baret在Allen平板模型基础上发展起来的辐射传输模型,其通过模拟叶片在400～2 500 nm上行和下行辐射通量来计算叶片半球反射率和透射率[7-8]。 模型假设叶片由N层同性质的平板堆叠而成,由N-1层空气隔开,光线的非漫射特性只存在于最顶层。 在模型中,散射过程用折射指数n和叶片结构参数N来表达,而吸收则表示为叶片组分含量及相应吸收系数的线性组合。 最早的PROSPECT模型只包含三个参数:N、 叶绿素和等效水厚,随后更多的叶片组分(干物质、 类胡萝卜素、 花青素)先后被纳入模型,涌现出多个版本的PROSPECT模型(PROSPECT-4、 PROSPECT-5、 PROSPECT-D)。 PROSPECT-D是第一个同时包含了叶绿素、 类胡萝卜素和花青素的叶片辐射传输模型,相对于以前的版本,其对叶片光学特性的描述更加准确,不确定性降低[9-10]。

1.2 参数敏感性分析

敏感性分析是定性或者定量地分析评价模型参数对模型结果影响程度的一种方法。 通过调整模型输入参数来确定参数变化对模型输出的影响程度和范围,识别影响模型输出的关键参数,敏感性分析对模型的构建和优化有非常重要的意义[11-13]。 敏感性分析包括局部敏感性分析和全局敏感性分析,局部敏感性分析只计算单个参数变化对模型输出的影响,不考虑参数之间的相互作用; 全局敏感性分析考虑了参数之间的耦合作用,分析多个参数协同变化对模型输出的影响[14-16]。 相对于局部敏感性分析方法,全局敏感性分析更能全面地反映参数对模型输出的影响。

典型的全局敏感性分析方法包括FAST法[17]、 EFAST[18-19]、 Sobol法[20]、 改进的Sobol’法[21]等。 Sobol’法是最具有代表性的全局敏感性分析算法,由Sobol’提出,Saltelli等对该算法做了一定的改进[22]。 与原算法相比,改进的Sobol’算法具有更强的稳定性和更小的计算代价[22-23]。 Sobol’法基于分解的思想,将模型输出的总方差表示为单个参数和组合参数的函数

(1)

式(1)中,V(Y)为模型输出的总方差,k为模型输入参数个数,Vi为模型参数xi的方差贡献,Vij为参数xi通过参数xj贡献的方差(耦合方差),V1,…,k为参数为参数xi通过参数x1,x2,…,xk共同作用贡献的方差。

参数xi的一阶敏感性指数Si可以表示为

(2)

式(2)中,Si是参数xi对模型输出总方差的直接贡献,而参数xi的二阶、 三阶以及高阶敏感性指数Sij、Sijm、Si,…,k则反映了参数之间相互作用对模型输出的影响,其表达式分别为式(3)—式(5)

(3)

(4)

(5)

参数xi的总敏感性指数STi为各阶敏感性指数之和,反映的是参数xi的直接贡献及与其他参数的交互耦合作用间接对模型输出方差的贡献,其表达式为式(6)

(6)

若参数之间无耦合作用,则高阶敏感性指数为0,Si与STi值相等。 敏感性指数越大,表示参数对模型输出的贡献越多,模型对参数越敏感,反之亦然。 考虑到参数之间存在交互影响,各参数的STi之和会出现大于1的情况,故而对其进行归一化处理。

1.3 参数设置

模型参数敏感性分析需要设定模型参数上下边界和样本数。 PROSPECT-D模型输入参数叶绿素(Cab)、 胡萝卜素(Car)、 等效水厚(Cw)、 干物质(Cm)等生化组分含量边界主要根据LOPEX’93[24]和ANGERS数据集[25]来设定,花青素取值范围参考Gitelson等的研究[26](表1)。

表1 PROSPECT-D模型参数的设置Table 1 PROPSPECT-D parametersin the sensitivity analysis

改进的Sobol’法采用蒙特卡洛法对参数空间进行采样,样本数n与模型输入参数的个数k共同决定了敏感分析所需的运算次数:n(k+2)次[27-28]。 大样本运算代价高,但小样本生成的结果缺乏稳定性。 在一定程度上,n是分析精度和处理效率之间的权衡。 由于模型参数的数量和复杂程度不一,不同的研究采用了相异的n值[15,23,29]。 针对PROSPECT-D模型,研究采用逐步计算的方法,从小样本开始逐渐增加样本数量,将敏感性指数值开始达到稳定时的n值作为最终的样本数。

1.4 花青素反演方法

植被色素的遥感反演模型大致可分为两类: 经验模型和物理模型。 物理模型是严格按照电磁波谱和辐射传输特性经过推导得到的机理模型,物理意义明确; 经验模型是以统计学为基础的数据驱动模型,模型的精度和稳定性依赖于数据的质量和数量。 相对于物理模型,经验模型机理性不足,其优势在于参数少、 简便、 适用性强。 如果物理模型的参数不能被准确描述,物理模型的效果反而不如经验模型[29]。 主要讨论两种类型的参数——窄波段反射率和多波段组合的波谱指数反演Anth的可行性,并比较其差异。 首先基于敏感性分析筛选Anth的敏感波段,然后挑选不同波段的反射率构建具有一定物理指示意义的波谱指数,最后分别统计高光谱窄波段、 波谱指数与Anth含量之间的相关关系,建立Anth反演模型。 以波谱指数为参量的经验模型有时被称为半经验模型,研究主要考虑两种光谱指数与Anth的关系: ARI(anthocyanin reflectance index)和mARI(modified anthocyanin reflectance index)[30]

ARI=ρ(λ1)-1-ρ(λ2)-1

(7)

mARI=[ρ(λ1)-1-ρ(λ2)-1]×ρ(λ3)

(8)

式(7)和式(8)中,λ1为Anth的敏感波段,λ2为Anth不敏感但是对其他色素敏感的波段,λ3不受任何色素吸收的影响,主要说明叶片组织结构对叶片反射的影响。

2 结果与讨论

2.1 样本数量对敏感性指数的影响

图1 样本数对色素总敏感性指数的影响Fig.1 The impact of the number of samples on pigment global sensitivity index

2.2 叶片反射率对花青素的敏感性

在样本数为3 000时,PROSPECT-D模型参数在400～2 500 nm 范围的总敏感性指数如图2所示。 Anth对叶片反射率的影响范围为400～689 nm,在509 nm处有最高的总敏感性指数,86.64%; 在428～509 nm区间,Anth总敏感性指数随着波长的增加而增加,而在400～428和509～689 nm区间,情况则相反。 400～689 nm范围内,叶片反射特性受多个组分的影响,Cab、 Car、 Anth、 Cm、 Cw的累积总敏感性指数分别为9 386.16%、 2 844%、 12 284.85%、 606.82%、 0.000 065%。 Cab主要影响400～780 nm范围的叶片反射,在红外波段叶片反射特性不受Cab变化的影响。 在400～780 nm范围内,Cab总敏感性指数存在两个显著的波峰: 441 nm(24.16%)、 679 nm(88.74%),524 nm处总敏感性指数最低(0.09%)。 在598～689 nm区间,Cab比Anth有更高的敏感指数。 Car影响400～559 nm的叶片反射率,当波长大于560 nm时,Car总敏感性指数为0,在417～451 nm波段,Car比Anth有更高的敏感指数; Cm对叶片反射率的影响主要集中于400～432和700～2 500 nm,在可见光波段,其总敏感性指数随着波长的增加从迅速降低,432 nm以后,Cm总敏感性指数小于1%; Cw主要影响波长大于1 100 nm范围的反射率,在400～689 nm范围Cw的总敏感性指数非常小,对叶片反慑率的影响可忽略; 结构参数N主要影响红外波段的叶片反射率,在400～689 nm,平均总敏感性指数仅为13.37%。

图2 PROSPECT-D参数的总敏感性指数Fig.2 Total sensitivity indices of PROSPECT-D parameters along the 400～2 500 nm spectral range

局部敏感性分析可认识单一因素变化条件下模型输出的变化规律。 在其他参数设为基准值时,将Anth浓度从0 μg·cm-2按照5 μg·cm-2的步长逐渐增加到50 μg·cm-2,叶片反射光谱随Anth的变化特征如图3所示。 由图3可以看出,Anth的变化仅仅影响400～689 nm的叶片反射,尤其500～650 nm。 当叶片Anth含量为0时,即叶片无花青素时,叶片有最高的反射率,在531 nm有显著的波峰。 随着叶片中Anth的出现,叶片反射率开始降低,并且Anth浓度越高,叶片在400～689 nm范围的反射率越低,反射波谱曲线越平缓,反射峰值向长波方向移动。 Anth浓度为5 μg·cm-2时,反射峰位置和反射率分别为530 nm、 8.95%,当其浓度升高为15 μg·cm-2时,峰值降低为5.17%,峰值位置将移至602 nm处,Anth浓度到达50 μg·cm-2时,反射峰值仅为4.36%(633 nm处)。

图3 不同花青素浓度下叶片反射率Fig.3 Leaf reflectance with different anthocyanin concentrations

2.3 花青素的敏感波段

叶片尺度上,叶片组分对反射率的影响范围和程度存在显著差异,在同一波谱区间,叶片反射率可能对多个参数敏感。 寻找参数的敏感波段有助于参数波谱指数的设计,最大程度减少反演过程中其他参数的干扰,降低参数反演的不确定性。 对于一个参数,最佳敏感波段的反射率应对该参数敏感而对其他参数不敏感。 根据PROSPECT-D模型的敏感性分析结果,Anth对叶片反射率的影响范围窄,仅限于可见光波段(400～689 nm),但在该波段范围,除了Anth外,Cab、 Car、 Cm和N的变化都将引起叶片反射率的波动,尤其在400～405、 417～451和598～689 nm范围,其他参数甚至比Anth有更高的敏感指数。 显然,这些区间不能作为Anth的敏感区。

为了明确Anth的敏感区间,将N、 Cab、 Car、 Cm、 Cw总敏感性指数进行累加,并与Anth的总敏感指数进行比较(图4)。 叶片反射率仅在467～589 nm区间对Anth的敏感度高于其他参数。 这意味着在467～589 nm,Anth的浓度是决定着叶片反射特征的主要因子。 但应注意的是,N、 Cab、 Car的总敏感性指数并不为0,它们的变化仍旧会对叶片波谱产生一定影响。 参数N是模型的一个假设参数,无法实测,主要通过迭代优化或经验计算公式来计算[8,25]。 对于同一种类植物,可认为是定值[31]。 由于所感兴趣的是可测的生化组分,因此主要讨论Cab、 Car对Anth的干扰程度。 在467～589 nm范围内,随着波长的增加Cab的总敏感指数先下降后增加,而Car敏感度主要呈下降趋势。 根据叶片反射对Cab、 Car敏感度的变化特征,对Anth敏感区域进行进一步划分。 在467～505 nm区间,Cab、 Car平均敏感性指数分别为7.13%、 13.97%,其变动对叶片反射率的影响都较大,不可忽视,可认为是Car、 Cab、 Anth共同影响区; 在506～541 nm波谱范围内,Cab总敏感性指数均低于1%,平均为0.32%,意味着Cab的变化对叶片反射仅产生细微的影响,可认为叶片反射特征只受Anth、 Car的影响; 在542～589 nm区间,Car的总敏感性指数小于1%,可视为Cab、 Anth影响区,严格意义上,Car的总敏感系数在波长大于560 nm时才为0。 以505和541 nm为节点,Anth的敏感区被划分为三个子敏感区: 467～505、 506～541和542～589 nm(图5),在这些区间,Anth是支配叶片反射特征变化的最主要的色素,但Cab、 Car在每个区域也会对叶片反射特征产生影响,当利用这些敏感区间的遥感数据反演Anth时,仍需要消除Cab、 Car的干扰。

图4 Anth与其他参数累计总敏感性指数对比Fig.4 The comparison of the total sensitivity index in anthocyanin and other parameters

图5 花青素敏感区间划分Fig.5 Division of sensitive spectral regions of anthocyanin

2.4 叶片花青素反演模型

波谱指数ARI、 mARI的计算首先需要确定Anth的敏感波段和非敏感波段以及不受色素影响的红外波段的波谱范围。 在植被生化参数遥感反演中,Sentinel-2和Landsat8数据是两种较为常用的卫星遥感数据。 Sentinel-2多光谱成像仪(MSI)的B2和B3波段(B2: 458～523 nm,B3: 543～578 nm)处于Anth的敏感区间,在B2波段,Cab、 Car对叶片反射率的贡献较大,而B3波段,只需要考虑Cab的干扰,因此,B3波段反射率更适宜计算ARI、 mARI; B5波段(698～713 nm)反射率对Cab敏感,而对Anth、 Car不敏感,满足要求; B7(773～793 nm)和B8(785～900 nm)为近红外波段,叶片反射率对三种色素都不敏感,主要依赖于叶片密度和厚度,因此选择Sentinel-2 MSI的B3、 B5、 B7(或B8)分别代表来计算ARI、 mARI。 相比于MSI传感器,Landsat8卫星搭载的陆地成像仪(OLI)波宽大,其B2(450～515 nm)、 B3(525～600 nm)反射率对Anth敏感,与MSI相似,B3更适于计算花青素波谱指数,但是Landsat缺少红边波段,只能选择Cab敏感的B4(630～680 nm)。 但在B4红波段,除了对Cab敏感外,Anth也会对叶片反射率产生一定的影响; B5(845～885 nm)对色素不敏感,因此对于Landsat8,研究选择B2、 B3和B4分别对应计算波谱指数。

模型构建所采用的数据为实地采集的51个红端木(Cornus alba)叶片样本的反射率(400～800 nm)和Anth含量数据,植物叶片反射率采用积分球测定,Anth含量则采用分光光度法测定。 首先应用OLI、 MSI波谱响应函数将地面实测叶片高光谱反射率转换为宽波段反射率,在此基础上计算ARI和mARI两个指数。 同时,根据三种色素的敏感性特征,研究选取了四个高光谱窄波段的反射率数据: 509 nm(Anth敏感度最高)、 560 nm(Car不敏感)、 524 nm (Cab不敏感)、 467 nm(Cab、 Car影响大),统计量化Anth与它们之间的关系。

根据图6,在叶片尺度上,Anth与560 nm处的反射率倒数有最高的相关系数(r=0.948),其次为524 nm的反射率倒数,而对Anth最敏感的509 nm叶片反射率倒数与Anth之间的相关系数却仅有0.774 6,Anth与467 nm处的反射率倒数相关性最低。 Anth与ARI、 mARI之间同样存在显著的线性关系,但是Sentinel-2和Landsat8数据有不同的表现[图7(a,b)]。 Sentinel-2 ARI、 mARI与Anth的相关系数分别为0.953 2,0.953 6,线性拟合模型的残差平方和为分别307.97(ARI)、 305.45(mARI)。 无论是相关系数还是残差平方和,mARI稍高于ARI但并不显著。 对于Landsat8 数据,当Anth含量低于10 μg·cm-2时,绿波段反射率高于红波段,导致ARI和mARI为负值。 说明基于Landsat数据的ARI、 mARI难以用于反演花青素浓度较低的叶片。 对Anth浓度高于10 μg·cm-2的叶片样本数据进行单独统计,Landsat8 ARI、 mARI与Anth的相关系数分别为0.804 7、 0.814 7,低于Sentinel-2。 可能来自两方面的影响,

图6 叶片花青素含量与高光谱窄波段反射率的关系Fig.6 Relationship between anthocyanin concentration and highspectral bands

图7 叶片花青素含量与ARI、 mARI的关系(a): Sentinel-2; (b): Landsat8Fig.7 Relationship between anthocyanin concentration and ARI,mARI(a): Sentinel-2; (b): Landsat8

一是Landsat8的波宽大,波谱反射率受多种色素的共同影响,第二个原因则是Landsat8 缺乏仅对叶绿素敏感的红边波段,而替代的红波段反射率又受叶绿素和花青素两种色素的影响,这些因素降低了波谱指数与Anth的相关性。 因此,相比于Landsat8,Sentinel-2数据更适于反演叶片Anth。

相比于ARI,mARI引入近红外波段反射率,同样条件下,其数值范围小于ARI。 近红外波段数据引入的目的是消除叶片结构差异的影响,尽管Anth与ARI、 mARI相关性差异不明显,但Anth总是与mARI有更高的相关系数,这说明近红外波段的引入在一定程度上可提高花青素反演模型的精度。 同时,相对于单一窄波段,基于Sentinel-2的ARI和mARI,与Anth表现出了更好的相关性。 Anth的敏感波段范围内,叶片在高光谱窄波段的反射率仍旧不可避免地受Cab或Car的影响,而对于ARI和mARI,它们通过数学运算的方式最大程度上消除或减弱了其他色素的干扰,使得波谱指数和Anth相关性更好。 无论是波谱指数还是高光谱窄波段,它们与Anth浓度之间的相关系数都可通过0.01水平的显著性检验。 最佳高光谱窄波段和波谱指数与Anth的定量关系可分别表示为式(9)—式(11)

(9)

Anth=1.918 1×mARI-0.808 1(Sentinel-2)

(10)

Anth=1.009 9×ARI-1.131 6(Sentinel-2)

(11)

式中,Anth为花青素浓度(μg·cm-2),R560表示560 nm处叶片反射率。

3 结 论

应用辐射传输模型PROSPECT-D和敏感性分析算法,系统模拟和分析了Anth对叶片反射波谱的影响,确定了Anth的敏感波段,并研究了不同参数反演叶片Anth的可能性,结论如下:

(1)Anth对叶片光学反射特性的影响限定于400～689 nm。 在509 nm处叶片反射率对Anth最敏感,叶片中Anth含量的增加将会降低叶片在可见光波段的反射率。

(2)467～589 nm为Anth的敏感波段,在此区间叶片反射率主要受Anth支配。 根据Cab、 Car对叶片反射特性的影响程度,Anth敏感区间可划分三个部分: 467～505 nm(Car、 Cab、 Anth共同影响区)、 506～541 nm(Anth、 Car影响区)、 542～589 nm(Cab、 Anth影响区)。

(3)叶片尺度上,560 nm是反演Anth的最佳高光谱窄波段,融合了多个波段的ARI、 mARI与Anth有更好的相关性,且应用mARI反演Anth能得到更好的结果,而卫星传感器的波段设置会影响ARI和mARI反演叶片Anth的适用性,波宽更窄、 包含红边波段的Sentinel-2较Landsat 8数据更适于反演叶片Anth。

花青素作为植物三大色素之一,其动态变化监测的重要性不言而喻,发展基于光学的Anth反演方法,对于未来快速便捷非破坏性地获取植被花青素具有重要的实际意义。 然而由于其含量相对较低、 吸收区间与其他色素重叠等问题,从光学角度开展其研究的范围和深度远不如叶绿素。 本研究从辐射传输理论出发,对花青素的光学影响进行研究,确定了花青素的敏感波段,并初步构建了叶片尺度的Anth遥感反演模型,为未来更大范围的花青素遥感估算奠定重要的理论基础。 在花青素的遥感反演模型构建时,采用了一种植物叶片,对于其他植物种类叶片是否存在相似的关系,有待研究进一步开展。