基于MODIS数据的三种模型对区域玉米生产力的估算效果*

钱 娅，郭建茂，2**，李 羚，郭彩云，刘俊伟

（1.南京信息工程大学应用气象学院，南京 210044；2.南京信息工程大学无锡研究院，无锡 214100）

陆地生态系统在地−气系统的物质、能量和动量交换以及全球气候变化中扮演着重要的角色。定量估计陆地生态系统的光合作用，不仅是估计作物产量和生产力的基础，而且对理解区域乃至全球的碳收支和全球气候变化机理及影响具有重要意义[1]。国内外学者已建成不同尺度下，基于遥感参数估算生态系统GPP 的模型，如1972年Monteith[2]利用光能利用率原理，首次提出利用植被吸收的光合有效辐射（APAR）和光能利用率（LUE）估算陆地净初级生产力的方法[3]。20 世纪90年代Gamon 等[4]证明了531nm 的光谱反射率变化能够明显反应叶黄素循环色素之间的相互转换，建立了光化学植被指数PRI（Photochemical Reflectance Index，PRI），并认为叶黄素循环是监测叶片LUE 波动的一种有效手段。1995年美国马里兰大学地理系构建了全部用卫星遥感数据获取APAR 以及光能利用率的GLO-PEM 模型[5]，该模型应用的PAR 数据和气候数据均由卫星遥感观测得到，实现了GPP/NPP 遥感模型的全遥感化[6]。2004年Xiao 等[7]建立了VPM 模型，并利用SPOT-4 卫星的VEGETATION（VGT）数据和MODIS数据基于VPM 模型估算了常绿针叶林、落叶阔叶林、常绿阔叶林和高山草原的GPP[8]。2007年李世华结合GLO-PEM 和VPM 模型，建立了REG-PEM 模型，分别估算了江西省和黑河流域的初级生产力，并对模型参数进行了优化[9−13]。

现有的基于LUE（光能利用率）估算GPP（总初级生产力）的模型中，大多数是利用查表法，确定植被最大光能利用率，进行温度、湿度订正，求算实际LUE 后得到GPP 的值。LUE 是估算GPP 的模型的重要输入，目前的估算方法过于简单粗糙，计算公式复杂且精度较差，可能导致生态系统光能利用率进而影响GPP 的估算精度。

农田作为陆地生态系统的重要组成部分，在全球能量平衡和温室气体的收支方面起到非常重要的作用，研究农田生态系统的碳收支和验证大尺度卫星遥感模型的适用性具有典型代表性。目前，大多研究对玉米农田系统生产力的估算仅仅是基于单一的估算模型，为进一步探究遥感参数估算生态系统GPP 的效果，本研究以GPP/NPP 遥感模型建模思想为基础，将卫星遥感植被指数作为参数引入模型，分别利用APAR 模型、PRI 模型和REG-PEM 模型估算不同时间尺度下玉米农田系统的GPP 值，并通过与锦州生态系统野外观测站观测的GPP 值进行比较，以检验各卫星遥感植被指数及其模型在估算玉米农田生态系统生产力方面的效果。

1 资料与方法

1.1 研究区概况

锦州市地处欧亚大陆东部，属于典型的温带季风型气候，年平均气温为7.8～9.0℃，其中最冷月为1月，平均气温−8.0℃，最热月为7月，平均气温24.4℃，年无霜期为144～180d，年平均降水量567mm。锦州玉米农田生态系统野外观测站位于东北玉米产区带，地处辽宁省锦州市太和区英屯村玉米地（41°49'N，121°12'E，海拔17m），隶属于中国气象局沈阳大气环境研究所，周边为大面积平坦均一的农田。

1.2 研究数据

1.2.1 通量数据

随着微气象理论发展及气象观测仪器、数据采集和计算机存储、数据分析和自动传输等技术的不断进步，涡度相关技术在实际应用中也取得了长足的发展和进步。通量数据主要包括半小时序列的净生态系统碳交换速率（Net Ecosystem Exchange,NEE）、光合有效辐射（Photosynthetically Active Radiation，PAR）、距地面2m 处的空气温度（Air Temperature，Ta）、摩擦风速（u*）等。由于仪器故障、标定、维护或雨雪的影响，观测系统的实测数据存在异常值或缺失值，需要对通量数据进行数据插补、整合等预处理。采用平面拟合法对数据进行坐标旋转，通过WPL 校正[14]消除水热传输造成的CO2通量变化的密度效应，用差分法[15]剔除异常值，利用平均值比较法[16]确定u*C=0.18m·s−1，筛去无效数据。预处理时发现，仅约70%的CO2通量数据是有效数据。因此，分别用平均昼夜变化MDV 插补[15]和Van't Hoff 呼吸方程[17]插补，针对日间和夜间的CO2通量数据进行校正。用插补后的通量数据建立夜间生态系统呼吸与近地层气温间的函数关系，计算式为

式中，Renight为夜间生态系统呼吸速率，以CO2或C 物质的量计，单位为μmolC·m−2·s−1，等于通量站观测的夜间净生态系统CO2交换量NEE（Net Ecosystem Exchange，NEE）；Ta为气温。夜间GPP值为

夜间由于植被无法进行光合作用，观测的NEE（净生态系统CO2交换量）主要为生态系统呼吸Renight，因此夜间GPPnight为0；拟合后得到玉米生态系统中，A 值设为0.1729，B 值设为0.158；假定日间生态系统呼吸（Reday）对气温的响应与夜间相同，可根据式（1）计算日间生态系统呼吸量。日间生态系统光合作用和呼吸作用共存，因此，日间GPPday应为观测的NEE 与生态系统呼吸Reday的合项。其计算式为

1.2.2 MODIS 数据

MODIS 数据具有存储高效、信息丰富、获取速度快、覆盖范围广等优点，对作物生长信息提取有显著优势。MODIS 数据由陆地过程分布式数据档案中心LP DAAC 提供，包括2013−2014年覆盖锦州通量站点L1B 级辐射率数据（MOD/MYD021KM）、8d合成的地表反射率数据（MOD09A1）。由于传感器故障或云雾的影响和输入数据缺失等原因使得MODIS 影像产生的异常值，显示了地物反射率的错误信息，以此错误信息计算的结果将对模型的精度产生一定影响，因此需要对MODIS 数据进行校正、筛选等预处理。利用MCTK 投影插件工具对MODIS数据重投影和几何校正，通过FLAASH 大气校正获得地物反射率和辐射率、地表温度等真实物理模型参数，根据相邻域的像元进行线性插补剔除异常的像元信息。预处理后，以传感器和太阳的几何角度为标准分为全向散射和后向散射，筛选后满足条件的MODIS 反射率影像共38 景，其中后向散射方向的24 景。

1.3 总初级生产力（GPP）估算方法

1.3.1 光合有效辐射吸收项模型（APAR 模型）

植被光合作用的能量来源于太阳辐射，植物对光合辐射PAR（Photosynthetically Active Radiation，PAR）或APAR（光合有效辐射）的截获与利用，转换为自身组织和器官生长发育所需能量和干物质。与PAR 相比，APAR 能够反映植被生物量、繁茂程度和冠层结构等特征，因此针对季节性变化明显的农作物而言，APAR 通过反映作物生物量的变化能够较好地监测作物生产力[18]。利用遥感数据计算的APAR 来估算GPP，可以满足不同区域尺度、不同时间频度的连续监测要求[19]。利用卫星遥感数据监测APAR 是反演生态系统光合潜力、潜在产量、作物生长模拟研究以及全球碳储量的有效手段之一。APAR与PAR 和光合有效辐射吸收比（Fraction of absorbed Photosynthetically Active Radiation，FPAR）之间存在如下关系，即

本研究主要是针对生态系统，故直接使用通量塔观测的PAR 数据，使得参数PAR 的精度更高。式中，FPAR 取决于植被类型和植被冠层特征，而植被冠层是由光合有效植被（PAV，绿叶等，可进行植物光合作用）和非光合有效植被（NPV，衰老的叶片、枝条和茎等，不能进行植物光合作用）两部分组成[20]，因此相应地FPAR 也由两个部分组成，即

其中只有FPARPAV用来计算光合作用。

归一化植被指数（Normalized Difference Vegetation Index，NDVI）对植物的生长状况及生理活性变化的反映具有很强的敏感性，但相关研究证明其易受土壤背景和大气的影响[21]，Huete 等[22]在大气抵抗指数和植被指数的基础上提出了增强型植被指数EVI[23]（Enhanced Vegetation Index，EVI），其计算式为

式中，ρNIR和ρRED分别表示近红外和红光的地表反射率，可利用MODIS 第2 波段（841−870nm）和第1 波段（620−670nm）反射率代替，ρBLUE表示蓝光波段的地表反射率，可利用MODIS 第3 波段（459−479nm）反射率代替；G 表示增益因子，设为2.5；L 为土壤调节参数，设为1；C1和C2为气溶胶抵抗系数，设为6 和7.5。

通过基于 EVI[24]的线性方程来估算FPARPAV，系数a=1.0。

植被总初级生产力GPPAPAR的计算式为

式中，LUE 为通量观测数据计算得到的光能利用率，APAR 由式（4）计算得到。

1.3.2 光化学反射植被指数−光能利用率模型（PRI 模型）

由于玉米的植被特征与草原、灌木和森林等植被特征差异较大，需要予以区别，因此选取单一均匀种植的玉米农田生态系统，进行区域反演。郭建茂等[25]在探究光化学植被指数与光能利用率之间的相关性时发现，在PRI488、PRI551、PRI667中，PRI551与LUE 的相关性最强，最适合反演LUE，并得到PRI-LUE 之间的线性关系。与APAR 相比，PRI 并不反映植被生物量、覆盖度和冠层结构等内在特征，PRI 与LUE 的内在生理机制决定了PRI 在各尺度上在监测植被LUE 方面的优势。

在Penuelas 等[26]对PRI 公式修正的基础上，可将PRI 定义为

式中，ρ531和ρλref分别表示测量波段和参照波段处的反射率。MODIS 第11 波段（526−536nm）的波段中心位于531nm，且该波段的宽度仅为10nm，符合计算PRI 的测量波段要求，因此将MODIS 第11波段的反射率作为计算 PRI 的测量波段。由于MODIS 数据缺少计算PRI 需要的理想参照波段，因此用MODIS 第10 波段（483−493nm）、第12 波段（546−556nm）和第13 波段（662−672nm）作为参照波段进行测试。虽然MODIS 第一波段（620−670nm）和第4 波段（545−565nm）的波谱也较接近，但是这两个波段的宽度都比计算PRI 要求的窄波段宽很多，其中MODIS 第1 波段的波谱宽度为50nm，第4 波段的波谱宽度为20nm。因此，将ρ531设为MODIS影像中第11 波段的反射率，第10（波段中心为488nm）、12（551nm）和13 波段（667nm）作为替代参照波段进行计算。

根据拟合出的PRI-LUE 之间的相关关系[25]（表1），分别计算上述3 个波段下的光能利用率LUE488、LUE551、LUE667。

式中，GPPλref为不同波段下的生态系统总生产力，单位为μmolC·m−2·d−1；λref 为参照波段处的反射率，APAR 由式（4）计算而得。

表1 PRI-LUE 各拟合线对应方程的参数Table 1 The parameters of the corresponding equations of every PRI-LUE fitting line

1.3.3 区域生产力模型（REG-PEM 模型）

REG-PEM 模型[9]是根据全球GPP 估算模型GLO-PEM 框架构建的一个由遥感数据驱动的区域陆地生态系统初级生产力模型，其中有关LUE 的计算式为

式中，LUEmax为潜在光能利用率，玉米的LUEmax为2.76g·MJ−1，即0.051μmolC·μmol−1photons；f(T)和f(W)分别为温度和水分对光能利用率的胁迫，即

式中，f(T)中植被进行光合作用的最低(Tmin)、最高（Tmax）和最适温度（Topt）分别定义为−1℃、50℃和22℃，当气温高于最高温度或低于最低温度时，f(T)为0[27]；f(W)中水分对植被光合作用的影响程度，用卫星遥感植被指数LSWI 计算。

根据短波红外波段对植被含水量和土壤湿度较敏感的特征，Xiao 等[28]构建了表示地表水分含量的陆表水分指数（Land Surface Water Index，LSWI），并将其定义为近红外波段和短波红外波段之间的归一化计算，即

式中，ρSWIR表示短红外波段的地表反射率，可利用MODIS 第6 波段（1628−1652nm）反射率代替。LSWI 植被指数用ENVI 软件中提取，利用IDL 程序提取站点周边3 ×3 像元的反射率均值来计算。

结合光能利用率模型，得到基于区域生产力模型估算作物生产潜力的计算式为

2 结果与分析

2.1 光合有效辐射吸收模型（APAR）的估算效果

由图1可见，在两个分析年度内（2013年和2014年），与锦州玉米农田生态系统野外观测站2013−2014年的实测GPP 数据相比，光合有效辐射吸收模型（即APAR 模型）估算的玉米生长季逐日GPP 值的变化特点基本一致，季节变化趋势基本相同，但总体上波动相对平缓。具体地，2013年GPP平均估算值为357.2μmolC·m−2·d−1，与GPP 观测值的平均相对误差为68.1%，生长季中期（第162−232天，以1月1日为1）的GPP 估算峰值低于观测值，但在生长季前期（第91−161 天）和后期（第233−303天）存在高估现象；2014年GPP 平均估算值为238.8μmolC·m−2·d−1，与GPP 观测值的平均相对误差为59%，同样存在峰值低估、谷值高估的现象。

图1 锦州生态系统观测站APAR 模型模拟的玉米生长季（4−10月）逐日GPP 值与实测值的比较Fig.1 Comparison of daily GPP estimation based on APAR and GPP observation at Jinzhou station

为探究APAR 在小时尺度上与GPP 的相关性，通过筛选、分析MODIS 影像过境时段内一小时APAR 与GPP 的关系（图2），发现GPP 与APAR 呈现经典的光饱和曲线，弱光条件或植被覆盖度较低时，GPP 随着APAR 的增大而显著提高，但当APAR增大到高值时（1000～1200μmol photons·m−2·s−1），GPP 趋近于饱和并有下降的趋势，此时APAR 对植被光合作用的响应能力减弱。

图2 观测站小时尺度MODIS 影像过境时段APAR 与实测值的关系Fig.2 The relationship between hourly APAR and GPP

2.2 光化学植被指数模型（PRI）的估算效果

从图3可以看出，基于模型估算的GPP488和GPP667与生态系统野外观测站的GPP 观测值呈正线性关系，而因为PRI551值为负数，表现为GPP551与GPP 观测值呈负线性关系，其中GPP551与GPP 观测值的相关性稍好于GPP488和GPP667，但GPP488和GPP667与观测值的相关性均未达到显著水平。

图3 观测站PRI 模型估算的小时尺度GPP 与实测值的比较Fig.3 Comparison of hourly GPP estimation based on PRI and GPP observation at Jinzhou

2.3 区域生产力模型（REG-PEM）的估算效果

图4所示，通过比较REG-PEM 模型估算的逐日GPP 值与涡度相关观测的GPP 值发现，GPP估算值与观测值在研究时间内的季节性规律一致，呈现由增到减的变化规律，峰值主要集中在第 200−250 天，较为符合玉米作物实际生长规律。具体地，2013年玉米生长季内GPP 观测值平均为357.2μmolC·m−2·d−1，估算值平均为415.9μmolC·m−2·d−1，二者平均相对误差为81.3%，GPP 估算值存在明显高估。在玉米生长季前期（第91−161 天）和后期（第233−303 天），GPP 估算值高于观测值，但在生长季中期（第162−232 天），GPP估算值和观测值同时到达峰值，二者较为接近，说明GPP 高估的现象主要发生在低值区；2014年玉米生长季内GPP 观测值平均为238.8μmolC·m−2·d−1，估算值平均为266.5μmolC·m−2·d−1，平均相对误差为54.2%，与2013年略有不同，GPP 估算值峰值低于实际，但在GPP 低值区仍然有明显高估。可能是以下原因：一方面是LUE 的误差，LUEmax给定值稍大，或气温和水分因子函数f(T)和f(W)均低估了温度和水分对实际植被LUE 的影响；另一方面，模型中存在FPAR 和LSWI 遥感参数的误差。

图4 观测站基于REG-PEM 模型的逐日GPP 估算值与实测值的比较Fig.4 Comparison of daily GPP estimation based on REG-PEM and GPP observation

2.4 三种方法的结果比较

通过检验皮尔逊相关系数来定量评价APAR 模型、PRI 模型和REG-PEM 模型在不同时间尺度上估算GPP 的效果，结果如表2。

由表2可知，逐日尺度上，REG-PEM 模型和APAR 模型都能较好地估算玉米农田生态系统的全天GPP。其中，锦州生态系统观测站2014年基于REG-PEM 估算的玉米农田生态系统总初级生产力GPP 与实测值的相关性（R2=0.71,P＜0.01）稍优于2013年（R2=0.58,P＜0.01）；利用APAR 模型估算GPP发现，2014年APAR 与GPP 间的相关性(R2=0.65，P＜0.01)稍好于2013年(R2=0.51，P＜0.01)，并且估算的GPP 与实际变化规律基本一致。小时尺度上，APAR 模型和PRI 模型对生态系统GPP 的估算效果明显减弱，且小时尺度上APAR 模型模拟值与GPP 观测值的相关关系要弱于逐日尺度，尤其在MODIS 影像集中的中午时段，植被叶片易出现由强光引起的光饱和现象或气温过高出现的午休现象，削弱了APAR对植被光合作用监测的敏感性。MODIS 影像提取的三种波段下的植被指数PRI估算的GPP与观测值的相关性由强到弱依次为，PRI551、PRI667、PRI488，这与MODIS PRIs 与LUE 的相关性强弱表现一致。

表2 基于三种模型估算GPP 与GPP 观测值的拟合方程Table 2 Estimation of fitting equations of GPP and GPP observations based on three models

3 结论与讨论

3.1 结论

以锦州站点周边为研究区，利用生态系统观测站观测的通量数据和卫星遥感植被指数（EVI、LSWI、PRI）相结合的方法，分别以APAR 模型、PRI 模型和REG-PEM 模型估算不同尺度下玉米生态系统的总初级生产力GPP，研究结果表明，逐日尺度上，REG-PEM 模型和APAR 模型都能较好地模拟农田GPP 基本变化规律，但二者普遍存在GPP 峰值低估、谷值高估的现象，但在不同年份高估或低估程度有所不同，其中2014年的表现好于2013年，APAR 模型的相对误差小于REG-PEM 模型。小时尺度上，APAR 模型和PRI 模型对GPP 的估算能力明显减弱。通过比较APAR 与GPP 的相关性发现，小时尺度上，其相关关系明显弱于逐日尺度上二者的相关关系，这主要是受到植被光合作用日变化规律的影响。提取MODIS 影像第10 波段（中心波段488nm）、第12 波段（中心波段551nm）和第13 波段（中心波段667nm）的PRI 值，分别计算光能利用率LUE，并估算GPP 值，PRI 与GPP 的相关性由强到弱依次为PRI551、PRI667、PRI488，这和MODIS PRIs 与LUE 的相关性强弱表现一致。

3.2 讨论

光能利用率模型结构简单，所需驱动参数少，由遥感参数直接获取，使用范围广，已成为陆地生态系统生产力模拟研究的重要发展方向，是研究农田生态系统生产力与碳循环的重要依据，但不同光能利用率模型的模拟存在很大的差异。REG-PEM 模型估算GPP 时给定了LUEmax，而植被LUE 受大气透射、植被冠层结构、非光合组织呼吸作用以及自身光化学转化效率等因素的影响，在同类型植被覆盖区域差异不明显，但在低植被覆盖度的混合作物种植区容易被高估[29]；气温和水分因子函数f(T)和f(W)可能低估了温度和水分对实际植被LUE 的影响。APAR 模型虽然建立在已有数据的基础上估算GPP，估算值与观测值的相关性从侧面说明，包含植被生物量、冠层结构及覆盖度信息的遥感因子APAR具有估算生态系统逐日GPP 的潜力；弱光或低植被覆盖度时，GPP 随APAR 的增大而增大，而中午时段，太阳辐射、气温将达到一天中的最高值，冠层叶片出现光饱和现象，气温升高，土壤水分蒸发，造成光合午休现象，植被自身光合作用效率降低[30−31]，从而削弱了APAR 对植被光合作用监测的敏感性，因此，表现出APAR 与GPP 的相关性在小时尺度上较差。与APAR 相比，PRI 通过对植被LUE 的影响机制，避免了因植被冠层结构特征变化不明显以及光合作用抑制现象造成的敏感性降低的现象。但不可否认的是，PRI 与GPP 的相关关系不够显著，PRI 在估算玉米生态系统GPP 方面的应用还需更深入的研究，本文认为利用PRI 估算GPP 有一定的参考价值。

由于数据的限制，本文只研究了玉米农田生态系统，随着更多生态系统通量观测站的建立，需要更全面的研究来验证遥感数据在估算农田生态系统光能利用率LUE 和总初级生产力GPP 方面的可行性。