基于Sentinel-2影像与PROSAIL模型参数标定的玉米冠层LAI反演

苏 伟， 邬佳昱， 王新盛， 谢茈萱， 张 颖， 陶万成， 金 添

1. 中国农业大学土地科学与技术学院， 北京 100083 2. 农业部农业灾害遥感重点实验室， 北京 100083

引 言

叶片是玉米光合作用的主要器官， 玉米的干物质积累大多来自叶片， 因此准确描述玉米叶片生长状态对于玉米长势监测、 灾害胁迫监测、 产量预测等具有重要意义[1]。 叶面积指数是单位土地面积上的植物叶片表面积的一半[2]， 与植被各生理生化过程有密切联系， 是表征玉米长势的一种重要参数。 遥感数据具有空间上的全覆盖性、 时间上的连续性和数据类型的多样性， 能够以不同的时空尺度获取多种作物冠层信息。 基于遥感影像反演玉米冠层叶面积指数， 可为玉米长势监测、 病虫害胁迫监测、 产量预测提供有效依据[3]。

PROSAIL辐射传输模型综合考虑了作物冠层结构、 生长状况以及遥感观测环境的影响， 能够准确模拟作物冠层反射率， 由于其高敏感性与完善的辐射传输机理， PROSAIL模型被广泛应用于叶面积指数(leaf area index, LAI)反演中。 近年来， 众多学者针对PROSAIL模型的LAI反演进行了诸多研究， 反演方法包括查找表法[4]、 经验线性回归、 机器学习[5]以及模型耦合[6]等， 取得了显著成果。 其中查找表法是传统的反演算法， 通过建立大量情况下模型参数与冠层反射率之间的映射关系， 使用代价函数进行迭代运算， 将计算所需的时间转移至反演前， 因此具有广泛的应用[7]。 机器学习等方法则弱化了模型的机理过程， 训练模型更为灵活， 近些年发展迅速。 然而， 由于PROSAIL模型的未知参数较多且测量费时、 费力， 无论是经验统计方法或机理方法， 均鲜有PROSAIL模型参数标定的应用， 在反演过程中缺少作物参数先验信息， 导致参数的取值范围宽泛， 反演的计算量较大， 部分反演结果出现不符合实际的异常值[8]。

所谓参数标定， 就是调整模型的参数范围使其适用特定的研究区和作物类型[9]。 对于模型参数标定， 传统方法是使用试错法或设置经验值， 设定的参数取值范围较为粗略， 难以表达研究区内作物参数的时空差异。 此外， 模型中存在着不同参数可能会带来同一模拟结果的“异参同效”现象， 在反演计算中将不符合实际情况的参数作为计算结果， 导致反演结果与实测值的偏差较大。 贝叶斯理论中马尔可夫链蒙特卡洛方法能够通过构建模型参数的马尔科夫链， 不断进化使模型模拟结果达到观测值的附近， 获取观测值不确定性范围内的模型参数后验分布， 其在各类模型的参数标定中已经有了广泛应用。

为优化LAI反演过程提高反演精度， 本研究对如何进行PROSAIL模型参数标定进行了探索研究。 在利用Sobol法进行PROSAIL模型参数敏感性分析基础上， 使用马尔可夫链蒙特卡洛方法(Markov Chain Monte Carlo， MCMC)方法对PROSAIL模型进行参数标定， 为研究区内不同地块的玉米种植区提供准确的模型输入参数， 并应用于LAI反演中， 通过参数标定优化前后的LAI反演精度比较分析验证基于MCMC方法参数标定的可行性及有效性。

1 研究区概况与数据源

1.1 研究区概况

研究区位于河北省廊坊市安次区， 地理范围为116°35′—116°51′E， 39°20′—39°33′N， 地理位置如图1所示。 研究区地处中纬度地带， 属暖温带大陆性季风气候， 主要地形为平原， 平均海拔约13 m， 年平均气温8～19 ℃， 年平均降水量550 mm。 主要农作物有冬小麦、 夏玉米、 春玉米、 马铃薯等， 研究区夏玉米的生长时期一般为6月底至9月中下旬。

图1 研究区位置及采样点空间分布

1.2 数据源

1.2.1 Sentinel-2A影像

Sentinel-2环境监测卫星是欧盟委员会(EC)和欧洲航天局(ESA)共同倡议的全球环境与安全监测系统—“哥白尼计划”中的第二颗卫星， Sentinel-2A、 B星分别于2015年6月23日、 2017年3月7日发射成功。 Sentinel-2A遥感影像共有十三个波段， 包括可见光、 近红外、 红边、 水汽、 卷云和短波红外波段， 各波段空间分辨率分别为10， 20和60 m， 影像产品为经大气底部反射率校正的1C级影像， 可以通过ESA的数据共享网站下载(https://scihub.copernicus.eu)。 Sentinel-2A卫星时空分辨率较高， 具有红边波段， 红边波段是指植被可见光谱与近红外光谱之间反射率增大最快的波段， 与植被覆盖度、 LAI及各类参数密切相关， 被广泛用于农作物长势监测中。 本研究用于模型参数标定及LAI反演的波段为可见光和近红外波段， 空间分辨率为10 m， 波谱范围如表1所示。 影像获取时间2019年8月18日， 影像获取时太阳天顶角为30.38°， 太阳方位角为144.1°， 平均观测天顶角为8.15°， 平均观测方位角为104.08°。

表1 Sentinel-2A影像部分波段信息

1.2.2 地面测量数据

与遥感影像对应的地面测量时间为2019年8月16日至22日， 由于不同地区的种植时间差异， 玉米的生育期为抽穗期至乳熟期。 玉米冠层LAI的测量仪器为美国LI-COR公司生产的植物冠层分析仪LAI-2200C(LI-COR, Lincoln, NE, USA)， 仪器设置为一个天空光、 四个测量目标值， 探头佩带45°张角的镜头盖， 四个叶下测量点分别位于垄上、 两垄之间1/4处、 两垄中间、 距离垄3/4处。 LAI测量的同时使用中绘i80智能 RTK测量系统进行定位样点。 在研究区内不同地块设置128个采样点， 每个采样点采集两次， 取每一地块内所有LAI的均值作为该地块的测量结果。

2 模型与方法

2.1 PROSAIL辐射传输模型

PROSAIL模型由PROSPECT叶片光学模型和植被冠层二向反射率SAIL模型耦合而成， 是目前应用最广泛的植被冠层辐射传输模型之一。 PROSPECT模型将叶片结构假设为多层表面粗糙的平板， 假设光线各向同性平行， 通过输入叶片理化参数， 模拟光线在平板之间的朗伯散射， 得到400～2 500 nm叶片的反射率与透射率[10]。 SAIL模型假设植被冠层是均一的、 无限延展的混合介质， 并且叶片各向同性， 以辐射传输方程为理论依据， 模拟植被冠层内部辐射传输过程， 得到冠层尺度的反射率。 PROSAIL模型将PROSPECT模型输出的叶片反射率透射率作为SAIL模型的输入进行模型耦合， 考虑了植被对太阳辐射的吸收、 二向反射、 叶片反射透射率、 叶片结构参数、 地表状况等因素的影响， 能较为真实地模拟植被冠层情况， 并且具有很高的稳定性， 经常用于植被参数的遥感定量反演中。

PROSAIL模型的输入参数共包括四类： 叶片理化参数、 冠层结构参数、 地面参数以及角度信息。 用于描述叶片生理生化属性的参数有叶片结构系数N、 叶绿素含量cab(μg·cm-2)、 类胡萝卜素含量car(μg·cm-2)、 干物质含量cm(g·cm-2)、 等效水厚度cw(cm)； 冠层结构参数包括LAI和叶倾角分布参数Lidfa； 地面参数包括热点系数hspot和土壤光谱； 角度信息包括遥感影像拍摄时的太阳天顶角tts、 观测天顶角tto和相对方位角psi。

2.2 MCMC参数标定方法

贝叶斯理论可以根据模型参数的先验知识和模型输出的相应观测值， 实现模型参数的后验估计。 马尔可夫链蒙特卡洛方法在贝叶斯理论框架下， 将马尔科夫(Markov)过程引入到通过计算机进行模拟的蒙特卡洛方法(Monte Carlo)中， 实现抽样分布随模拟的进行而改变的动态模拟， 构造合适的马尔科夫链进行抽样而使用蒙特卡洛方法进行积分计算， 即马尔科夫链可以收敛到平稳分布。 贝叶斯基本公式如式(1)

(1)

式(1)中，θ为PROSAIL模型的输入参数，y为模型模拟反射率，p(θ/y)为参数的后验概率密度函数，f(y/θ)为观测值似然函数，g(θ)为参数的先验分布。

本研究中采用差分进化马尔科夫链算法(differential evolution Markov chain, DE-MC)， 通过多条马尔科夫链的并行更好地搜索参数空间， 并且为提高采样效率， Braak等加入snooker更新算法部分代替平行方向上的更新， 克服原DE-MC算法并行链数必需大于空间维数的限制， 该算法在模型参数标定中已有应用[9]， 算法的具体原理见参考文献[11]。

2.3 Sobol模型敏感性分析方法

敏感性分析是指定性或定量地描述模型参数变化对输出结果的影响程度， 是模型应用过程中的基础和重要环节。 模型的敏感性分析根据作用范围可分为局部敏感性分析和全局敏感性分析， 全局敏感性分析方法除了分析单参数取值变化对模型结果的影响外， 将各参数间的相互作用也考虑进来， 分别得到参数的各阶及全局敏感性， 常用的敏感性分析方法有E-FAST， GLUE， Sobol和Morris等。

Sobol法是一种基于方差的典型全局敏感性分析方法， 本研究采用Sobol法进行PROSAIL模型主要参数的敏感性分析。 Sobol法通过输出结果的方差分解和在输入参数范围内的随机抽样取值， 通过输入参数对输出方差的贡献定量测量参数敏感性[12]。 本研究使用python下SALib敏感性分析环境提供的Sobol法进行PROSAIL模型的敏感性分析， 目标参数Y为Sentinel-2A四个波段波长的模拟反射率。

2.4 LAI反演方法

本研究中PROSAIL模型进行LAI反演的方法为查找表法， 查找表法是通过建立离散的参数与反射率的对应关系表， 使用代价函数遍历迭代计算得到目标参数， 适用于参数多且复杂时无法求出参数反函数的情况。 根据PROSAIL模型参数敏感性分析结果， 将敏感参数作为查找表中的可变参数， 通过设定可变参数的取值范围和变化步长， 其余参数设为定值， 将模型多次运行的输入与输出进行聚合建立LAI反演的查找表。 根据所建立的模型输入参数-冠层反射率查找表， 逐像素分析四个波段反射率与查找表中反射率的匹配情况， 寻找像元反射率与模拟反射率最接近的一组参数。 研究中利用小二乘法构建代价函数， 计算模拟反射率与影像反射率之间的绝对差值和并最小化， 代价函数的建立方法如式(2)

(2)

式(2)中，n为参与反演的影像波段数，Rmodel为模型模拟反射率，RS2A为Sentinel-2A影像反射率， 通过多次迭代计算得到X取最小值时对应的LAI值。

3 结果与讨论

3.1 PROSAIL模型的敏感性分析

PROSAIL模型进行敏感性分析， 首先需要明确输入参数的取值范围， 模型中的角度参数根据影像头元文件得到， 叶倾角分布函数通过地基激光雷达数据进行构建， 本研究中采用已经建立的对应生育期玉米冠层叶倾角函数[13]， 对余下的六个主要参数进行敏感性分析， 各参数范围根据经验值及相关文献获取如表2所示。

表2 PROSAIL模型主要参数取值范围

使用SALib中的Saltelli采样器， 采样次数N为2000， Saltelli采样器生成的参数样本数由采样次数N和参数个数D共同决定， 在只输出一阶敏感性和全局敏感性的情况下， 样本数n=N×(D+2)。 通过每次在参数范围内的随机采样， 计算各参数对模型模拟光谱中波长为490， 560， 665和842 nm的反射率的敏感性大小， 敏感性分析结果如图2所示。

从图2(a—b)中可以看出， 各参数对各波段(波段2， 波段3， 波段4， 波段8)反射率的敏感性各不相同， Sobol法认为敏感性指数>0.05的参数对结果表现出敏感性， 敏感性指数>0.2的参数为强敏感参数， 全局敏感性由于考虑了参数间的相互作用， 一般高于一阶敏感性指数。 从分析结果可以看出： LAI的敏感性最高， 对各波段都表现为强敏感性， 对植被冠层反射率的影响程度最大，N对绿波段和近红外波段敏感， 且在绿波段表现出强敏感性[图2(b)]。 cab的敏感性主要体现在绿波段附近， 对其余波段影响程度较小， 而另一色素参数car对各波段均不敏感。 cm和cw分别反映玉米生长过程中干物质和水分的积累情况， 对可见光波段不敏感， 在近红外波段[图2(d)]仅有cm表现出敏感， 且敏感程度较低。 因此， 基于模型敏感性分析结果， 在LAI反演的查找表中， 将LAI，N和cab作为可变参数， 变化步长分别LAI为1，N为0.1， cab为2， 其余参数根据经验值在建立查找表过程中设为定值。

图2 PROSAIL模型敏感性分析结果(波段2， 3， 4， 8)

3.2 基于MCMC的参数标定结果与分析

鉴于图2的PROSAIL模型的敏感性分析， 由于car和cw对研究中所用的波段范围均不敏感， 在本研究中不进行标定， 仅对LAI，N， cab和cm四个输入参数进行标定。 除LAI的标定范围根据历年采样数据及生育期经验值设定外， 其余参数的先验取值范围与查找表中相同。 本研究认为冠层反射率服从以影像观测值为期望的高斯分布， 使用取对数计算建立似然函数

logLref=-0.5(x-xobs)TΣ-1(x-xobs)-

0.5Klog(2π)-log(detΣ)

(3)

式(3)中，L为似然函数的返回值； T为矩阵转置值； 矢量x和xobs分别表示波段反射率的模型模拟值和观测值，Σ表示反射率观测值的协方差矩阵， 本研究中认为波段反射率观测值相互独立， 方差为观测值的5%；K为空间维数， 即反射率观测值的个数； detΣ表示Σ的行列式值。

算法运行过程中并行链数为4， 每5次采样进行一次链更新， 每10 000次迭代通过比差法进行一次收敛性判断， 计算诊断指标R， 当连续超过3次诊断指标R<1.03， 则认为马尔科夫链达到收敛， 收敛后舍弃前期样本， 保留参数的后验样本。 本研究中的采样点共分布在25个地块中， 理论上同一地块仔的玉米生长状态相似， 因此以地块为单位进行参数标定， 取地块内所有玉米冠层反射率的平均值为该地块参数标定的输入， 通过参数马尔科夫链的进化达到收敛， 得到各参数的后验分布与范围。 图3中列出了两次地块参数标定后各参数的先验后验分布， LAI和Cab在两次定标中都得到高斯分布的定标结果。

注： X轴表示参数取值， Y轴表示概率密度的相对大小， 绿色柱状图表示参数后验概率分布， 面积和为1

从图4(a,b)中可以看出， 各参数先验分布为取值范围内的均匀分布(蓝色线)， 经参数标定后各参数的取值及概率分布均发生了较大变化。 两次参数标定后， 各参数的取值范围均得到缩小， 的确只有LAI和cab在新的取值范围内呈现出高斯分布。N和cm两参数在实际测量中难以获取， 给反演中的参数设置带来了困难， 参数标定后两参数的取值集中在右边界， 说明对于玉米及这一生育期而言，N和cm的取值波动较小， 集中在某一定值附近。 参数标定根据各地块的反射率观测值， 通过马尔科夫参数链的进化得到了观测值不确定性范围内的参数取值分布， 并提供了每一参数概率最大的取值， 由此各地块的玉米生长差异得以表达， 参数的后验取值范围替代原先LAI反演时查找表内参数较为宽泛的取值， 提高同一地块内LAI反演的稳定性， 避免出现极端值的情况。

3.3 模型参数标定对LAI反演精度的优化分析

根据2.1中的敏感性分析结果， 首先建立LAI反演的参数-反射率查找表， 各参数的取值范围为默认， 输入Sentinel-2A影像的观测反射率， 使用最小化代价函数进行计算反演得到各采样点的LAI。 其次， 根据各地块参数标定结果， 使用后验分布范围优化参数设置， 进一步缩小各地块LAI的可能取值范围建立新的查找表， 并使用每一地块标定后验分布中LAI的最高概率密度取值作为该地块的参考值， 与查找表反演的LAI结果各取50%的可信度进行计算得到新的LAI反演结果。 分别将优化前的LAI反演和加入参数标定结果优化的LAI反演结果与实测值进行精度评价， 评价指标包括决定系数R2、 均方根误差RMSE、 平均偏差Bias和估算精度EA， Bias和EA的计算方法如式(4)和式(5)

(4)

(5)

式(4)和式(5)中，yg和ys分别表示参数的模型模拟值和实测值， Mean表示实测值的均值， 提取对应实测点位置的LAI反演结果与实测数据进行精度验证， LAI反演的精度验证结果如图4(a,b)和表3所示。

图4 参数标定前后LAI反演精度比较分析

表3 参数标定前后LAI反演精度验证结果

从图4(a,b)中可以看出， 两次反演的LAI与实测LAI值拟合效果较好，R2分别达到0.73和0.76， 反演结果没有明显的系统偏差， 拟合线与1∶1线较为接近， 证明查找表反演方法的有效性。 传统的查找表法进行反演时， 查找表中参数的范围通常根据经验值或默认范围设定， 是一种先验信息较少的反演方法， 参数范围较大导致反演计算量的加大以及异常值的出现。 参数标定根据观测反射率获取不确定性范围内的参数取值， 排除不合先验的取值， 缩小了各地块内LAI的可能范围， 优化了原查找表反演中的设定， 使得各地块玉米生长状况的差异得以表达， 加快反演速度的同时有效减少了偏差较大的病态反演出现。 表3列出了加入参数标定结果前后的反演精度差异， RMSE表示反演结果的离散程度， 参数标定优化前后的RMSE分别为0.79和0.32， 有较为明显的提高， 反演结果均匀地分布在一阶拟合线附近。 Bias是反演值与实测值的平均偏差， 描述反演结果偏大或偏小的程度， 在实际应用中这一指标同样较为重要， 原反演结果的Bias为0.67， 约为实测均值的20%， 参数标定优化后的Bias降低至0.26， 约为实测均值的8%， 同时估算精度提升至90%， 体现了参数标定在LAI反演中的优化效果， 有效地降低了LAI反演值的偏差， 满足实际应用需求。

4 结 论

以Sentinel-2A影像为数据源， 首先使用Sobol方法对PROSAIL模型进行参数敏感性分析， 量化各参数对可见光及近红外波段的敏感性， 为LAI反演提供基础， 其次基于MCMC方法对PROSAIL模型进行了参数标定， 目的是提供区域内玉米各参数的取值范围和分布， 获取更为丰富的参数信息， 最终对比分析了参数标定优化前后LAI反演精度差异， 探索参数标定在LAI反演中的应用潜力， 主要研究结论如下：

(1)PROSAIL模型中LAI对可见光和近红外光均表现出强敏感性， 满足反演的需求， 同时叶绿素含量Cab对绿波段表现出较强敏感性， 结构系数N对绿波段和近红外波段较为敏感， 将这三个参数作为LAI反演的查找表中的可变参数， 其余主要参数对可见光和近红外波段不敏感， 在反演过程中设为固定参数。

(2)基于MCMC方法的参数标定能够获取玉米生长区域内各参数的后验概率分布， 表达区域间的玉米生长状况差异， 缩小区域内LAI反演时中各参数的取值范围， 为参数反演提供先验信息。 参数标定优化后的查找表有效降低了LAI反演偏差， RMSE和Bias分别降低了0.47和0.41， 估算精度由原先的76%提升至90%， 提高了基于PROSAIL模型的LAI反演在实际应用过程中的准确性， 能够为作物长势监测提供更为准确的数据支持。

本研究使用MCMC方法对PROSAIL模型进行过了参数标定， 对LAI的反演精度有所提高， 但是参数标定的应用仅限于为查找表及反演提供参考。 在今后的研究中， 将着重挖掘参数标定在模型中的应用潜力， 与机器学习等方法结合， 进一步优化参数反演过程。