技术领域
[0001] 本
发明属于生态学
生物生产力技术领域，具体涉及一种基于降尺度的高分辨率植被生产力遥感估算方法。
背景技术
[0002] 植被是陆地
生态系统的主体，植被生产力是
碳循环和碳收支研究的重要组成部分，反映了植被对大气中CO2的固定能力。总初级生产力(Gross Primary Productivity,GPP)是指生态系统中绿色
植物通过光合作用吸收
太阳能同化CO2制造的有机物，净初级生产力(Net Primary Production，NPP)是指GPP中再减去异养呼吸所消耗的光合产物。GPP和NPP不仅可以表征植被的生长状况及生长过程，还能直观地反映不同生态系统对全球变化的响应，同时也是判定生态系统碳源汇的重要手段。
[0003] 从遥感的
角度分析，尺度具体包括
时间分辨率、空间分辨率和
光谱分辨率。尺度转换是指将某一尺度上所获得的信息和知识扩展到其他尺度上的过程。通过尺度转换，可以使不同分辨率的数据具备相同的
时空分辨率，增强数据之间的可比性和相关性。遥感数据的尺度转换主要包括两个方面：升尺度和降尺度。升尺度是指从高分辨率数据向低分辨率数据转换；降尺度是指从低分辨率数据向高分辨率数据的转换。
[0004] 当前，不同分辨率的卫星
传感器，例如AVHRR、MODIS、TM/ETM+等，提供了大量
覆盖整个陆地生态系统的多尺度的数据源，从而形成了不同时空分辨率的影像数据层次体系。这些多尺度的海量遥感数据也被陆续应用到陆地植被生产力的估算以及植被生产力产品的生产中。当前高时间分辨率的遥感数据的空间分辨率较低，而高空间分辨率遥感数据由于受重访周期和
云的影响，限制了其在植被生产力产品的生产上的应用。在这种情况下，充分发挥不同分辨率遥感数据的优势，融合高、低分辨率的遥感数据，发展植被生产力产品降尺度方法，生产时间连续的高分辨率植被生产力数据产品，对于植被生产力的精确估算、研究植被生产力的变化状况及驱动力、评估区域碳收支与碳交换、监测及评估生态环境的变化趋势等，具有重要意义。
[0005] 植被生产力的估算经历了简单的统计模型、基于资源平衡理论的光能利用率模型、描述植物生理生态学特征的过程模型等多个发展阶段。遥感降尺度的方法大致可以分为：基于统计的降尺度、空间插值、
机器学习和基于光谱混合模型的降尺度等。基于统计模型的降尺度的基本假设是尺度因子与高、低分辨率影像之间的关系，不会随着影像分辨率的变化而变化。这种方法通过建立尺度因子与低分辨率影像的像元之间的统计关系，然后将该统计关系应用到高分辨率影像上，从而实现低分辨率影像的降尺度；空间插值是一种将离散的数据转换为连续的曲面的方法，即对一组已知的分区数据或者离散点，从这些已知数据中找出一个拟合函数，使该函数最好的逼近已知的空间数据，并能够根据该函数推导出区域范围内其他任意点或者任意区域的值；常用的机器学习方法，如
人工神经网络(Artificial Neural Network,ANN)、
支撑向量机(Support Vector Machine,SVM)、
决策树等机器学习技术，也逐渐应用于遥感降尺度；基于光谱混合模型的遥感降尺度的基本原理是根据线性光谱混合模型，关联高、低分辨率的影像，进而求解高分辨率的影像。
[0006] 现有的植被生产力降尺度方法仍较少，大多以植被指数为切入点：首先基于高、低分辨率的遥感数据和数据融合模型获得高分辨率的植被指数，如归一化植被指数NDVI、增强型植被指数EVI等；然后根据植被指数与光合有效
辐射吸收比(Fraction of Absorbed Photosynthetically Active Radiation，FPAR)之间的统计关系估算高分辨率FPAR；最后结合光能利用率模型、统计模型等，估算高分辨率的植被生产力。
[0007] 由于数据获取时间、成像视角等的差异、云和大气条件等的影响，FPAR和植被指数之间的关系会发生变化，所以，将植被指数作为融合切入点会导致FPAR估算结果的不确定性，从而给植被生产力的估算结果带来误差。所以，这种植被生产力降尺度方案具有一定的区域局限性。
发明内容
[0008] 本发明提供一种基于降尺度的高分辨率植被生产力遥感估算方法，充分发挥不同分辨率遥感数据的优势，融合高、低分辨率的遥感数据，构建一种
精度更高、普适性更强的高分辨率植被生产力产品的技术方案，解决现有的植被生产力降尺度方案精度不高、普适性不足的问题。
[0009] 具体的技术方案为：
[0010] 基于降尺度的高分辨率植被生产力遥感估算方法，包括以下步骤：
[0011] (1)将植被生产力估算的因子LAI/FPAR进行降尺度，得到时间序列的高分辨率的LAI/FPAR；
[0012] (2)对
温度进行空间插值和地形校正，得到高分辨率的温度因子；利用山地小
气候模型模拟高分辨率的太阳短波辐射；将降尺度后的高分辨率的LAI/FPAR、经过地形校正后高分辨率的温度因子和太阳短波辐射数据输入植被生产力模型MuSyQ-NPP，得到连续时间序列的高分辨率植被生产力。
[0013] 具体的，步骤(1)LAI/FPAR降尺度的具体过程包括：
[0014] 首先对低分辨率的MODIS LAI/FPAR产品进行数据
质量检测，选择数据质量高的像元；然后结合Landsat影像的地表反射率，利用变异系数的
阈值选择出地表异质性低的纯像元；最后，将聚合后的Landsat地表反射率与纯像元的LAI/FPAR作为样本，训练Cubist回归树模型，并将模型运用到高分辨率的Landsat地表反射率数据上，即可得到高分辨率Landsat分辨率的LAI/FPAR；
[0015] 将t1时刻低分辨率LAI/FPAR，高分辨率Landsat LAI/FPAR数据，以及t2时刻的低分辨率LAI/FPAR数据，输入遥感数据融合模型STARFM，即可得到降尺度后的高分辨率LAI/FPAR。
[0016] 其中，所述的纯像元的选择方法为：选择纯像元时，依据大像元内的Landsat像元的变异系数，计算一个波段内变异系数的均值，然后设定阈值来选择：
[0017]
[0018] 式中，CVij为变异系数，σij为MODIS大像元内的Landsat像元反射率的标准差，μij为MODIS大像元内的Landsat像元反射率的平均值；
[0019] 在选择纯像元的时候，除去第Lanssat TM数据第6波段以外，第1、2、3、4、5、7波段均参与纯像元的选择；
[0020] CVij描述了MODIS大像元内的Landsat小像元反射率的变化情况，值越小，该像元被认为越纯；
[0021] 然后计算这6个波段的CVij的均值CVmean。将CVmean的值排序，选择CVmean值最小的10％的MODIS像元，认为是纯像元。
[0022] 其中，步骤(2)高分辨率植被生产力估算，包括高分辨率温度、高分辨率辐射、高分辨率植被生产力的计算。具体的：
[0023] 高分辨率
温度计算，具体方法包括以下步骤：首先进行海平面校正，由垂直递减率和气象
站点海拔高度将站点温度校正为海平面的温度，即海拔高度为零的温度；然后采用克里金插值的方法进行插值，模拟海拔高度为零处的温度，并利用数字高程模型数据和气温垂直递减率，对温度进行地形校正，从而得到高分辨率的温度。
[0024] 高分辨率辐射计算，具体方法包括以下步骤：根据日最高气温、日最低气温等气象数据，利用山地小气候模型模拟高分辨率的太阳短波辐射。
[0025] 高分辨率(30m)GPP/NPP计算，具体方法包括以下步骤：采用MuSyQ-NPP植被生产力模型，GPP基于光能利用率原理来估算，自养呼吸通过生态过程方法来模拟，NPP为GPP和自养呼吸的差，具体公式如下：
[0026] GPP＝εg×FPAR×PAR×f1(T)×f2(β) (2)
[0027] NPP＝GPP-Ra (3)
[0028] 其中，εg为最大光能利用率，它的取值由植被覆盖类型决定；PAR为光合有效辐射；FPAR为植物可吸收的光合有效辐射分量，是植被冠层吸收的光合有效辐射占入射辐射的比例；f1(T)为温度胁迫因子；f2(β)为
水分胁迫因子；Ra为自养呼吸。
[0029] 其中自养呼吸的计算方法为：
[0030] Ra＝Rm+Rg＝∑(Rm，i+Rg，i) (17)
[0031] 式中，Ra为植被的自养呼吸，Rm为植被的维持呼吸，Rg为植被的生长呼吸，i代表植被的不同部分，i＝1、2、3，分别代表植物的叶、茎干和根；
[0032] 在计算Rm时，将土地覆盖类型分为林地和非林地；在非林地，Rm根据LAI和比叶面积计算：
[0033]
[0034] 在林地，Rm分为三部分，叶维持呼吸，茎干维持呼吸和根部维持呼吸；Rm是各部分的生物量与温度的函数：
[0035]
[0036] 其中，Mi是植被i部分参与呼吸的生物量，rm,i是植被i部分的维持呼吸系数，Q10,i是植被i部分的温度敏感因子，Tb为基温，T为日平均气温；
[0037] 假设植物各组织部分的生长呼吸系数γ相同，总生长呼吸的计算表示为：
[0038] Rg＝γ(GPP-Rm) (20)。
[0039] 本发明提供的基于降尺度的高分辨率植被生产力遥感估算方法精度较高，普适性较强。
附图说明
[0041] 图2为MuSyQ植被生产力模型流程图。
具体实施方式
[0042] 结合
实施例说明本发明的具体技术方案。
[0043] 本实施例选择在中国黑河流域验证区(38°10′～39°35′N，99°57′～101°46′E)。
[0044] 基于降尺度的高分辨率植被生产力遥感估算方法流程如图1。叶面积指数(Leaf Area Index,LAI)和FPAR是植被生产力模型MuSyQ-NPP的两个关键输入，所以，该植被生产力降尺度技术方案中，首先将植被生产力估算的因子LAI/FPAR进行降尺度，得到时间序列的高分辨率的LAI/FPAR；对温度进行空间插值和地形校正，得到高分辨率的温度因子；利用山地小气候模型(MT-CLIM)模拟高分辨率的太阳短波辐射；将降尺度后的高分辨率的LAI/FPAR、经过地形校正后高分辨率的温度因子和太阳短波辐射等数据输入植被生产力模型MuSyQ-NPP，得到连续时间序列的高分辨率(30m)GPP/NPP。具体的步骤为：
[0045] 1.LAI/FPAR降尺度。首先对低分辨率(1km)的MODIS LAI/FPAR产品进行数据质量检测，选择数据质量高的像元；然后结合Landsat影像的地表反射率，利用变异系数的阈值选择出地表异质性低的纯像元；最后，将聚合后的Landsat地表反射率与纯像元的LAI/FPAR作为样本，训练Cubist回归树模型，并将模型运用到高分辨率(30m)的Landsat地表反射率数据上，即可得到高分辨率Landsat分辨率的LAI/FPAR。将t1时刻低分辨率(1km)LAI/FPAR，高分辨率(30m)Landsat LAI/FPAR数据，以及t2时刻的低分辨率(1km)LAI/FPAR数据，输入遥感数据融合模型STARFM，即可得到降尺度后的高分辨率(30m)LAI/FPAR。
[0046] 本发明设计了纯像元的选择方法。选择纯像元时，依据大像元内的Landsat像元的变异系数，计算一个波段内变异系数的均值，然后设定阈值来选择：
[0047]
[0048] 式中，CVij为变异系数，σij为MODIS大像元内的Landsat像元反射率的标准差，μij为MODIS大像元内的Landsat像元反射率的平均值。在选择纯像元的时候，除去第Lanssat TM数据第6波段以外，第1、2、3、4、5、7波段均参与纯像元的选择。CVij描述了MODIS大像元内的Landsat小像元反射率的变化情况，值越小，该像元被认为越纯。然后计算这6个波段的CVij的均值CVmean。将CVmean的值排序，选择CVmean值最小的10％的MODIS像元，认为是纯像元。
[0049] 需要注意的是，在建立回归树模型的时候，除了Landsat TM数据的第1、6波段，第2、3、4、5、7波段均参与建立Cubist回归树模型，这是因为Landsat的第1波段经大气校正后
信噪比很低，第6波段为热红外波段，空间分辨率较低，且与其他波段分辨率不一致。
[0050] 2.高分辨率植被生产力估算。该部分分为高分辨率(30m)温度、高分辨率(30m)辐射、高分辨率(30m)GPP/NPP计算。具体方法如下：
[0051] (1)高分辨率(30m)温度计算：首先进行海平面校正，由垂直递减率和气象站点海拔高度将站点温度校正为海平面的温度，即海拔高度为零的温度。然后采用克里金插值的方法进行插值，模拟海拔高度为零处的温度，并利用数字高程模型(DEM)数据和气温垂直递减率，对温度进行地形校正，从而得到高分辨率(30m)的温度。
[0052] (2)高分辨率(30m)辐射计算：.根据日最高气温、日最低气温等气象数据，利用山地小气候模型(MT-CLIM)模拟高分辨率(30m)的太阳短波辐射。
[0053] (3)高分辨率(30m)GPP/NPP计算：MuSyQ-NPP植被生产力模型的流程如图2。该模型基于光能利用率和生态过程原理，考虑了植被自身的生物学特性和外界环境因子的共同影响。GPP基于光能利用率原理来估算，自养呼吸通过生态过程方法来模拟，NPP为GPP和自养呼吸的差。具体公式如下：
[0054] GPP＝εg×FPAR×PAR×f1(T)×f2(β) (2)
[0055] NPP＝GPP-Ra (3)
[0056] 其中，εg为最大光能利用率，它的取值由植被覆盖类型决定；PAR为光合有效辐射；FPAR为植物可吸收的光合有效辐射分量，是植被冠层吸收的光合有效辐射占入射辐射的比例；f1(T)为温度胁迫因子；f2(β)为
水分胁迫因子；Ra为自养呼吸。
[0057] ①光合有效辐射的计算
[0058] 通过转换系数由入射
太阳辐射得到光合有效辐射：
[0059] PAR＝0.50×R (4)
[0060] 式中，PAR为光合有效辐射，R为太阳短波辐射。
[0061] ②温度胁迫因子的计算
[0062] 温度是影响光能利用率的主要因素之一。温度胁迫因子的计算公式为：
[0063] f1(T)＝Ts1×Ts2 (5)
[0064] Ts1＝0.8+0.02Topt-0.0005(Topt)2 (6)
[0065]
[0066] 式中，Topt为最适宜生长温度，为区域年内植被长势最好时的月平均温度，取区域内LAI或NDVI最大时的月平均温度；T为日平均温度。Tε1反映了在高温和低温时，植物生化作用对植物光合作用产生抑制作用的强弱。Tε2表示周围
环境温度从最适宜温度Topt向高温或者低温变化时，植物的光能转化率逐渐变小的趋势。温度胁迫因子f1(T)的取值范围在0～1之间，温度越接近最适宜生长温度，f1(T)值越趋近于1，此时温度对植物生长的胁迫也越小；温度越偏离最适宜生长温度，f1(T)就越趋近于0，此时温度对植物生长的胁迫也就越大。
[0067] ③水分胁迫因子的计算
[0068] 水分胁迫因子采用实际蒸散E与潜在蒸散Ep进行计算。实际蒸散采用Penman-Monteith模型获得，潜在蒸散采用Priestley-Taylor方程计算。
[0069] f2(β)＝0.5+0.5(E/Ep) (8)
[0070] E＝Ecanopy+Esoil (9)
[0071]
[0072]
[0073]
[0074] Acanopy＝FPAR*A (13)
[0075] Asoil＝(1-FPAR)*A (14)
[0076] A≈Rns+Rnl (15)
[0077]
[0078] 其中，λ为
蒸发潜热系数，△为饱和水汽压曲线斜率，ρ为空气
密度，Cp为空气
比热容，VPD为水汽压差，ga为空
气动力学导度，γ为
干湿表常数，gsx为最大气孔导度，Qh为冠层顶部的可见光辐射，KQ、Q50、D50与植被类型有关，RH为空气
相对湿度，gtotc为修正后的总
空气动力学导度，为Prisestley-Taylor系数，Rns为短波净辐射，Rnl为长波净辐射。
[0079] ④自养呼吸的计算
[0080] 本发明设计了自养呼吸的计算方法。植物的自养呼吸可以分解为维持呼吸和生长呼吸。植被的维持呼吸与生物量和温度都有关系，而生长呼吸与温度无关，它只与总初级生产力成一定的比例关系。
[0081] Ra＝Rm+Rg＝∑(Rm，i+Rg，i) (17)
[0082] 式中，Ra为植被的自养呼吸，Rm为植被的维持呼吸，Rg为植被的生长呼吸，i代表植被的不同部分(i＝1、2、3，分别代表植物的叶、茎干和根)。在计算Rm时，将土地覆盖类型分为林地和非林地。在非林地，Rm根据LAI和比叶面积(Specific Leaf Area,SLA)计算：
[0083]
[0084] 在林地，Rm分为三部分，叶维持呼吸，茎干维持呼吸和根部维持呼吸。Rm是各部分的生物量与温度的函数：
[0085]
[0086] 其中，Mi是植被i部分参与呼吸的生物量，rm,i是植被i部分的维持呼吸系数，Q10,i是植被i部分的温度敏感因子，Tb为基温，T为日平均气温。
[0087] 假设植物各组织部分的生长呼吸系数(γ)相同，总生长呼吸的计算可表示为：
[0088] Rg＝γ(GPP-Rm) (20)
[0089] 本实施例在黑河流域验证区(38°10′～39°35′N，99°57′～101°46′E)，通过与地面观测数据的对比表明，该方案生产的植被生产力产品的平均相对误差(Mean Relative Error，MRE)小于10％。
