CN106446564A - 一种植被净初级生产力遥感估算方法 - Google Patents

Abstract

本发明公开了一种植被净初级生产力遥感估算方法，包括步骤：对遥感数据进行数据预处理；提取遥感数据中MODIS产品的LAI数据，计算FPAR；根据研究区的地理位置信息，包括经纬度和日期，计算天文辐射总量Q₀；根据研究区域过往的总辐射和日照百分率数据，提取经验系数；使用Q₀和经验系数计算太阳总辐射量SOL；根据FPAR和SOL计算光合有效辐射APAR；提取遥感数据中的土壤含水量和经验参数，计算当前的叶水势，并结合以往的经验系数计算水分胁迫因子提取气象数据中研究区的温度数据，计算温度胁迫因子f₁(T_a)；根据f₁(T_a)、和最大光能利用率，计算光能利用率ε；根据APAR和ε，计算净初级生产力NPP。本发明具有估算精度高的优点。

Description

一种植被净初级生产力遥感估算方法

技术领域

本发明涉及遥感分析、农业遥感应用研究领域，特别涉及一种植被净初级生产力遥感估算方法。

背景技术

植被净初级生产力(NPP，Net Primary Production)是绿色植物单位时间、单位面积所积累的有机物数量，是生态系统的基本数量特征之一，其大小直接反应生态系统自身的生物学特征和系统结构特征，以及环境因子对森林生长的综合影响。准确估算植被净初级生产力，对全球碳循环、陆地生态系统演化和生态系统碳汇以及生态调节过程的判定有着重要的作用及意义。

近年来，随着遥感技术的发展，遥感技术逐渐被用于对植被初级生产力的估算，遥感在植被净初级生产力上的应用可以总结为三类模型：1、统计模型，如1975年Leith提出的Miami模型，1985年Uchijima提出的Chikugo模型等，该模型是通过测量样本区的真实NPP值，从遥感影像上获取NDVI/LAI值，建立NPP与NDVI/LAI间得相关关系，此方法属于经验模型，具有地域性特点，难以应用于别的区域；2、过程模型，如1987年Parton提出的Century模型，1996年Hunt提出的BIOME-BGC模型等，该模型获取的参数都是从植被生理角度出发，精度高，过程复杂，许多参数无法从遥感学角度获取，难以推广应用；3、光能利用率模型，如1993年Potter提出的CASA模型等，此模型是目前应用最多的模型，但是其水分胁迫系数和FPAR属于经验值，误差大、难以推广。

现有技术中，周广胜、郑元润等利用卫星影像提取NDVI或LAI，建立与地面NPP实测值之间的回归分析模型对中国陆地植被的NPP进行了估算。这种方法只考虑了反映植被自身情况的一个参数(NDVI或LAI)，没有考虑气候、土壤等外在条件的影响，仅从遥感数据出发，造成估算的误差较大。而周才平、江洪等从植物的生态机理过程模型角度出发，利用MODIS遥感影像，采用BEPS过程模型对福建森林生态系统的NPP进行了评价。但是过程模型从植物的生态机理出发，结构复杂，所需的参数较多，难以得到推广。朴世龙、陈正华、董丹等根据不同的卫星遥感数据，采用CASA模型，估算了中国不同区域森林净初级生产力，验证了CASA模型的输入参数完全可以通过遥感手段直接获得，无需野外测定。但是CASA模型是针对于北美地区建立的，在其他地区的应用中，其模型参数是否有效及如何修正是一个难题。

为此，寻求一种通过遥感手段能够有效获取植被净初级生产力的估算方法具有重要研究意义。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种植被净初级生产力遥感估算方法，该方法在对植被光合有效辐射的吸收比例(FPAR)和水分胁迫系数进行修正的基础上，建立适用于研究区的植被净初级生产力遥感估算模型，模拟研究区植被净初级生产力，解决野外实测NPP的困难，具有估算精度高的优点。

本发明的目的通过以下的技术方案实现：一种植被净初级生产力遥感估算方法，包括步骤：

(1)对遥感数据进行数据预处理；

(2)提取遥感数据中MODIS产品的LAI数据，计算植被光合有效辐射吸收比例参数FPAR；

(3)根据研究区的地理位置信息，包括经纬度和日期，计算天文辐射总量Q₀；根据研究区域过往的总辐射和日照百分率数据，提取经验系数；使用天文辐射总量Q₀和经验系数计算太阳总辐射量SOL；

(4)根据步骤(2)得到的FPAR和步骤(3)得到的太阳总辐射量SOL计算光合有效辐射APAR；

(5)提取遥感数据中的土壤含水量和经验参数，计算当前的叶水势，并结合以往的经验系数计算水分胁迫因子

(6)提取气象数据中研究区的温度数据，计算温度胁迫因子f₁(T_a)；

(7)根据温度胁迫因子f₁(T_a)、水分胁迫因子和最大光能利用率，计算光能利用率ε；

(8)根据光合有效辐射APAR和光能利用率ε，计算净初级生产力NPP；

NPP(x,t)＝APAR(x,t)×ε(x,t)；

式中，x表示空间位置，t表示时间。

优选的，所述遥感数据包括“高分一号”卫星影像数据、Landsat7卫星影像数据和MODIS产品MOD13A2级NDVI植被指数数据、MOD15A2级FPAR、LAI数据。

优选的，步骤(1)中，对遥感数据进行数据预处理，包括辐射定标、大气校正以及几何精校正；

所述辐射定标计算公式如下：

L＝Gain×DN+Bias；

式中，L为辐亮度，单位为W·m^-2·sr^-1·um^-1，Gain为增益，单位为W·m^-2·sr^-1·um^-1，DN为卫星载荷观测值，Bias为偏置，单位为W·m^-2·sr^-1·um^-1；

大气校正采用ENVI软件下的FLAASH模块；

几何精校正采用图像对地图的方式进行，遥感数据中原始的MODIS数据为HDF格式，正弦投影(Sinusoidal projection)，将其进行重投影、裁剪研究区等预处理，得到Geographic-WGS84投影的IMGAE文件。

优选的，步骤(2)中，根据研究区LAI数据，降尺度处理后做基于像元的相关性分析，得到FPAR-LAI的相关性，根据相关性得到FPAR的计算公式如下：

FPAR＝A·LAI^B；

其中，A、B均为通过拟合得到的参数。

优选的，步骤(3)中，天文辐射总量的计算公式如下：

式中，Q₀单位为MJ·m^-2·d^-1，T为周期，I₀为太阳常数13.67×10^-4MJ·m^-2·s^-1，ρ为日地相对距离，ω_θ为日落时角(rad)，为地理纬度，δ为太阳赤纬(rad)；

日地相对距离ρ的计算公式为：

日落时角ω_θ的计算公式为：

太阳赤纬δ的计算公式为：δ＝0.409sin(0.0172J-1.39)，J为年内的天数，从1月1日的0到12月31日的364；

太阳总辐射量SOL的计算公式如下：

SOL(x,t)＝Q₀(a+bs)；

式中，a,b为经验系数，s为日照百分率，s＝n/N，n为实际日照实数，N为最大日照实数，N＝24/π×ω_θ。

更进一步的，根据研究区历年各月的总辐射和日照百分率数据，做最小二乘法拟合，拟合公式中各月的经验系数即为经验系数a,b。

优选的，步骤(4)中，光合有效辐射APAR的计算公式如下：

APAR(x,t)＝FPAR(x,t)×SOL(x,t)×0.5；

其中，0.5表示植被所能利用的太阳有效辐射(波长为0.4-0.7um)占太阳总辐射的比例。

优选的，步骤(5)中，水分胁迫因子的计算公式如下：

式中，为当前的叶水势，b为常数，取值0.009，w_s为土壤水含量，为气孔半闭时的叶水势,取经验值-0.63Mpa，n为形状经验参数，取值5，r_{s min}取经验值1，r_a相对较小，可忽略不计，r_m取值4.5。

优选的，步骤(6)中，温度胁迫因子f₁(T_a)的计算公式如下：

式中，T_a为气温，T_max、T_min、T_opt分别为最高温度、最低温度和最适宜温度，最适宜温度采用生长旺季的平均温度。

优选的，步骤(7)中，光能利用率ε的计算公式如下：

式中，ε^*为理想条件下的最大光能利用率，取常值0.389g.MJ^-1。

本发明与现有技术相比，具有如下优点和有益效果：

本发明以“高分一号”卫星可见光、近红外和TM热红外为数据源，借鉴CASA模型的基本架构，在对植被光合有效辐射的吸收比例(FPAR)和水分胁迫系数进行修正的基础上，建立适用于研究区的植被净初级生产力遥感估算模型，模拟研究区植被净初级生产力，解决野外实测NPP的困难，具有估算精度高的优点。

附图说明

图1是本实施例方法的流程示意图。

图2是本实施例方法中各个参数的处理流程图。

图3(a)、(b)是FPAR与NDVI的拟合关系示意图。

图4(a)、(b)是FPAR与LAI的拟合关系示意图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1

为了提高植被净初级生产力的估算精度，本发明在前人研究的基础上，对光能利用率模型(CASE)中的植被光合有效辐射吸收比例(FPAR)和水分胁迫系数进行了修正，以便建立适用于不同研究区的植被净初级生产力遥感估算模型。方法流程参见图1、2。

本实施例以广州市增城区为研究区，首次运用2013年10月的高空间分辨率的国产“高分一号”可见光-近红外数据产品和TM热红外数据反演修正后的CASE模型参数，估算陆地植被净初级生产力，是一种全新的研究手段，下面根据该实例对本实施例植被净初级生产力遥感估算方法的步骤进行具体说明，并将其模拟结果与其他模型模拟的NPP进行比较。

一、数据提取

本实施例所使用的遥感数据包括“高分一号”卫星影像数据、Landsat7卫星影像数据和MODIS产品MOD13A2级NDVI植被指数数据、MOD15A2级FPAR、LAI数据。“高分一号”卫星影像数据获取时间为2013年10月1日，空间分辨率为8m。Landsat7卫星影像数据获取时间为2013年10月4日，空间分辨率为30m。MOD13A2级、MOD15A2级数据获取时间为2013年10月8日，空间分辨率为1KM。

气象数据来源于中国气象科学数据共享服务网，时间为2013年10月，数据内容包括日照实数、降雨量、气温以及各站点的经度、纬度和海拔高度。根据经纬度挑选出位于研究区的站点实测值。

植被净初级生产力实测数据采用野外调查的方法获得，调查主要针对研究区森林植被，由于缺少影像同步的野外观测数据，考虑到森林植被净初级生产力年季节变化不大，将调查时间调整为2014年10月。由于地形等外在因素影响，本实施例选择了10个样方，考虑到高分卫星影像的空间分辨率为8m，每个样方设置为8m*8m的大小。用手持式GPS进行样方定位，记录样方的中心点坐标，测出样方内数木的树高H，胸径D，根据曾伟生在全国立木生物量方程建模方法研究中建立的生物量方程计算出生物量，再除以样方内树木的平均生长周期得到净初级生产力实测值。

二、数据预处理

遥感数据预处理主要包括辐射定标、大气校正和几何精纠正等。辐射定标即将遥感卫星影像各载荷的通道观测数字量化输出值DN值转换为卫星载荷入瞳处等效表观辐亮度数据，计算式如下：

L＝Gain×DN+Bias(1)

式中，L为辐亮度(W·m^-2·sr^-1·um^-1)，Gain为增益(W·m^-2·sr^-1·um^-1)，DN为卫星载荷观测值，Bias为偏置(W·m^-2·sr^-1·um^-1)。辐射定标系数取值如表1所示。本实施例利用中国资源卫星应用中心提供的“高分一号”卫星影像定标系数和NASA官网提供的Landsat7卫星影像定标系数将影像进行辐射定标处理。

表1高分一号及Landsat TM7影像辐射定标系数

为消除大气对影像的影响，获得真实的地表反射率，需对影像进行大气校正。本实施例采用ENVI软件下的FLAASH模块进行大气校正。为了消除投影差和几何畸变，采用图像对地图的方式进行几何精纠正。原始的MODIS数据为HDF格式，正弦投影，将其进行重投影、裁剪研究区等预处理，得到Geographic-WGS84投影的IMGAE文件。

三、构建模型

本实施例所使用的NPP估算模型是由Monteith(1977)提出的CASA模型，该模型从资源平衡的观点出发，考虑水分、温度和养分的胁迫建立起来的。在模型中，植被的净初级生产力由植被所吸收的光合有效辐射(APAR)与光能转化率(ε)两个变量来确定，公式如下：

NPP(x,t)＝APAR(x,t)×ε(x,t)(2)

式中，x表示空间位置，t表示时间，APAR表示植被吸收的光合有效辐射，ε表示光能利用率。结合图2，下面对上述公式中的各个参数进行计算。

1、植被吸收的光合有效辐射APAR参数计算

植被吸收的光合有效辐射APAR是植物光合作用的驱动力，表示为：

APAR(x,t)＝FPAR(x,t)×SOL(x,t)×0.5 (3)

上式中，0.5表示植被所能利用的太阳有效辐射(波长为0.4-0.7um)占太阳总辐射的比例，SOL为太阳总辐射量，FPAR为植被光合有效辐射的吸收比例，它是获得植被生产力的关键变量。太阳总辐射量SOL的计算公式如下：

SOL(x,t)＝Q₀(a+bs) (4)

式(4)中，a,b为经验系数，s为日照百分率，Q₀为天文总辐射。经验系数a,b采用杜尧东等(2003)根据广州、汕头的历年各月的总辐射和日照百分率数据，用最小二乘法拟合出的公式中各月的经验系数(如表2所示)。

表2广州各月经验系数和统计检验

Table 1 Experiential coefficients and statistical tests for differentmonths in Guangzhou

式(4)中天文辐射总量的计算公式如下：

式中，Q₀为天文辐射总量，单位为MJ·m^-2·d^-1，T为周期(24*60*60s)，I₀为太阳常数(13.67×10^-4MJ·m^-2·s^-1)，ρ为日地相对距离，ω_θ为日落时角(rad)，为地理纬度，δ为太阳赤纬(rad)。日地相对距离的计算如式(6)，太阳赤纬计算如式(7)，日落时角计算如式(8)，日照百分率计算如式(9)：

δ＝0.409sin(0.0172J-1.39)(7)

s＝n/N (9)

上式中，J为年内的天数，从1月1日的0到12月31日的364，n为实际日照实数，N为最大日照实数，计算公式如下：

N＝24/π×ω_θ (10)

大多数的研究表明，FPAR与归一化植被指数(NDVI)间存在着线性关系，与叶面积指数(LAI)间存在着指数的关系。这些通过分析FPAR与NDVI/LAI间的关系建立起来的FPAR估算模型具有一定的地域差异和时间差异，因此，本实施例利用研究区2013年10月MODIS1KM分辨率的FPAR、NDVI、LAI产品数据，降尺度处理后做基于像元的相关性分析，图3表示本实施例FPAR与NDVI的两种拟合关系图，图3(a)中，图中实线Linear表示线性分析曲线，该拟合直线公式为y＝0.8126x+0.272，图3(b)中，实线Power表示拟合的曲线，该曲线的公式为y＝0.9379x^0.512。图4表示本实施例FPAR与LAI的两种拟合关系图，图4(a)中，图中实线Linear表示线性分析曲线，该拟合直线公式为y＝3.2469x-0.4463，图4(b)中，实线Power表示拟合的曲线，该曲线的公式为y＝0.4503x^0.6169。通过上述计算，可以发现FPAR-LAI中曲线拟合的相关性相较于FPAR-NDVI更显著，R²达到0.919，因此式(2)中FPAR的计算公式如下：

FPAR＝0.4503LAI^0.6169 (10)

2、光能利用率ε参数计算

生长季节内光能利用率ε是随环境(土壤湿度、大气水汽压和叶子水势等)和气候变化的量，在没有任何因素影响的条件下植被光能利用率能达到最大，实际条件下，光能利用率ε主要受温度和水分胁迫因子的影响，计算公式如下：

式中，ε^*为理想条件下的最大光能利用率，取常值0.389g.MJ^-1，f₁(T_a)是温度胁迫因子，是水分胁迫因子。其中，温度胁迫因子采用Cao MK等(2004)提出的方法，公式如下：

式中，T_a为气温，T_max、T_min、T_opt分别为最高温度、最低温度和最适宜温度，最适宜温度可以通过生长旺季的平均温度来替代。

公式(11)中，鉴于以往水分胁迫因子多采用建立与区域实际蒸散量和潜在蒸散量的相关关系来计算，如周广胜-张新时模型。但大量的研究表明，土壤水分的增减是通过植物叶水势来影响光合作用的，因此，本实施例利用植物的叶水势来计算水分胁迫因子，公式如下：

式中，为当前的叶水势，b为常数，取值0.009，w_s为土壤水含量，为气孔半闭时的叶水势,取经验值-0.63Mpa，n为形状经验参数，取值5，r_{s min}取经验值1，r_a相对较小，可忽略不计，r_m取值4.5。具体方法可参考刘岩著的《半干旱草地净第一性生产力遥感模型研究》。

植被所吸收的光合有效辐射跟植被本身有很大的关系，本实施例根据公式，计算出植被光合有效辐射的吸收比例FPAR中，林地的FPAR值最高，达1.892，表明林地的光和有效辐射吸收能力最强，草地次之，建筑的FPAR值最小，接近0。理论上，建筑的FPAR值应为0，本实施例模拟出的FPAR值最小为0.01，分析原因，建筑物周边有绿化，对其有一定的影响。

考虑到LandsatTM7卫星影像的热红外数据是30m分辨率，先降尺度处理到8m分辨率后，再分别使用单窗算法和基于地表温度-植被指数特征空间的方法来反演地表温度和土壤水分，进而计算水分胁迫系数，根据水分胁迫系数的计算结果，高值区为林地和草地，其水分胁迫系数均大于1.21，低值区为建筑区，不透水率高，土壤水分低，水分胁迫系数低。

根据公式，计算出研究区10月的NPP值，同时为了验证本实施例模型的可靠性及合理性，用CASA模型模拟出研究区的NPP，通过比较二者结果可发现，本实施例模拟出的NPP值要大于未改进模拟的值。分析原因可知：建筑用地不透水面比例高，绿地面积相对较低，是净初级生产力的低值区域，在221.626gc/m².a以下；林地的净初级生产力最高，在221.626gc/m².a-660.593gc/m².a之间，这跟山区受人类活动干扰较小以及植被的浓密程度有关，植被净初级生产力随着海拔的增高而增高，即从建筑区到山区，净初级生产力逐渐增大。植被净初级生产力是绿色植被呼吸后所剩下的单位面积时间内所固定的能量或生产的有机物质，因此，NPP值越高的地方，代表森林有机物质越高。

将本实施例计算结果与常规的CASE模拟结果比较，精度得到很大的提高，随机抽取10个样点的结果作比较，9个样点的误差都减少了，改进的模型比常规模型误差减少了4％---15％，结果如表3所示。

表3模拟值与实测值得比较(gc/m².a)

注：误差的计算方法为：

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种植被净初级生产力遥感估算方法，其特征在于，包括步骤：

(1)对遥感数据进行数据预处理；

(2)提取遥感数据中MODIS产品的LAI数据，计算植被光合有效辐射吸收比例参数FPAR；

(3)根据研究区的地理位置信息，包括经纬度和日期，计算天文辐射总量Q₀；根据研究区域过往的总辐射和日照百分率数据，提取经验系数；使用天文辐射总量Q₀和经验系数计算太阳总辐射量SOL；

(4)根据步骤(2)得到的FPAR和步骤(3)得到的太阳总辐射量SOL计算光合有效辐射APAR；

(5)提取遥感数据中的土壤含水量和经验参数，计算当前的叶水势，并结合以往的经验系数计算水分胁迫因子

(6)提取气象数据中研究区的温度数据，计算温度胁迫因子f₁(T_a)；

(7)根据温度胁迫因子f₁(T_a)、水分胁迫因子和最大光能利用率，计算光能利用率ε；

(8)根据光合有效辐射APAR和光能利用率ε，计算净初级生产力NPP；

NPP(x,t)＝APAR(x,t)×ε(x,t)；

式中，x表示空间位置，t表示时间。

2.根据权利要求1所述的植被净初级生产力遥感估算方法，其特征在于，所述遥感数据包括“高分一号”卫星影像数据、Landsat7卫星影像数据和MODIS产品MOD13A2级NDVI植被指数数据、MOD15A2级FPAR、LAI数据。

3.根据权利要求1所述的植被净初级生产力遥感估算方法，其特征在于，步骤(1)中，对遥感数据进行数据预处理，包括辐射定标、大气校正以及几何精校正；

所述辐射定标计算公式如下：

L＝Gain×DN+Bias；

式中，L为辐亮度，单位为W·m^-2·sr^-1·um^-1，Gain为增益，单位为W·m^-2·sr^-1·um^-1，DN为卫星载荷观测值，Bias为偏置，单位为W·m^-2·sr^-1·um^-1；

大气校正采用ENVI软件下的FLAASH模块；

几何精校正采用图像对地图的方式进行，遥感数据中原始的MODIS数据为HDF格式，正弦投影，将其进行重投影、裁剪研究区，得到Geographic-WGS84投影的IMGAE文件。

4.根据权利要求1所述的植被净初级生产力遥感估算方法，其特征在于，步骤(2)中，根据研究区LAI数据，降尺度处理后做基于像元的相关性分析，得到FPAR-LAI的相关性，根据相关性得到FPAR的计算公式如下：

FPAR＝A·LAI^B；

其中，A、B均为通过拟合得到的参数。

5.根据权利要求1所述的植被净初级生产力遥感估算方法，其特征在于，步骤(3)中，天文辐射总量的计算公式如下：

式中，Q₀单位为MJ·m^-2·d^-1，T为周期，I₀为太阳常数13.67×10^-4MJ·m^-2·s^-1，ρ为日地相对距离，ω_θ为日落时角，为地理纬度，δ为太阳赤纬；

日地相对距离ρ的计算公式为：

日落时角ω_θ的计算公式为：

太阳赤纬δ的计算公式为：δ＝0.409sin(0.0172J-1.39)，J为年内的天数，从1月1日的0到12月31日的364；

太阳总辐射量SOL的计算公式如下：

SOL(x,t)＝Q₀(a+bs)；

式中，a,b为经验系数，s为日照百分率，s＝n/N，n为实际日照实数，N为最大日照实数，N＝24/π×ω_θ。

6.根据权利要求5所述的植被净初级生产力遥感估算方法，其特征在于，根据研究区历年各月的总辐射和日照百分率数据，做最小二乘法拟合，拟合公式中各月的经验系数即为经验系数a,b。

7.根据权利要求5所述的植被净初级生产力遥感估算方法，其特征在于，步骤(4)中，光合有效辐射APAR的计算公式如下：

APAR(x,t)＝FPAR(x,t)×SOL(x,t)×0.5；

其中，0.5表示植被所能利用的太阳有效辐射占太阳总辐射的比例。

8.根据权利要求1所述的植被净初级生产力遥感估算方法，其特征在于，步骤(5)中，水分胁迫因子的计算公式如下：

式中，为当前的叶水势，b为常数，取值0.009，w_s为土壤水含量，为气孔半闭时的叶水势,取经验值-0.63Mpa，n为形状经验参数，取值5，r_smin取经验值1，r_a相对较小，可忽略不计，r_m取值4.5。

9.根据权利要求8所述的植被净初级生产力遥感估算方法，其特征在于，步骤(6)中，温度胁迫因子f₁(T_a)的计算公式如下：

$f_1\left(T_a\right)=\frac{(T_{min}-T_a)(T_{max}-T_a)}{(T_{min}-T_a)(T_{max}-T_a)-{(T_{opt}-T_a)}^2};$

式中，T_a为气温，T_max、T_min、T_opt分别为最高温度、最低温度和最适宜温度，最适宜温度采用生长旺季的平均温度。

10.根据权利要求9所述的植被净初级生产力遥感估算方法，其特征在于，步骤(7)中，光能利用率ε的计算公式如下：

式中，ε^*为理想条件下的最大光能利用率，取常值0.389g.MJ^-1。