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]}。

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

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

V(Y)=∑i=1kVi+∑i=1k∑j=i+1kVij…+V1, …, k(1)

式(1)中, V(Y)为模型输出的总方差, k为模型输入参数个数, V_i为模型参数x_i的方差贡献, V_ij为参数x_i通过参数x_j贡献的方差(耦合方差), V_{1, …, k}为参数为参数x_i通过参数x₁, x₂, …, x_k共同作用贡献的方差。

参数x_i的一阶敏感性指数S_i可以表示为

Si=ViV(Y)(2)

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

Sij=VijV(Y)(3)

Sijm=VijmV(Y)(4)

Sij…k=Vij…kV(Y)(5)

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

STi=Si+∑j≠iSij+∑m≠j≠iSijm+Sij…k(6)

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

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

表1

Table 1

表1(Table 1)

表1 PROSPECT-D模型参数的设置 Table 1 PROPSPECT-D parametersin the sensitivity analysis 结构参数 叶绿素a+b /(μ g· cm-2) 类胡萝卜素 /(μ g· cm-2) 花青素 /(μ g· cm-2) 等效水厚度 /cm 干物质 /(g· cm-2) 最大值 1 0 0 0 0.001 0.001 最小值 4 100 50 50 0.05 0.05 基准值 1.5 30 10 5 0.02 0.02

表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值作为最终的样本数。

植被色素的遥感反演模型大致可分为两类: 经验模型和物理模型。 物理模型是严格按照电磁波谱和辐射传输特性经过推导得到的机理模型, 物理意义明确; 经验模型是以统计学为基础的数据驱动模型, 模型的精度和稳定性依赖于数据的质量和数量。 相对于物理模型, 经验模型机理性不足, 其优势在于参数少、 简便、 适用性强。 如果物理模型的参数不能被准确描述, 物理模型的效果反而不如经验模型^[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)中, λ ₁为Anth的敏感波段, λ ₂为Anth不敏感但是对其他色素敏感的波段, λ ₃不受任何色素吸收的影响, 主要说明叶片组织结构对叶片反射的影响。

色素主要影响叶片在可见光波段的光学特征, 色素总敏感性指数随样本数量个数的变化(1~10 000)特征如图1(a, b, c)所示。 Anth总敏感指数S_{Ti_Anth}在波长大于689 nm时恒为0, 即与样本数量大小无关, 样本数量的变化主要影响400~689 nm范围的S_{Ti_Anth}, 在小样本时, S_{Ti_Anth}随样本数的变化强度大于大样本。 如当样本数由1增加到10时, S_{Ti_Anth}明显增加, 尤其是在548 nm处, 增加量达到94.77%; 样本数从10增加到100, S_{Ti_Anth}在400~520 nm范围增加, 而在521~606 nm区间, 情况则相反。 随着样本数的逐渐增加, 不同样本条件下得到的结果之间差异逐渐缩小。 如当样本数从100到1 000, 平均每单位波长(nm)S_{Ti_Anth}的变化仅为236.56%。 当样本数量达到1 000以后, 不同样本数的S_{Ti_Anth}曲线开始聚集, 变化越来越小。 统计发现, 样本数量从1 000增加到2 000、 2 000增加到3 000, 所有波长(400~689 nm)S_{Ti_Anth}的累计变化量 ∑λ=400689ΔSTi_Anth分别为65.99%、 27.85%, 此后, 样本数每增加1 000, ∑λ=400689ΔSTi_Anth不超过11.02%, 平均每单位波长(nm)的变化仅为0.038%。 对于Cab和Car, 它们的总敏感性指数随样本数的变化与Anth类似, 在样本数小于1 000时, 总敏感性指数波动较大, 随着样本数的逐渐增加, 总敏感性指数趋于稳定。 意味着, 样本数为3 000时, 花青素总敏感性指数已经稳定, 与Morcillo-Pallaré s等的研究结果相似。

图1

Fig.1

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

在样本数为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

Fig.2

	Figure Option ViewDownloadNew Window
	图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

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

叶片尺度上, 叶片组分对反射率的影响范围和程度存在显著差异, 在同一波谱区间, 叶片反射率可能对多个参数敏感。 寻找参数的敏感波段有助于参数波谱指数的设计, 最大程度减少反演过程中其他参数的干扰, 降低参数反演的不确定性。 对于一个参数, 最佳敏感波段的反射率应对该参数敏感而对其他参数不敏感。 根据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

Fig.4

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

图5

Fig.5

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

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

图6

Fig.6

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

图7

Fig.7

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

相比于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)

Anth=0.9732×R560-1-4.3739(9)

Anth=1.9181×mARI-0.8081(Sentinel⁃2)(10)

Anth=1.0099×ARI-1.1316(Sentinel⁃2)(11)

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