CASA模型全部流程

CASA模型估算NPP的全部流程

CASA模型中NPP的估算可以由植物的光合有效辐射（APAR）和实际光能利用率（ε）两个因子来表示，其估算公式如下：

NPP(x,t)APAR(x,t)(x,t)

式中，APAR(x,t)表示像元x在t月吸收的光合有效辐射 （g C·m-2·month-1），ε(x,t)表示像元x在t月的实际光能利用率（g C·MJ-1）

太阳总辐射Total radiation归一化植被指数NDVI土地覆盖类型Vegetation MapNDVI最大与最小值Maxumum and Minimumvalue of NDVI最大光能利用率εmax截栏的光合有效辐射IPAR光合有效辐射吸收比例FPAR光合有效辐射APAR温度、水分胁迫因子Temperature and waterscale实际光能利用率ε净初级生产力NPP净初级生产力（NPP）估算模型总体框架Frame of net primary productivity (NPP) Estimation Model

APAR的估算

APAR的值由植被所能吸收的太阳有效辐射和植被对入射光合有效辐射的吸收比例来确定。

APAR(x,t)SOL(x,t)FPAR(x,t)0.5

式中，SOL(x,t)表示t月在像元x处的太阳总辐射量（ g C·m-2·month-1 ），FPAR(x,t)植被层对入射光合有效辐射的吸收比例，常数0.5表示植被所能利用的太阳有效辐射（波长为0.4-0.7μm）占太阳总辐射的比例。

FPAR的估算

在一定范围内，FPAR与NDVI之间存在着线性关系（Ruimy & Saugier，1994），这一关系可以根据某一直被类型NDVI的最大值和最小值以及所对 应的FPAR最大值和最小值来确定。

(NDVI(x,t)NDVIi,min)(NDVIi,maxNDVIi,min)FPAR(x,t)(FPARmaxFPARmin)FPARmin

式中，NDVIi,max和NDVIi,min分别对应第i种植被类型的NDVI最大和最小值。

FPAR与比值植被指数（SR）也存在着较好的线性关系（Field et al.，1995;Los et al.，1994），可由以下公式表示：

(SR(x,t)SRi,min)(SRi,maxSRi,min)FPAR(x,t)(FPARmaxFPARmin)FPARmin

式中，FPARmin和FPARmax的取值与植被类型无关，分别为0.001和0.95； SRi,max和SRi,min分别对应第i种植被类型NDVI的95%和5%下侧百分位数， SR(x,t)由以下公式表示：

SR(x,t)1NDVI(x,t)1NDVI(x,t)

通过对FPAR-NDVI和FPAR-SR所估算结果的比较发现，由NDVI所估算的FPAR比实测值高，而由SR所估算的FPAR则低于实测值，但其误差小于直接由NDVI所估算的结果[1]，因此我们可以讲二者结合起来，取其加权平均或平均值作为估算FPAR的估算值：

FPAR(x,t)FPARNDVI(1)FPARSR

[1] 朱文泉, 潘耀忠, 张锦水. 中国陆地植被净初级生产力遥感估算, 植物生态学报 , 2007, 31(3), 413-424.

光能利用率的估算

光能利用率是在一定时期单位面积上生产的干物质中所包含的化学潜能与同一时间投射到该面积上的光合有效辐射能之比。环境因子如气温、土壤水分状况以及大气水汽压差等会通过影响植物的光合能力而调节植被的NPP。

(x,t)T1(x,t)T2(x,t)W(x,t)max

式中，Tε1(x, t)和 Tε2(x, t)表示低温和高温对光能利用率的胁迫作用； Wε(x,t)为水分胁迫影响系数，反映水分条件的影响； εmax是理想条件下的最大光能利用率（g C/MJ）。

植被分类图Vegetation calssification map归一化植被指数NDVI月平均气温Mean monthlytemperature月总降水量Total monthlyprecipitation最适温度The suitability oftemperature最大光能利用率Maximum light useefficiency温度胁迫系数1Temperature scale 1温度胁迫系数2Temperature scale 2水分胁迫系数Water scale光能利用率Light use efficiency光能利用率估算流程图Flow chart of light use of efficiency estimation

温度胁迫因子的估算

Tε1(x, t)的估算：其反映在低温和高温时植物内在的生化作用对光合的而降低第一性生产力。

Tε1(x,t)0.80.02Topt(x)0.0005[Topt(x)]2

式中，Topt(x)为植物生长的最适温度，定义为某一区域一年内NDVI值达到最高时的当月平均气温(℃)；当某一月平均温度小于或等于-10℃ 时，其值取0。

Tε2(x, t)的估算：表示环境温度从最适温度Topt(x)向高温或低温变化时植物光能利用率逐渐变小的趋势，这是因为低温和高温时高的呼吸消耗必将会降低光能利用率，生长在偏离最适温度的条件下，其光能利用率也一定会降低。

T2(x,t)1.184/{1exp[0.2(Topt(x)10T(x,t))]}1/{1exp[0.3(Topt(x)10T(x,t))]}

当某一月平均温度T(x, t)比最适温度Topt(x)高10℃或低13℃时，该月的Tε2(x, t)值等于月平均温度T(x, t)为最适温度Topt(x)时Tε2(x, t)值的一半。

水分胁迫因子的估算：

水分胁迫影响系数Wε(x,t)反映了植物所能利用的有效水分条件对光能利用率的影响,随着环境中有效水分的增加,Wε(x,t)逐渐增大,它的取值范围为0.5(在极端干旱条件下)到1(非常湿润条件下)。

W(x,t)0.50.5EET(x,t)/EPT(x,t)

式中：

EET为区域实际蒸散量（mm）； EPT为区域潜在蒸散量（mm）

最大光能利用率的确定

月最大光能利用率εmax的取值因不同的植被类型而有所不同，在CASA模型中全球植被的最大光能利用率为0.3 g C·MJ-1

MODIS NDVI数据TM高空间分辨率数据太阳总辐射Total radiationNDVI最大与最小值Maximum and MinimumValue of NDVI融合后的归一化植被指数 NDVI月平均气温Mean monthly temperature月总降水量Total monthly precipitation植被分类图Vegetation classification map截拦的光合有效辐射 IPAR光合有效辐射吸收比例 FPAR最适温度The suitability of temperature温度胁迫系数1Temperature scale1温度胁迫系数2Temperature scale2最大光能利用率Maximum light use efficiency水分胁迫系数Water scale光合有效辐射APAR实际光能利用率 εLight use efficiency净初级生产力（NPP）净初级生产力（NPP）估算模型总体框架Frame of primary productivity(NPP)Estimation Model

MODIS NDVI数据TM高空间分辨率数据太阳总辐射Total radiationNDVI最大与最小值Maximum and MinimumValue of NDVI融合后的归一化植被指数 NDVI月平均气温Mean monthly temperature月总降水量Total monthly precipitation植被分类图Vegetation classification map截拦的光合有效辐射 光合有效辐射吸收比例 IPARFPAR最适温度The suitability of temperature温度胁迫系数1Temperature scale1温度胁迫系数2Temperature scale2水分胁迫系数Water scale最大光能利用率Maximum light use efficiency光合有效辐射APAR实际光能利用率 εLight use efficiency净初级生产力（NPP）净初级生产力（NPP）估算模型总体框架Frame of primary productivity(NPP)Estimation Model

CASA模型具有以下优点:

1) 这是基于植被的生理过程而建立的植被NPP机理模型,并在大尺度植被NPP研究和全球碳循环研究中被广泛应用，是目前国际上最通用的NPP模型之一 ；

2) 考虑了NPP计算的两个主要驱动变量,即植被所吸收的光合有效辐射(APAR)与光能利用效率(ε) ，而这两个变量又分别通过太阳辐射、NDVI、土壤水分、降水量、平均温度等指标来体现；

3) CASA模型相对于其他模型所需要的输入参数较少,就避免了由于参数缺乏而人为简化或者估计而产生的误差;

4) 模型采用的遥感数据覆盖范围广,时间分辨率高,能够实现对区域和全球NPP的动态监测。

CASA模型的不足:

1) 该模型是针对北美地区所有植被而建立的， 世界各地差异较大，模型参数的修改比较困难；

2) 模型仅仅是在FPAR 的估算过程中，比值植被指数最大值SRmax 的确定时考虑了不同植被类型，但不能很好的从本质上揭示植被类型与NPP 的关系；

3) 光能利用率的准确估算是利用CASA 模型模拟生产力的关键因素之一，模型作者提出在理想状态下植被存在着最大光能利用率，不同植被类型的月值为0.3gC MJ-1。事实上，不同植被类型的光能利用率存在着很大差异，受到温度、水分、土壤、植物个体发育等因素的显著影响，把它作为一个常数在全球范围内使用会引起很大的误差；

4) 模型在估算水分胁迫因子时用到了土壤水分子模型，过程比较复杂，其中涉及到大量的参数，包括降水量、田间持水量、萎蔫含水量、土壤粘粒和砂粒的百分比、土壤深度、

土壤体积含水量等，数据较难获取，且通常土壤参数都是由土壤分类图来确定的，其精度难以保证。

基于Landsat、MODIS的数据同化算法研究

美国陆地资源卫星Landsat数据已经在土地利用/覆盖监测方面取得了极其广泛的应用。经过大气校正后的地表反射率产品是获取连续的地表变化参数的重要途径，然而Landsat的16天回访周期以及云污染极大的了Landsat影像在快速地表变化监测中的应用。结合Landsat的空间分辨率与MODIS的时间分辨率为解决这一问题提供了有效手段。MODIS影像具有半天至一天的回访周期，其空间分辨率有250m、500m、1000m三种级别，其各个波段的中心波长与Landsat影像非常接近，只是带宽要小一些。将两种影像相结合，生成具有Landsat空间分辨率和MODIS时间分辨率的预测影像，对于土地利用/覆盖短周期监测具有实际意义。

① STARFM融合模型理论基础

在忽略配准误差、大气校正误差的情况下，单一地物类型区域tk时刻的MODSI影像可以由同时间的Landsat影像加权平均得到，考虑到两者不同的光谱和时间分辨率，以及成像时的光照条件的差异，它们的反射率会有一定的偏差εk，假设(xi,yi)处的地表覆盖类型和各项系统误差不随时间变化，那么εk也就恒定，这样tn时刻的Landsat影像反射率就能由该时刻的MODIS影像以及tk时刻的Landsat、MODIS影像求得，公式如下：

L(xi,yi,tn)M(xi,yj,tn)L(xi,yj,tk)M(xi,yj,tk) （1）

式中L和M分别代表Landsat和MODIS影像的像素值。

实际上地表地物类型往往较复杂，相对于Landsat的30m分辨率来说，许多MODIS像元是混合像元，同时，地表覆盖类型及其状态（物候现象）和由光照条件引起的双向反射分布函数（BRDF）也会发生变化，因此仅利用该处单对像元的信息预测不出tn时刻的Landsat反射率。

若将邻近的具有相同光谱特征的像元作为辅助信息，则能极大的提高预测准确度。选取以预测像素为中心的邻域窗口，并利用权重函数W对窗口内像素进行卷积运算，确定中心像素的预测值，在整幅影像上滑动卷积窗口，得到预测影像。计算公式如下：

L(xw/2,yw/2,tn)Wijk(M(xi,yj,tn)L(xi,yj,tk)M(xi,yj,tk))i1j1ww （2）

式中w为滑动窗口的宽度，L(xw/2,yw/2,tn)为tn时刻滑动窗口中心像元的像素值。权重函数具有重要的作用，它决定了滑动窗口内各像素对预测值贡献的大小。算法依据光谱距离、时间距离、空间距离三项来确定权重函数。

光谱距离定义为tk时刻给定位置处的Landsat像素值与MODIS像素值之差：

SijL(xi,yj,tk)M(xi,yj,tk) （3）

时间距离定义为给定位置处的tn时刻MODIS与tk时刻MODIS像素值之差：

TijM(xi,yj,tn)M(xi,yj,tk) （4）

空间距离定义为窗口内中心像元与邻域像元的欧式距离：

dij(xw/2xi)2(yw/2yj)2 （5）

②模型实现

算法的具体实现需要考虑的更多的问题，如与中心像元具有相似光谱特征的邻域像元选取，权重函数的实现形式，以及异常候选像元剔除等。

 候选像元选取

利用光谱特征相似的邻域像元才能确保选择正确的辅助信息。有两种获取相似光谱特征的邻域像元的方法。一种是对整幅Landsat影像做非监督分类，认为相同类别的地物具有相似的光谱特征，这种方法具有全局性，整幅影像用的是同一个标准。另一种方法是阈值法，将与中心像元的光谱之差小于阈值的邻域像元作为候选像元，具体实现过程为：计算各波段窗口内像素值的标准差，所有波段都满足式（6）的像元选为候选像元。

V(xi,yi)V(xw/2,yw/2)2(B)/n （6）

式中V表示像元值，σ(B)为B波段窗口内图像的标准差，n为窗口内地物类别数。阈值法具有局部适应性，选择标准随着滑动窗口的位置变化而不同，局部相似像元的选取失败并不会影像整幅影像预测的准确性，本文利用该方法选取光谱相似像元。

 异常像元剔除

基于相似光谱特征选择候选邻域像元后，还需要提出异常值。首先根据Landsat和MODIS反射率产品的质量评价数据掩膜，剔除质量差的候选像元；其次，候选像元提供

的辅助信息量如果低于中心像元，则剔除。合格的候选像元应符合下式：

Sijmax(L(xw/2,yw/2,tk,B)M(xw/2,yw/2,tk,B))Tijmax(M(xw/2,yw/2,tk,B)M(xw/2,yw/2,tn,B)) （7）

由于影像处理过程不同、大气校正误差、配准误差等，Landsat和MODIS反射率的值会存在一定的误差，假设分别为ξ1和ξm，所有放射率值都相互，则二者之间的光谱距离误差为：

2lml2m （8）

不同时刻的MODIS影像的时间距离误差为：

22mmmm2m （9）

将（8）和（9）两式代入（7）中有：

Sijmax(L(xw/2,yw/2,tk,B)M(xw/2,yw/2,tk,B))lmTijmax(M(xw/2,yw/2,tk,B)M(xw/2,yw/2,tn,B))mm （10）

 权重函数实现

在计算滑动窗口的权重函数之前，先要利用公式（11）将空间距离转换为与中心像元的相对距离。

Dij1dij/(w/2) （11）

从上式可以看出相对距离Dij的变动范围是[1,1+√2]。

在确定权重函数时，有如下假设：1）同一地物类别区域的MODIS影像反射率随时间变化的值与Landsat相同；2）在从基准时间到预测时间这段时期内，变化越小的像元提供的辅助信息越可靠；3）距离中心像元越近的像元提供的辅助信息越可靠。基于以上三种假设，定义组合距离

CijSijTijDij （12）

组合距离越大的像元权重越小，因此定义函数权重为组合距离的倒数，归一化后得权重函数表达式：

Wij(1/Cij)/(1/Cij)i1j1ww （13）

③技术流程

图3为STARFM模型的技术流程。需要注意的是首先要将MODIS影像投影到与Landsat相同的坐标系下，即WGS84椭球，UTM投影。

图3 STARFM模型技术流程图