CN104166782A - 一种林地土壤有机碳遥感估测的方法 - Google Patents

Abstract

本发明提供了一种林地土壤有机碳遥感估测的方法，由于林地植被覆盖区遥感图像像元则表现为植被的光谱特征，遥感影像植被指数NDVI是目前应用最为广泛和适用的一个植被指数，广泛应用于植被性第一性生产力的估测研究，为此，在本发明中，通过CASA生态过程模型，在应用NDVI遥感数据(利用高分辨遥感影像反映区域差异性)，解决传统土壤有机碳估测模型中，一些具有空间异质性的模型参数被简化为常数，而降低了估测精度的不足，以及土壤有机碳储量有密切关系的土壤基础呼吸系数，建立土壤有机碳储量与土壤基础呼吸的关系模型，实现林地土壤有机碳的遥感估测研究。

Description

一种林地土壤有机碳遥感估测的方法

技术领域

本发明属于林地土壤估测技术，尤其涉及一种林地土壤有机碳遥感估测的方法。

背景技术

自上世纪50年代研究土壤有机碳以来，土壤有机碳估测方法在不断的改进、提升与完善，主要有植被类型法、土壤类型法、生命带法、模型法等4种类型。其中，土壤类型法要求具备完备、较为完整的土壤理化性质数据(包括土层厚度等)、分布图，在实际操作中难以做到，而且受不同国家土壤分类系统的不统一等因素影响，类型法的估测差异性较大。在人类对社会开发显多元化发展和急剧式变化的影响下，基于植被与土壤类型一一对应的理论假设为基础的生命带法与植被类型法的估测结果带来不确定因素，计算误差也会较大。模型法利用剖面数据实现相似土壤、生态区域土壤碳储量的外推与计算，已发展有相关关系模型、机理过程模型，也有基于实测样本和遥感数据耦合的数学模型，良好地解决尺度转换的技术问题，但由于缺乏大量相关和连续观测的数据，使得模型的参数量化存在不稳定因素，影响估测效果。随着地统计学、空间信息技术(遥感技术、地理信息系统技术)的发展，进一步推动了土壤有机碳尺度性估测与分布的研究，尤其是遥感以空间尺度性和时间动态性的时空技术优势，为资源环境的监测提供了重要技术方法，使得基于野外调查取样的地统计学估测、GIS空间插值估测、遥感尺度性反演等测定方法，成为当前土壤有机碳的常用的测定方法。

应用遥感估算土壤有机碳或其他土壤养分的长期目标是确定和推导植被冠层覆盖下的土壤养分特征，土壤反射光谱包含着丰富的土壤信息，从中可获取土壤有机质含量、氧化铁含量、土壤质地、主导粘土矿物类型等多种有用信息。据此，人们利用遥感技术测定土壤有机碳含量，并取得了一定的成果。但由于林地上的森林覆被物丰富，植被冠层覆盖下的土壤理化性质(例如土壤有机碳、有机质、养分含量和排水类型)不能直接使用卫星或航空影像数据来观测和计算。

到目前为此，利用遥感技术测定土壤有机碳含量或其他土壤养分特征，基本上是基于裸土地或农业用地上的研究和应用。选择合适的替代法是实现基于野外调查样本取样的林地土壤尺度性遥感估测的第一步。当前用于替代法估测土壤有机碳的生态过程模型较多，如何选择适合区域的，通过遥感可较易实现的，是研究实现土壤有机碳的关键技术。

发明内容

本发明的目的在于提供一种林地土壤有机碳遥感估测的方法，旨在借助于遥感与碳循环过程模型，通过替代法思维，实现林地土壤有机碳的的时空领域估测与反演，并解决传统土壤有机碳估测模型中，一些具有空间异质性的模型参数被简化为常数，而降低了估测精度的不足的问题。

本发明是这样实现的，一种林地土壤有机碳遥感估测的方法，包括以下步骤：

S1、在野外林地中取土壤并测定土壤有机碳含量；

S2、利用遥感数据和相关的气象数据，根据碳循环过程CASA模型，提取NPP的空间分布专题信息；

S3、依据生态系统碳循环处于平衡状态，植被净第一性生产力固定的有机碳，与通过土壤异养呼吸释放的有机碳相等的基础理论，即在测算出NPP的基础上，测算土壤基础呼吸，完成区域的土壤基础呼吸空间反演；

S4、利用反演出土壤基础呼吸的空间分布，结合野外所采集实验测定的有机碳含量值进行拟合，构建林地土壤有机碳遥感估测模型。

优选地，所述步骤S1具体为：在野外林地中选择具有代表性的地形和植被条件设采样点，根据分层采样法，在每个样点在剖面深度0～20cm、20～40cm、40～60cm取土壤均匀混合后用重络酸钾氧化-稀释热比色法测定土壤有机碳含量。

优选地，在步骤S2中，所述采用CASA模型用函数定义为：

式(1)中，NPP为植被性第一生产力；PAR光合有效辐射能，是从地表太阳辐射估算得出PAR＝0.5R，R为月平均太阳辐射观测数据；f(NDVI)为基于遥感影像标准化准差值植被(NDVI)的函数；ε^*光能利用率，是具体的常数；g(T)为温度胁迫因子；为水分胁迫因子，两者分别为湿度和温度的函数，表示温度和水分对土壤植被净初级生产力的限制影响。

优选地，所述光合有效辐射能PAR用函数定义为：

PAR＝SOL(x，t)×0.5     (2)

式(2)中，SOL(x，t)是t月份象元x位置的太阳总辐射量(MJ·m^-2)，常数0.5表示能够被植被所利用的光合有效辐射量约占太阳总辐射的一半比例，有效辐射波长大体位于0.4～0.7μm的区间段。

优选地，所述光能利用率ε^*用函数定义为：

ε(x，t)＝g(T)×h(w)×ε_max        (3)

式(3)中，ε(x，t)为光能利用率，g(T)为温度胁迫因子，h(w)为水分胁迫因子，ε_max为最大光能利用率。

优选地，所述最大光能利用率ε_max为BGC生态生理过程模型模拟的结果，将研究区森林类型图中的树种进行合并，赋予不同的植被类型以不同的最大光能利用率，得到研究区最大光能利用率专题信息。

优选地，所述温度胁迫因子g(T)的确定包括以下步骤：

g(T)是温度的函数，由构成。

Tε1(x)＝0.8+0.02Topt(x)-0.0005[Topt(x)]²      (4)

Tε2(x，t)＝1.1814/{1+e^{[0.2(Topt(x)-10-T(x，t))]}}/{1+e^{[0.3(-Topt(x)-10+T(x，t))]}}  (5)

式(4)以及(5)中，Topt(x)表示区域内月平均NDVI值达到最高时所对应的月平均气温，也被称为最适宜温度；T(x，t)表示某一个月的平均气温。

优选地，所述水分胁迫因子的确定包括以下步骤：

(1)建立局地潜在蒸散量Ep的多元回归模型，用函数定义为：

Ep＝2037.98-18.8308LAT-4.5801LONG-0.157861ALT   (6)

式(6)中，LAT为纬度，LONG为经度，ALT为海拔；

(2)通过陆地表面净辐射与潜在蒸散和降水量关系函数式估算地面净辐射量RN；

(3)测算区域实际蒸散量EET，用函数定义为：

$\mathrm{EET}=\frac{r\·\mathrm{Rn}(r^2+{\mathrm{Rn}}^2+r\·\mathrm{Rn})}{(r+\mathrm{Rn})\·(r^2+{\mathrm{Rn}}^2)}---\left(7\right)$

式(7)中，r为降水量；Rn为净辐射量；

(4)测算区域潜在蒸散量PET，用函数定义为：

Ep+EET＝2PET     (8)

(5)以水分胁迫因子变量作为环境水分因子的函数，由水分胁迫影响系数Wε表示，用函数定义为：

W_ε＝0.5+(EET/PET)    (9)。

优选地，在步骤S3中，所述区域的土壤基础呼吸的测定用函数定义为：

$A=\frac{\mathrm{NPP}}{e^{(b\×T)}\×y}---\left(10\right)$

式(10)中：A为土壤基础呼吸，NPP为植被性第一生产力，y为水分对土壤呼吸的限制因子，T为平均气温，b为温度敏感性系数。

优选地，在步骤S3中，所述土壤基础呼吸限制因子y用函数定义为：

$y=\frac{1}{1+30.0\×e^{-8.5x}}---\left(11\right)$

式(11)中，PPT为年降水量，PET为年潜在蒸散量；

所述温度敏感性系数b用函数定义为：

b＝lnQ10/10     (12)

式(12)中，Q10为土壤呼吸的温度敏感性，指温度升高10℃时，土壤呼吸强度增大的倍数。

本发明克服现有技术的不足，提供一种林地土壤有机碳遥感估测的方法，由于林地植被覆盖区遥感图像像元则表现为植被的光谱特征，遥感影像植被指数NDVI是目前应用最为广泛和适用的一个植被指数，广泛应用于植被性第一性生产力的估测研究，为此，在本发明中，通过CASA生态过程模型，在应用NDVI遥感数据(利用高分辨遥感影像反映区域差异性)，解决传统土壤有机碳估测模型中，一些具有空间异质性的模型参数被简化为常数，而降低了估测精度的不足，以及土壤有机碳储量有密切关系的土壤基础呼吸系数，建立土壤有机碳储量与土壤基础呼吸的关系模型，实现林地土壤有机碳的遥感估测研究。

附图说明

图1是本发明实施例中研究区地理位置示意图；

图2是本发明林地土壤有机碳遥感估测的方法的步骤流程图；

图3是本发明本发明实施例中野外样点分布图；

图4是本发明实施例中永安市1999年温度胁迫因子专题图；

图5是本发明实施例中永安市2008年温度胁迫因子专题图；

图6是本发明实施例中水分胁迫因子求取过程图；

图7是本发明实施例中土壤基础呼吸-有机碳拟合方程图；

图8是本发明实施例中1999年永安市林地有机碳分布专题图；

图9是本发明实施例中2008年永安市林地有机碳分布专题图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

在本发明实施例中，选择具体研究区进行说明。研究区概况：永安市地处福建省闽中偏西的三明市(图1)，东经116°56′～117°47′，北纬25°33′～26°12′，东西宽约82km，南北长约71km，土地总面积294110hm²。永安素有“九山半水半分田”之称，境内群山叠拥，山地面积占全市面积的59％，达到174900hm2，丘陵面积约111000hm²，占全市土地面积38％，拥有着较早上市的永安林业(集团)股份有限公司，山地资源开发利用程度总体较高，作为福建省重要经济发展区域，是我国南方集体林重要改革试验区、用材林重要产区。气候为亚热带季风性湿润气候，降水较为丰沛，年均降雨量可达1490～2050mm，光照充足，年均气温可达14.3～19.2℃。

一种林地土壤有机碳遥感估测的方法，如图2所示，包括以下步骤：

S1、在野外林地中取土壤并测定土壤有机碳含量；

在步骤S1中，更具体的，野外林地实测土壤有机碳含量，采样时间段为2010年8月15日至9月10日。为了减少降雨等因素的影响，保证研究区野外采样前三日无降水。依据研究区二类调查固定样地分布点，选择具有代表性的地形和植被条件的样点116个(图3)，剔除记录错误的数据，剩余104个数据用于建模与验证。根据分层采样法，每个样点在剖面深度0～20cm、20～40cm、40～60cm共取1kg土壤均匀混合。实验室内利用重络酸钾氧化-稀释热比色法测定土壤有机碳(王冬梅主编.土壤理化分析实验指导[M].北京：中国林业出版社，2002.)。

S2、利用遥感数据和相关的气象数据，根据碳循环过程CASA模型，提取NPP的空间分布专题信息

在步骤S2中，遥感数据研究采用与调查时间同步的ALOS影像。利用ERDAS9.2对遥感图像进行几何精度校正，在此基础上利用ENVI面向对象法结合小班图进行监督分类，得到混交林、常绿阔叶林、常绿针叶林的分布信息。气象数据主要包括月平均气温，年均降水量和月平均太阳辐射数据。利用GIS克里格空间插值法将其数据化，获得与遥感图像投影相同、输出象元大小一致的气象要素栅格图。

更具体的技术思路如下所示：

遥感与生态系统碳平衡方程相结合的方法能够有效地反演土壤异养呼吸的空间分布，进而提高土壤有机碳模型的预测精度。当大气、植被、土壤处于平衡状态下时，通过植被净第一性生产力固定的有机碳与土壤异养呼吸释放的有机碳(土壤基础呼吸)相等，作为生态系统碳循环处于平衡状态下的理论值。此时可以根据公式(周涛，史培军，王绍强.气候变化及人类活动对中国土壤有机碳储量的影响[J].地理学报，2003，58(5)：727-734)来计算土壤基础呼吸。

采用CASA模型作为估算植被净第一性生产力的方法。模型的主要参数如下表1所示：

表1各参数的确定方法

CASA模型是从植被的生理过程出发而建立的植被净第一性生产力机理模型，该模型考虑了植被净第一性生产力计算的两个主要驱动变量：植被所吸收的有效光合辐射与光能利用率。而这两个变量又可以通过NDVI、土壤水分、降水量、平均温度等指来体现。大部分过程模型比较复杂，参数比较多，由于一些变量的定量化较困难，通常须人为指定某些参数或对模型进行简化。相比之下，CASA模型所需要的参数相对于其他模型来说较少，这样就避免了由于参数缺乏而人为简化模型或者随研究地区差异而产生的误差。

(1)利用遥感数据(NDVI)和相关的气象数据(温度、降水、太阳辐射)，根据碳循环过程CASA模型，实现NPP的空间分布专题信息提取：

式中：NPP为净初级生产力；PAR光合有效辐射能，是从地表太阳辐射(R)估算得出(PAR＝0.5R)，R为月平均太阳辐射观测数据；f(NDVI)为基于遥感影像标准化准差值植被(NDVI)的函数；ε^*光能利用率，是具体的常数；g(T)为温度胁迫因子；为水分胁迫因子，两者分别为湿度和温度的函数，表示温度和水分对土壤植被净初级生产力的限制影响。

(2)当生态系统碳循环处于平衡状态时，植被净第一性生产力固定的有机碳，与通过土壤异养呼吸释放的有机碳相等，则可在测算出NPP的基础上，测算土壤基础呼吸，实现区域的土壤基础呼吸空间反演，土壤基础呼吸计算用函数定义为：

$A=\frac{\mathrm{NPP}}{e^{(b\×T)}\×y}$

式中：A为土壤基础呼吸，NPP为植被性第一生产力，y为水分对土壤呼吸的限制因子，T为平均气温。

(3)林地NPP的估测

A、太阳辐射参数(PAR)的确定

太阳辐射能是森林植被进行光合作用能量来源，植物利用其中部分波长区段的辐射能合成自身有机物，这部分被植物利用的辐射能称为光合有效辐射(Photosynthetically Active Radiation)，简称PAR，利用以下公式求得：

PAR＝SOL(x，t)×0.5

式中，SOL(x，t)是t月份象元x位置的太阳总辐射量(MJ·m^-2)，常数0.5表示能够被植被所利用的光合有效辐射量约占太阳总辐射的一半比例(MeCree，1981)，有效辐射波长大体位于0.4～0.7μm的区间段。

目前，国内外由于对PAR的测量和研究较为薄弱，其气候学计算方法的理论基础还不够扎实，根据经验公式间接推算，常用的方法便是采用总辐射乘以光合有效辐射系数的经验公式。目前国内学者利用CASA模型进行NPP估算时，均采用0.5作为光合有效辐射系数的取值，为此研究确定以0.5为光合有效辐射系数的取值。

B、植被对光合有效辐射吸收比例的确定

绿色植物对光能的吸收利用必须依靠叶片的光合作用来完成。f(PAR)表示植被对入射光合有效辐射(PAR)的吸收比例，是由植被指数构成的函数。植被指数能够简单、有效度量地表的植被状态，成为估测植被净第一性生产力研究的重要基础。常用的植被指数如下表2所示：

表2常见的植被指数及其计算公式

其中：NIR为近红外波段，R为红外波段，G为绿光波段。分别对应ALOS遥感影像的第4、3、2波段。

目前国内外的研究基本采用归一化植被指数NDVI(Normalized DifferenceVegetation Index)与SR指数(基于归一化植被指数NDVI指数的比值植被指数)，研究前后综合尝试了以上9种植被指数，除了GEMUI(全球环境监测植被指数)的参数可能不适合我国南方植被研究以外，与其他6种植被指数相比而言，NDVI指数与SR指数确实更为适合于项目研究。同时考虑影像质量的影响，研究采用95％的置信率，对提取的植被指数图进行处理，以增强数据质量的稳定性和可靠性。为此，先根据f(NDVI)函数式求出由NDVI直接估算的f(NDVI)值和f(SR)函数式求出由SR直接估算的f(SR)值，进而得到参数f(PAR)：

$f\left(\mathrm{NDVI}\right)=\frac{[\mathrm{NDVI}(x,t)-\mathrm{NDVI}(i,\min )\×(\mathrm{FPAR}\max -\mathrm{FPAR}\min )]}{[\mathrm{NDVI}(i,\max )-\mathrm{NDVI}(i,\min )]}+\mathrm{FPAR}\min$

式中：i表示植被类型，x表示像元，t表示月份；FPARmin和FPARmax的取值与植被类型无关，分别为0.001和0.95。

再由求出SR(比值植被指数)，进而求出由NDVI间接估算的f(NDVI)值。

$\mathrm{FPAR}(x,t)=\frac{[\mathrm{SR}(x,t)-\mathrm{SR}(i,\min )]\×[\mathrm{FPAR}\max -\mathrm{FPAR}\min ]}{\mathrm{SR}(i,\max )-\mathrm{SR}(i,\min )}+\mathrm{FPAR}\min$

研究表明，由NDVI直接估算的f(PAR)值比实测值高，而SR估算的f(PAR)值低于实测值。考虑到这种情况，参考前人研究将这两种方法结合起来：

$f\left(\mathrm{PAR}\right)=\∂\·f\left(\mathrm{NDVI}\right)+(1-\∂)\·f\left(\mathrm{SR}\right)$

式中，值是待定值，目前国内较多采用取均值的方法，若取0.5，表示二者的均值，研究对永安地区试验，发现永安取0.5结果较为理想。

C、光能利用率的确定

光能利用率是指植被将叶片吸收的有效光合辐射转化为有机碳的效率。在理想条件下植被应当具有最大光能利用率ε^*，而在实际条件下，受环境因素的影响，植被往往达不到最大光能利用率。现实条件下的植被实际光能利用率主要由温度和水分胁迫因子来调整最大光能利用率而确定。因此，在计算植被实际光能利用率时，需要考虑温度和水分胁迫因子对光能利用率的影响。

ε(x，t)＝g(T)×h(w)×ε_max

式中，ε(x，t)为光能利用率，g(T)为温度胁迫因子，h(w)为水分胁迫因子，ε_max为最大光能利用率。

(最大光合利用率)是指植被在没有任何限制的理想条件下对光合有效辐射的利用率，它是植被本身的一种生理属性，由其内在的生物学机制形成，其取值因不同的植被类型而有所不同，针对某一类型的植被而言是常数。通过试验比较分析，研究中最大光能利用率的取值采用Running等学者(参考文献：Running S W，Hunt E R.Generalization of a forest ecosystem process model forother biomes，Biome-BGC，and an application for global-scale models[M]||Ehleringer J R，Field C B.Scaling Physiological Processes：Leaf to Globe.SanDiego：Academic Press，1993.)根据BGC生态生理过程模型模拟的结果，将永安森林类型图中的树种进行合并，赋予不同的植被类型以不同的最大光能利用率，得到研究区最大光能利用率专题信息，具体取值如下表3所示：

表3不同植被的最大光能利用率

D、温度胁迫因子g(T)

温度胁迫因子g(T)，表示温度对最大光合利用率的影响，是对最大光合利用率的修正乘积指数。g(T)是温度的函数，由构成。

Tε1(x)＝0.8+0.02Topt(x)-0.0005[Topt(x)]²

Tε2(x，t)＝1.1814/{1+e^{[0.2(Topt(x)-10-T(x，t))]}}/{1+e^{[0.3(-Topt(x)-10+T(x，t))]}}

Topt(x)表示区域内月平均NDVI值达到最高时所对应的月平均气温，也被称为最适宜温度。T(x，t)表示某一个月的平均气温。在使用公式时，若某一月平均温度T(x，t)比最适宜温度Topt(x)高10℃或低13℃时，该月的值应计算为月平均温度T(x，t)为最适宜温度Topt(x)时值的一半。

利用GIS空间插值得到1月和9月的月均温栅格图(1月为图像遥感拍摄时刻；9月为最适宜温度产生的月份)。由于最适温度是NDVI出现最大值时所对应的月平均温度，因而确定最适温度的时间为NDVI最大值出现的时间。NDVI最大值多出现在6～9月份，8月最多，占34.6％；余下依次为7月(18.7％)、9月(13.8％)和6月(6.7％)。从植被分布的类型来看，耕地的NDVI最大值多出现在7、8月份；草原的NDVI最大值多出现在8月份；全国大部分地区的常绿针叶林和常绿阔叶林的NDVI最大值一般出现在9月份，但西南的部分地区多出现在11月份；落叶针叶林和落叶阔叶林的NDVI最大值多出现在6月份。福建属于中国东南地区，其气候类型决定了植被类型以常绿阔叶林和常绿针叶林为主，为此研究确定福建省最适宜的植被生长温度为9月份。

利用相关研究基础资料，根据上述式建模计算，永安市的取值范围在0.987～1，具体如图4～5所示。

E、水分胁迫因子

h(w)变量作为环境水分因子的函数，由水分胁迫影响系数Wε表示。

W_ε＝0.5+(EET/PET)

式中，EET(x，t)为区域实际蒸散量(mm)，根据下式测算：

$\mathrm{EET}=\frac{r\·\mathrm{Rn}(r^2+{\mathrm{Rn}}^2+r\·\mathrm{Rn})}{(r+\mathrm{Rn})\·(r^2+{\mathrm{Rn}}^2)}$

式中：r为降水量(mm)；Rn为净辐射量(mm)。当区域潜在蒸散量小于月平均降水量时，区域实际蒸散量EET(x，t)＝区域潜在蒸散量PET(x，t)。

PET(x，t)为区域潜在蒸散量(mm)，可根据下式的互补关系求得：

Ep+EET＝2PET

式中，EET为下垫面的区域实际蒸散量，Ep为局地潜在蒸散量(即小块充分湿润地面的蒸散量)。区域潜在蒸散量PET为大面积充分湿润陆面的蒸散量。当区域趋于愈来愈湿润的状态时，EET和Ep都以PET为限。

应用并求解互补关系时，可以根据由小蒸发面的资料归纳出来的Penman公式计算Ep，也可以借鉴植被-气候关系模型来计算Ep，植被-气候关系模型认为Ep作为热量指标，与我国的纬度(LAT)、经度(LONG)及海拔(ALT)存在高度相关，其多元回归模型如下(R＝91.7％，N＝671)：

Ep＝2037.98-18.8308LAT-4.5801LONG-0.157861ALT

式中，LAT为研究地的纬度、LONG为研究地的经度，ALT为研究地的海拔，依据上式，在我国Ep随经度每向西一度则约增加4.6mm，随纬度每向北一度则减少18.8mm，且随海拔升高而降低。Ep向北与向高海拔的趋低特征与热量递减是相关的。

Rn为地面净辐射量，由于大多数气象观测站均不进行地表净辐射观测工作，计算地表净辐射需要的气象要素也不易获取。因此，本研究利用式14中的陆地表面净辐射与潜在蒸散和降水量关系函数式，估算地面净辐射量：

Rn＝(Ep·r)^0.5·[0.369+0.598·(Ep/r)^0.5]

综上，水分胁迫因子的具体求取流程，如图6所示。

由于研究地永安市区域范围较小，经纬度按栅格取值没有意义(其极端误差不到1％)。因此经纬度取研究区的平均经纬度值作为Ep公式中的经纬度，其中永安为北纬25.88°、东经117.35°.而ALT由DEM图层直接提取。

在此基础上，根据CASA模型(即)，通过建模，获取研究区与遥感影像相对应的植被净第一性生产力：永安市1月的植被净第一性生产力。

S3、依据生态系统碳循环处于平衡状态，植被净第一性生产力固定的有机碳，与通过土壤异养呼吸释放的有机碳相等的基础理论，即在测算出NPP的基础上，测算土壤基础呼吸，完成区域的土壤基础呼吸空间反演

在步骤S3中，土壤呼吸的测定用函数定义为：

$A=\frac{\mathrm{NPP}}{e^{(b\×T)}\×y}$

式中，A为土壤基础呼吸(gm^-2)。y为水分对土壤呼吸的限制因子，介于(0.03，1)。当降水量大于或等于潜在蒸散量时，y≈1。T为气温温度，b为温度敏感性系数。温度敏感性系数与温度的乘积作为温度敏感度。y根据以下式子求得：

$y=\frac{1}{1+30.0\×e^{-8.5x}}$

式中，PPT为年降水量(r)，PET为年潜在蒸散量(即Ep0)。三明市各县年降水量平均在1400～2000mm之间，均属福建省丰水区。研究区永安2008年降水量充沛，y基本在0.99～1之间。

b根据以下式子求得：

b＝lnQ10/10

式中：Q10表示土壤呼吸的温度敏感性。它是指温度升高10℃时，土壤呼吸强度增大的倍数，是生物地球化学模型描述碳循环过程的一个重要参数。土壤呼吸的温度敏感性是模拟全球变暖与生态系统碳释放之间反馈强度的重要参数。虽然研究表明Q10值具有明显的空间异质性，但由于其空间分布格局的定量数据的缺乏，目前绝大多数生物地球化学模型将其简化成一个常数。综合国内外许多研究成果，从植被类型来看，森林植被土壤呼吸的温度敏感性最大，标准差相对较小(2.11±0.43)，草地土壤呼吸的温度敏感性较低，具有相对较高的标准差(2.01±0.61)；农田土壤的Q10值介于森林土壤与草地土壤之间(2.08±0.45)。对整个中国而言，平均Q10值为1.80。从不同土壤类型来看，Q10均值波动范围在1.09～2.38之间。目前中国范围内对不同土壤类型的温度敏感性Q10值的主要值见表4。

表4不同土壤类型的Q10取值

全国红壤的Q10平均值为1.87，符合研究区永安市所有林地土壤类型的全国土壤Q10平均值的变化范围在1.87～2.1之间。但以红壤为主福建省地处中国东南部，是典型的亚热带气候，植被覆盖率较高，其土壤温度敏感性应大于全国红壤平均水平1.87，两个研究区的气候均属于典型的亚热带气候，植被覆盖率较高，其土壤温度敏感性应大于全国平均水平。结合福建临近经纬度的Q10观测值加以综合考虑，取森林土壤Q10均值2.11，作为永安地区平均Q10的取值。该取值也在福建临近经纬度的Q10观测值(1.88～2.24)的范围内(表5)，符合众多学者在福建临近经纬度范围内的Q10观测结果。

表5福建临近经纬度的Q10观测值

综合上述获取各专题信息基础上，实现永安市研究区土壤基础呼吸的测定与遥感反演。

S4、利用反演出土壤基础呼吸的空间分布，结合野外所采集实验测定的有机碳含量值进行拟合，构建林地土壤有机碳遥感估测模型

在步骤S4中，利用上述研究反演出土壤基础呼吸的空间分布，结合野外所采集实验测定的有机碳值进行拟合，构建林地土壤有机碳遥感估测模型。考虑阴影区和非阴影区光谱值差异较大，直接构建模型误差大。为此，分非阴影区和阴影区分别建立林地土壤有机碳估测模型，在实现非阴影区和阴影区林地土壤有机碳估测基础上，叠加得到研究区的林地土壤有机碳专题图。

永安市非阴影区林地土壤有机碳拟合方程为：

Y＝0.7533X+2.4854(R＝0.76)

式中，Y为估测有机碳，X为土壤基础呼吸。

利用该方程，对剩余15个点进行验证，平均精度为83％，如图7所示。

阴影区中NDVI值并不与植被状况成正比或反比，同时受阴影深浅的影响，为此通过试验再分析，在土壤基础呼吸与实测有机碳拟合时加入坡向地形因子(坡向分为8个等级)，通过数量化理论拟合得到阴影区的林地土壤有机碳估测方程：

y＝18.14791855+1.4085772191X1-0.9017582558X2(R＝0.75)

式中，X1表示土壤基础呼吸值，X2表示坡向。

利用上述非阴影区林地拟合的林地土壤有机碳估测方程反演得到永安市的林地土壤有机碳专题图，计算平均精度81.2％，如下表6所示：

表6阴影区模型精度验证表

相比与现有技术的缺点和不足，本发明具有以下有益效果：

(1)本发明摆脱树种覆盖的影响，利用遥感与碳循环过程模型结合，采用替代法思维实现了林下土壤有机碳的估测，为遥感估测其他林地土壤属性特征提供了技术借鉴和理论参考。

(2)本发明利用遥感的时空技术优势，可实现林地土壤有机碳的历史反演，为森林经营过程对土壤的影响与测定、分析提供了重要的基础数据。

(3)本发明利用遥感的空间尺度良好地反映了林地土壤有机碳的空间差异性，解决传统土壤有机碳估测模型中，一些具有空间异质性的模型参数被简化为常数，而降低了估测精度的不足。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种林地土壤有机碳遥感估测的方法，其特征在于包括以下步骤：

S1、在野外林地中取土壤并测定土壤有机碳含量；

S2、利用遥感数据和相关的气象数据，根据碳循环过程CASA模型，提取NPP的空间分布专题信息；

S3、依据生态系统碳循环处于平衡状态，植被净第一性生产力固定的有机碳，与通过土壤异养呼吸释放的有机碳相等的基础理论，即在测算出NPP的基础上，测算土壤基础呼吸，完成区域的土壤基础呼吸空间反演；

S4、利用反演出土壤基础呼吸的空间分布，结合野外所采集实验测定的有机碳含量值进行拟合，构建林地土壤有机碳遥感估测模型。

2.如权利要求1所述的林地土壤有机碳遥感估测的方法，其特征在于，所述步骤S1具体为：在野外林地中选择具有代表性的地形和植被条件设采样点，根据分层采样法，在每个样点在剖面深度0～20cm、20～40cm、40～60cm取土壤均匀混合后用重络酸钾氧化-稀释热比色法测定土壤有机碳含量。

3.如权利要求1所述的林地土壤有机碳遥感估测的方法，其特征在于，在步骤S2中，所述采用CASA模型用函数定义为：

式(1)中，NPP为植被性第一生产力；PAR光合有效辐射能，是从地表太阳辐射估算得出PAR＝0.5R，R为月平均太阳辐射观测数据；f(NDVI)为基于遥感影像标准化准差值植被(NDVI)的函数；ε^*光能利用率，是具体的常数；g(T)为温度胁迫因子；为水分胁迫因子，两者分别为湿度和温度的函数，表示温度和水分对土壤植被净初级生产力的限制影响。

4.如权利要求3所述的林地土壤有机碳遥感估测的方法，其特征在于，所述光合有效辐射能PAR用函数定义为：

PAR＝SOL(x，t)×0.5       (2)

式(2)中，SOL(x，t)是t月份象元x位置的太阳总辐射量(MJ·m^-2)，常数0.5表示能够被植被所利用的光合有效辐射量约占太阳总辐射的一半比例，有效辐射波长大体位于0.4～0.7μm的区间段。

5.如权利要求3所述的林地土壤有机碳遥感估测的方法，其特征在于，所述光能利用率ε^*用函数定义为：

ε(x，t)＝g(T)×h(w)×ε_max      (3)

式(3)中，ε(x，t)为光能利用率，g(T)为温度胁迫因子，h(w)为水分胁迫因子，ε_max为最大光能利用率。

6.如权利要求5所述的林地土壤有机碳遥感估测的方法，其特征在于，所述最大光能利用率ε_max为BGC生态生理过程模型模拟的结果，将研究区森林类型图中的树种进行合并，赋予不同的植被类型以不同的最大光能利用率，得到研究区最大光能利用率专题信息。

7.如权利要求3或5所述的林地土壤有机碳遥感估测的方法，其特征在于，所述温度胁迫因子g(T)的确定包括以下步骤：

g(T)是温度的函数，由构成；

Tε1(x)＝0.8+0.02Topt(x)-0.0005[Topt(x)]²      (4)

Tε2(x，t)＝1.1814/{1+e^{[0.2(Topt(x)-10-T(x，t)]}}/{1+e^{[0.3(-Topt(x)-10+T(x，t))]}}  (5)

式(4)以及(5)中，Topt(x)表示区域内月平均NDVI值达到最高时所对应的月平均气温，也被称为最适宜温度；T(x，t)表示某一个月的平均气温。

8.如权利要求3或5所述的林地土壤有机碳遥感估测的方法，其特征在于，所述水分胁迫因子的确定包括以下步骤：

(1)建立局地潜在蒸散量Ep的多元回归模型，用函数定义为：

Ep＝2037.98-18.8308LAT-4.5801LONG-0.157861ALT      (6)

式(6)中，LAT为纬度，LONG为经度，ALT为海拔；

(2)通过陆地表面净辐射与潜在蒸散和降水量关系函数式估算地面净辐射量RN；

(3)测算区域实际蒸散量EET，用函数定义为：

$\mathrm{EET}=\frac{r\·\mathrm{Rn}(r^2+{\mathrm{Rn}}^2+r\·\mathrm{Rn})}{(r+\mathrm{Rn})\·(r^2+{\mathrm{Rn}}^2)}---\left(7\right)$

式(7)中，r为降水量；Rn为净辐射量；

(4)测算区域潜在蒸散量PET，用函数定义为：

Ep+EET＝2PET      (8)

(5)以水分胁迫因子变量作为环境水分因子的函数，由水分胁迫影响系数Wε表示，用函数定义为：

W_ε＝0.5+(EET/PET)   (9)。

9.如权利要求1所述的林地土壤有机碳遥感估测的方法，其特征在于，在步骤S3中，所述区域的土壤基础呼吸的测定用函数定义为：

$A=\frac{\mathrm{NPP}}{e^{(b\×T)}\×y}---\left(10\right)$

式(10)中：A为土壤基础呼吸，NPP为植被性第一生产力，y为水分对土壤呼吸的限制因子，T为平均气温，b为温度敏感性系数。

10.如权利要求9所述的林地土壤有机碳遥感估测的方法，其特征在于，在步骤S3中，所述土壤基础呼吸限制因子y用函数定义为：

$y=\frac{1}{1+30.0\×e^{-8.5x}}---\left(11\right)$

式(11)中，PPT为年降水量，PET为年潜在蒸散量；

所述温度敏感性系数b用函数定义为：

b＝lnQ10/10     (12)

式(12)中，Q10为土壤呼吸的温度敏感性，指温度升高10℃时，土壤呼吸强度增大的倍数。