CN101908196B - 基于植被-土壤水分响应关系的植被盖度估算方法 - Google Patents

Abstract

本发明公开了基于NDVI和NDVI对TVDI响应敏感性的植被-土壤水分响应盖度估算方法。植被盖度包括裸地盖度fc_b、草地盖度fc_g、林地和灌丛混合类型盖度fc_f/s。fc_f/s包括林地盖度fc_forest和灌丛盖度fc_shrub。该方法包括步骤：确定研究区的每一像元的归一化差值植被指数NDVI和TVDI值；确定研究区的每个像元的NDVI值对TVDI敏感性α值；基于研究区中的每一个像元的α值和NDVI值构建所述研究区域的三角形散点图；基于三角形散点图确定三角形的每一个端点的端元α值和NDVI值；利用最小二乘法，基于所确定的每一个三角形端点的端元的α值和NDVI值并且在约束条件0≤fc_b≤1，0≤fc_g≤1，0≤fc_f/s≤1下，对研究区内每一像元进行线性分解以获得研究区内每一像元的fc_b、fc_g、fc_f/s。

Description

基于植被-土壤水分响应关系的植被盖度估算方法

技术领域

本发明总的来说涉及植被盖度，具体地说涉及基于植被-土壤水分响应关系估算植被盖度的方法。

背景技术

较早利用遥感数据估算植被盖度的模型是经验模型。一种经验模型是利用野外实测植被盖度与遥感信息建立线性或非线性的经验关系模型，然后将该模型应用到整个研究区域。此类经验模型都是简单地在地面实测数据与遥感数据之间建立一定的关系，将地面的实测数据简单地扩展到遥感信息的二维空间中。由于反射光谱的线性和非线性组合而形成的能反映绿色植被生长状况和分布的特征指数(植被指数)可以用来监测植物的各种生长参数，并且植被盖度与植被指数表现出显著的相关性。因此，近年大多数研究根据地表实测植被盖度与基于遥感影像不同波段构建的植被指数之间的显著相关性，建立另一种植被盖度与植被指数经验关系模型。样方尺度、关系模型的选择(线性模型与多项式模型)对于植被盖度经验模型精度产生影响较大。经验模型对特定区域的地表实测数据具有依赖性，当研究区尺度较小时，测量结果具有一定的精度；当研究区尺度较大时，精度就会大大降低。

可以用来动态监测植被生长状况的多时相遥感影像空间分辨率较低，图像中的一个像元可由多个组分构成，即存在混合像元，每个组分对传感器所获取的信息都有贡献，每个像元的遥感信息可以分解为多个组分，一个像元中所有组分通过线性和非线性组合，构成像元的遥感信息。因此，可以把每一像元信息进行线性或非线性分解，建立像元分解模型，并以此模型估算植被覆盖度。针对中低分辨率遥感数据发展了两类像元分解模型：对遥感影像的不同波段数据和基于不同波段构成的植被指数进行线性和非线性的混合分解。基于遥感数据不同波段的像元分解方法主要是光谱分解模型。线性和非线性光谱混合模型是利用一个线性或非线性关系表达遥感系统中一个像元内各地物的类型、比例与地物的光谱响应。像元在某一光谱波段的反射率是由构成像元的基本组分的反射率及其所占像元面积比例为权重系数的线性或非线性组合。与传统的分类方法相比，混合像元分解后，由不同地物类型组成的一个混合像元被分解成像元内各个组分，各个植被组分占像元的百分比即为每种植被类型盖度。可以看出，像元分解方法提高了用一般植被盖度估算方法估算植被盖度的精度。近年来在光谱分解模型的基础上，发展了利用植被指数的分解模型，应用较早的是亚像元分解模型，根据遥感影像像元的特点，分析亚像元结构的分布特征，针对不同的亚像元结构，建立不同的植被盖度模型。亚像元分解模型为大面积植被盖度的估算提供了一种有效的途径，基本能够满足大尺度生态及气候模型研究的要求，由于亚像元分解模型是在若干假设和近似下建立起来的模型，而植被指数与植被盖度、叶面指数之间又存在着较为复杂的关系，因此难以利用亚像元分解模型获得高精度的植被覆盖率

已有的植被盖度估测模型考虑了遥感影像对植被光谱特征的综合反映，通过对各像元植被类型及分布特征分析，建立不同的植被盖度模型，大部分模型只是建立光谱特性对植物响应的基础上，没有考虑气候因素或其他因素对植被盖度的影响，但由于受地表蒸发、植物蒸腾等因素的影响，土壤水分空间变异较大，将降水量作为分解模型的状态变量，降低了模型估算精度。土壤水分是地表土壤水分影响到植物和作物的生长，是监测土地退化的一个重要指标，是气候、水文、生态，农业等领域的主要参数，在地表与大气界面的水分和能量交换中起重要作用。土壤水分降低导致植物生长量明显减退。因此，本方法将土壤水分作为影响植被盖度因子之一，引入植被盖度的估算模型，提出一种新的线性分解模型(植被-土壤水分响应模型)，将土壤水分作为模型的状态变量，建立新的植被盖度分解模型，估测研究区不同植被类型的盖度。此模型考虑了多年植被变化对土壤水分的敏感性，克服了以前的光谱混合模型依靠经验获取训练数据集的缺点，并在一定程度上消除了大气、传感器等引起的噪声点影响。此分解模型不仅可以用来进行当前的植被盖度估算，而且可以用来模拟未来特定条件下的植被盖度。

发明内容

根据本发明，在Ts/NDVI特征空间的概念分析基础上，计算地表干燥度指数(TVDI)，结合NDVI和NDVI对TVDI敏感性，提出基于NDVI和NDVI对TVDI响应敏感性的植被-土壤水分响应盖度估测模型(Vegetation-soil MoistureResponse Model，VSMRM)，从而估算像元内各植被类型组分的盖度。

根据本发明，提供了一种基于植被-土壤水分响应关系来估算一区域内各个像元中的植被盖度的方法，所述植被盖度包括裸地盖度fc_b，草地盖度fc_g、林地和灌丛混合类型盖度fc_f/s以及林地盖度fc_forest和灌丛盖度fc_shrub，该方法包括步骤：

(1)确定所述研究区的每一像元的归一化差值植被指数NDVI和TVDI值；

(2)确定所述研究区的每个像元的NDVI值对TVDI敏感性α值；

(3)基于所述研究区中的每一个像元的α值和NDVI值构建所述研究区域的三角形散点图，其中所述α值为自变量，NDVI值为因变量；

(4)基于三角形散点图确定三角形的每一个端点的端元的α值(α_f/s，α_g，α_b)和NDVI值(NDVI_f/s，NDVI_g，NDVI_b)；

(5)利用最小二乘法，基于所确定的每一个三角形端点的端元的α值和NDVI值并且在约束条件0≤fc_b≤1，0≤fc_g≤1，0≤fc_f/s≤1下，对研究区内每一像元进行线性分解以获得研究区内每一像元的fc_b、fc_g、fc_f/s。

在步骤(2)中，根据公式(11)-(13)确定每个像元的NDVI值对TVDI敏感性α值，

$\mathrm{\α}=\mathrm{dNDVI}/d\hat{T}---\left(11\right)$

$\hat{T}=(\mathrm{TVDI}-\stackrel{\‾}{T})/\mathrm{\σ}---\left(12\right)$

$\mathrm{\σ}=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{TVDI}}_i-\stackrel{\‾}{T})}^2/(n-1)}---\left(13\right)$

式中，T为所述研究区的每一个像元的多年TVDI时相平均值，σ为TVDI的标准差。

其中，每一个像元的植被盖度满足公式(8)-(10)，

fc_b+fc_g+fc_f/s＝1                       (8)

NDVI_bfc_b+NDVI_gfc_g+NDVI_f/sfc_f/s＝NDVI   (9)

α_bfc_b+α_gfc_g+α_f/sfc_f/s＝α    (10)

式中，NDVI_b，NDVI_g，NDVI_f/s分别表示每像元内裸地、草地和森林与灌丛混合类型对应的NDVI值，α_b，α_g，α_f/s分别表示每像元内裸地、草地、森林和灌丛混合类型对应的NDVI值对TVDI敏感性。

其中，在步骤(4)中，植被对干燥度指数敏感性较小、NDVI值最大的三角形端点为具有α_f/s，NDVI_f/s的林地和灌丛混合类型的端元，植被对气候干燥度指数敏感性最大、NDVI值较大的三角形端点为具有α_g，NDVI_g的草地的端元，植被对干燥度指数敏感性较小，NDVI值最小的三角形端点为具有α_b，NDVI_b的裸地的端元，

所述步骤(5)包括：

将所确定的每一像元的α值，NDVI值，以及所确定的α_f/s，NDVI_f/s，α_g，NDVI_g，α_b，NDVI_b的裸地的端元代入公式(8)-(10)中，

调用MatLab软件中的最小二乘函数对公式(8)-(10)进行线性分解以获得每一像元的裸地盖度fc_b，草地盖度fc_g以及林地和灌丛盖度fc_f/s。

所述基于植被-土壤水分响应关系来估算一区域内各个像元中的植被盖度的方法进一步包括对林地和灌丛混合类型盖度fc_f/s进一步分解，得到每一像元内林地盖度fc_forest和林地盖度fc_shrub步骤，该步骤包括：

(6)根据所述研究区的植被类型图，从该植被类型图的分类图像中选取一定数量的全为林地覆盖的纯像元，并且确定与这些全为林地覆盖的纯像元对应的NDVI值并且求取这些NDVI值的平均值，作为全为林地覆盖的纯像元的NDVI_f ^*值；

(7)根据所述研究区的植被类型图，从该植被类型图的分类图像中选取一定数量的全为灌丛覆盖的纯像元，确定与这些全为灌丛覆盖的纯像元对应的NDVI值并且求取这些NDVI值的平均值，作为全为灌丛覆盖的纯像元的NDVI_s ^*值；

(8)根据公式(15)-(17)确定每一像元的林地盖度fc_forest和灌丛盖度fc_shrub，

${\mathrm{fc}}_{\mathrm{forest}}=\frac{\stackrel{\‾}{{\mathrm{NDVI}}_{f/s}}{\mathrm{fc}}_{f/s}-{\mathrm{NDVI}}_s^{*}{\mathrm{fc}}_{f/s}}{{\mathrm{NDVI}}_f^{*}-{\mathrm{NDVI}}_s^{*}}---\left(15\right)$

$\stackrel{\‾}{{\mathrm{NDVI}}_{f/s}}={\mathrm{NDVI}}_{f/s}+\hat{T}*\mathrm{\α}---\left(16\right)$

fc_shrub＝fc_f/s-fc_forest    (17)

所述步骤(1)中，

确定所述研究区的每一像元的归一化差值植被指数NDVI值包括：

确定辐射校正参数和大气校正参数；

利用确定的确定辐射校正参数和大气校正参数以及大气校正软件6S对表观反射率进行校正，得到TM第3、第4波段的地表反射率；

根据公式(3)计算每一像元的NDVI值，

$\mathrm{NDVI}=\frac{R_{\mathrm{NIR}}-R_{\mathrm{RED}}}{R_{\mathrm{NIR}}+R_{\mathrm{RED}}}---\left(3\right)$

式中，R_NIR为TM第4通道近红外波段的光谱反射率，R_RED为TM第3通道红光波段的光谱反射率；

确定所述研究区的每一像元的TVDI值包括：

确定每一像元的地表温度T_s，

根据公式(7)计算每一像元的TVDI值，

$\mathrm{TVDI}=\frac{T_s-T_{S_{\min }}}{T_{S_{\max }}-T_{S_{\min }}}---\left(7\right)$

其中，Ts_min表示最小地表温度，对应的是湿边，Ts_min＝a₁+b₁NDVI，a₁、b₁是湿边拟合方程的系数；Ts是任意像元的地表温度；Ts_max为某一NDVI对应的最高温度，即干边，Ts_max＝a₂+b₂NDVI，a₂、b₂是干边拟合方程的系数。

其中所述确定辐射校正参数和大气校正参数包括：

利用公式(1)计算辐射校正参数L_λ，

$L_{\mathrm{\λ}}=\frac{L_{\max \mathrm{\λ}}-L_{\min \mathrm{\λ}}}{\mathrm{QCAL}-{\mathrm{QCAL}}_{\min }}\×(\mathrm{QCAL}-{\mathrm{QCAL}}_{\min })+L_{\min \mathrm{\λ}}---\left(1\right)$

式中，L_λ表示传感器所接收到的辐射强度，QCAL表示TM数据的像元灰度值(DN值)，QCAL_max为传感器所接收到的最大DN值，QCAL_min为传感器所接收到的最小DN值，Lmax_λ为传感器所接收到的的最大光谱辐射值，Lmin_λ为传感器所接收到的最小光谱辐射值；

利用公式(2)计算大气校正参数，

${\mathrm{\ρ}}_{\mathrm{App}}=\frac{{\mathrm{\πL}}_{\mathrm{\λ}}{d}^2}{{\mathrm{ESun}}_{\mathrm{\λ}}{\cos \mathrm{\θ}}_4}---\left(2\right)$

其中L_λ为传感器在某个波段内接收到的亮度，d为日地距离，ESun_λ为太阳光谱在某波段内的平均辐射照度，θ为太阳入射天顶角；

其中，所述确定每一像元的地表温度T_s包括：

根据公式(4)-(6)计算每一像元的地表温度T_s，

ε＝1.0094+0.047Ln(NDVI)    (4)

式中，ε为自然地表的比辐射率，

$L_{\mathrm{\λ}}\left(T_s\right)=\frac{(L_{\mathrm{\λ}}-L_{\mathrm{\λatm}}\↑)-\mathrm{\τ}(1-\mathrm{\ϵ})L_{\mathrm{\λatm}}\↓}{\mathrm{\τ\ϵ}}---\left(5\right)$

式中，L_λ(T_s)表示温度为T_s的黑体在热红外波段的辐射亮度，L_λatm↓表示大气向下辐射亮度，L_λatm↑表示大气向上辐射亮度，τ为大气在热红外波段的透过率，

$T_s=\frac{K_2}{\ln (\frac{K_1}{L_{\mathrm{\λ}}\left(T_s\right)}+1)}$

其中，K₁＝607.76(Wm^-2μm^-1Sr^-1)，K₂＝1260.56K。

附图说明

图1是Ts/NDVI特征空间图

图2是根据本发明的像元分解组分示意图。

图3是研究区的位置图。

图4a和4b分别是地表温度和地表干燥度指数(TVDI)示意图。

图5NDVI与地表干燥度指数的散点图。

图6a、6b、6c依次是地表裸地盖度、草地盖度、以及林地和灌丛混合类型盖度示意图。

图7a和7b分别是灌丛盖度和林地盖度示意图。

图8是土壤水分与TVDI的相关关系的示意图。

图9是灌丛和草地盖度实测值与估测值比较示意图。

图10是根据本发明的实施例的流程图。

具体实施方式

本发明的各个方面将依据附图1-9结合下述描述得以澄清。

1基于NDVI和TVDI的线性分解模型

1.1 NDVI和地表温度

利用Landsat TM第6波段(10.4～12.5μm)计算地表温度，其中TM是指美国陆地卫星4～5号专题制图仪(Thematic Mapper)所获取的多波段扫描影像。首先TM影像的第3、4、6波段进行辐射定标，把DN值转换为相应的表观辐射亮度(DN指的是图像的灰度值)。采用Landsat用户手册中提供的辐射校正方程(2)：

$L_{\mathrm{\λ}}=\frac{L_{\max \mathrm{\λ}}-L_{\min \mathrm{\λ}}}{\mathrm{QCAL}-{\mathrm{QCAL}}_{\min }}\×(\mathrm{QCAL}-{\mathrm{QCAL}}_{\min })+L_{\min \mathrm{\λ}}---\left(1\right)$

公式(2)中，L_λ表示传感器所接收到的辐射强度，QCAL表示TM数据的像元灰度值(DN值)，QCAL_max为传感器所接收到的最大DN值，QCAL_min为传感器所接收到的最小DN值，Lmax_λ为传感器所接收到的的最大光谱辐射值，Lmin_λ为传感器所接收到的最小光谱辐射值。

在假设地表为朗伯体的前提下，利用公式(2)求得表观反射率ρ_App。

${\mathrm{\ρ}}_{\mathrm{App}}=\frac{{\mathrm{\πL}}_{\mathrm{\λ}}{d}^2}{{\mathrm{ESun}}_{\mathrm{\λ}}{\cos \mathrm{\θ}}_4}---\left(2\right)$

其中L_λ为传感器在某个波段内接收到的亮度，d为日地距离(以日地平均距离为单位)，ESun_λ为太阳光谱在某波段内的平均辐射照度，θ为太阳入射天顶角。

利用大气校正软件6S对表观反射率进行校正，得到TM第3、第4波段的地表反射率，然后计算NDVI(公式3)。

$\mathrm{NDVI}=\frac{R_{\mathrm{NIR}}-R_{\mathrm{RED}}}{R_{\mathrm{NIR}}+R_{\mathrm{RED}}}---\left(3\right)$

R_NIR为TM第4通道近红外波段的光谱反射率，R_RED为TM第3通道红光波段的光谱反射率。NDVI值介于-1～1之间，正值的增加表示绿色植被增长，负值主要为积雪、水体、沙漠以及裸地。

Van通过实地同时测量一系列自然地表的热红外(8～14μm)比辐射率和归一化植被指数NDVI后发现，经过对数转换以后，它们之间具有良好的相关性，其相关系数为0.94。干燥的土壤比辐射率一般低于植被的比辐射率，这与植被叶子中含有较多的水分有关。根据Van的经验公式，NDVI较高的地表比辐射率相应较高，裸土的NDVI低，相应比辐射率也较低。这符合一般的经验，能在一定程度上反映地表比辐射率的变化规律。

ε＝1.0094+0.047Ln(NDVI)    (4)

ε为自然地表的比辐射率。水体在热波段范围内的比辐射率很高，接近于黑体，采用0.9925。

根据(1)式得到TM6表观辐射亮度L_λ，利用(5)式及大气参数，可以得到由T_s决定的辐射亮度L_λ(T_s)：

$L_{\mathrm{\λ}}\left(T_s\right)=\frac{(L_{\mathrm{\λ}}-L_{\mathrm{\λatm}}\↑)-\mathrm{\τ}(1-\mathrm{\ϵ})L_{\mathrm{\λatm}}\↓}{\mathrm{\τ\ϵ}}---\left(5\right)$

L_λ(T_s)表示温度为T_s的黑体在热红外波段的辐射亮度，L_λatm↓表示大气向下辐射亮度，L_λatm↑表示大气向上辐射亮度，τ为大气在热红外波段的透过率。由于没有卫星过境时的同步气象数据，不能很好的模拟当时的大气状况。这里我们参考中纬度夏季标准大气剖面，模拟得到各个大气参数，L_λatm↓为1.68Wm^-2μm^-1Sr^-1，  L_λatm↑为1.74Wm^-2μm^-1Sr^-1，τ为0.77。

利用普朗克公式的反函数，可以根据辐射亮度L_λ(T_s)求得地表真实温度T_s，简化公式如下：

$T_s=\frac{K_2}{\ln (\frac{K_1}{L_{\mathrm{\λ}}\left(T_s\right)}+1)}---\left(6\right)$

其中，K₁＝607.76(Wm^-2μm^-1Sr^-1)，K₂＝1260.56K。

1.2TVDI干燥度指数

植被覆盖度、蒸散量、地表热特性、净辐射、地表粗糙度和风速等因素影响地表温度，地表温度与土壤水分状况之间不存在直接的关系，但土壤水分是影响植被冠层温度的重要因素，从这个意义上，一定植被覆盖条件下的冠层温度能够间接反映土壤供水情况。研究表明，对于水分条件良好的地表，地表温度和NDVI的关系与地表土壤水分更为直接相关。基于Ts/NDVI特征空间，通过地表温度和植被指数散点构成的三角形进行了分析(参见图1)，提出了TVDI指数。

图1中，TVDI为A和B的比值。左侧边代表了不同湿度的裸土的温度变化，Ts和NDVI关系为随着植被绿度的上升，最大地表温度下降。在Ts/NDVI特征空间中定义等值线来代表不同的干旱程度。TVDI计算公式如(7)所示：

$\mathrm{TVDI}=\frac{T_s-T_{S_{\min }}}{T_{S_{\max }}-T_{S_{\min }}}---\left(7\right)$

Ts_min表示最小地表温度，对应的是湿边，Ts_min＝a₁+b₁NDVI，a₁、b₁是湿边拟合方程的系数；Ts是任意像元的地表温度；Ts_max为某一NDVI对应的最高温度，即干边，Ts_max＝a₂+b₂NDVI，a₂、b₂是干边拟合方程的系数。TVDI值为1，是干边(Dry edge)，代表有限的水分供应；TVDI值为0，则是湿边(Wet edge)，具有最大的土壤蒸发蒸腾总量和无限的水分供应。TVDI越大，土壤湿度越低，TVDI越小，土壤湿度越高。估计这些参数要求研究区域的范围足够大，地表覆盖从裸土变化到比较稠密的植被覆盖。

1.3线性分解模型

1.3.1分解方法

土壤水分是干旱、半干旱地区的主要生态限制因子，土壤水分降低导致植物生长量明显减退。可以看出，土壤水分对植被盖度有很大的影响。由于TVDI可以反映土壤的水分状况，因此，可以利用TVDI与NDVI结合估测植被盖度。根据研究区内的植被类型，把研究区内每30×30m范围内的植被盖度分为以下三种类型：(1)土地覆盖类型为裸地，其盖度表示为fc_b；(2)土地覆盖类型为草地，其盖度表示为fc_g；(3)除上述两种覆盖类型外，30×30m地区内的剩余覆盖类型为森林和灌丛混合类型，其盖度表示为fc_f/s，包括两部分：林地的盖度fc_forest和灌丛盖度fc_shrub(图2)。

根据线性分解模型，提出以下的假设：假定混合像元的NDVI和对TVDI的敏感性是该像元内的各类地物组分的线性组合，通过线性分离，可以确定每类地物所占比率。如果一个像元内仅包含一种地物，则称这个像元为典型像元，即端元，对应地物称为典型地物。

根据像元内植被盖度的定义，每个像元内所有组分的盖度之和应为100％。可以得到每30×30m像元内三种类型植被盖度总和为1，因此，可以用以下公式来表示：

fc_b+fc_g+fc_f/s＝1    (8)

将裸地盖度、草地盖度和森林与灌丛混合类型盖度分别作为30m像元内每种组分像元NDVI的权重系数，因此，可以得到公式(9)：

NDVI_bfc_b+NDVI_gfc_g+NDVI_f/sfc_f/s＝NDVI       (9)

公式(9)中，NDVI_b，NDVI_g，NDVI_f/s分别表示每像元内裸地、草地和森林与灌丛混合类型对应的NDVI值。NDVI表明了整个像元的植被生长状况，NDVI_b，NDVI_g，NDVI_f/s分别表明了每像元内裸地、草地和森林与灌丛混合类型的生长状况。

同样，将裸地盖度、草地盖度、森林和灌丛混合类型盖度分别作为像元内每种组分像元NDVI对TVDI敏感性的权重系数，可以得到公式(10)：

α_bfc_b+α_gfc_g+α_f/sfc_f/s＝α    (10)

公式(10)中，α表示每个像元NDVI对TVDI的敏感性，α_b，α_g，α_f/s分别表示每像元内裸地、草地、森林和灌丛混合类型对应的NDVI值对TVDI敏感性。每个像元的NDVI值对TVDI敏感性用公式(11)计算。

$\mathrm{\α}=\mathrm{dNDVI}/d\hat{T}---\left(11\right)$

公式(11)中

是标准化后的TVDI，其计算方法如公式(12)所示。

$\hat{T}=(\mathrm{TVDI}-\stackrel{\‾}{T})/\mathrm{\σ}---\left(12\right)$

公式(12)中，T为多年TVDI时相平均值，σ为TVDI的标准差，其计算方法如公式(13)所示。

$\mathrm{\σ}=\sqrt{\underset{i=1}{\overset{18}{\mathrm{\Σ}}}{({\mathrm{TVDI}}_i-\stackrel{\‾}{T})}^2/(n-1)}---\left(13\right)$

利用一元线性回归拟合模型计算NDVI对TVDI的敏感性，并进行了F检验，选取通过0.05显著水平检验的拟合直线的斜率作为NDVI对TVDI的敏感性。

利用最小二乘法对每一像元进行分解，分解的过程中，线性混合模型受到约束条件的限制。即：0≤fc_b≤1，0≤fc_g≤1，0≤fc_f/s≤1。

1.3.2线性分解模型端元

在线性分解模型中，像元的基本组分的类型应具有代表性，是研究区内多数像元的有效组成成分。而各个类型的端元直接影响线性模型分解的精度。在概念分析的基础上，依据植被类型分布地区植被对土壤水分的敏感性和可以反映不同植被类型生长状况的NDVI值的差异确定端元。一般而言，在同一像元内，林地和灌丛NDVI值大于裸地和草地的NDVI值，草地的NDVI值低于灌丛而高于裸地NDVI值。裸地对气候干旱响应不敏感，草地、林地和灌丛对气候干旱响应敏感，而草地、林地和灌丛相比较，由于林地和灌丛对干旱有的适应能力高于草地，草地对气候干旱响应更为敏感。以表示植被对气候土壤水分敏感性的TVDI为α为自变量，表明植被生长状况的NDVI为因变量绘制散点图，散点图构成一个三角形，三角形端点的坐标为线性分解模型端元值。在像元内三种覆盖类型中，植被对干燥度指数敏感性较小，NDVI值最大的三角形端点为林地和灌丛混合类型的端元(α_f/s，NDVI_f/s)，植被对气候干燥度指数敏感性最大，NDVI值较大的三角形端点为草地的端元(α_g，NDVI_g)，植被对干燥度指数敏感性较小，NDVI值最小的三角形端点为裸地的端元(α_b，NDVI_b)。

1.3.3林地、灌从的植被盖度估算

林地和灌丛混合类型的端元α_f/s，NDVI_f/s，草地的端元α_g，NDVI_g，裸地的端元α_b，NDVI_b确定后，利用最小二乘法，对研究区内每一像元进行线性分解，得到裸地盖度、草地盖度和森林与灌丛混合类型盖度三种覆盖类型的盖度分别为fc_b、fc_g、fc_f/s。由图2可知，fc_b包括林地的盖度fc_forest和灌丛盖度fc_shrub。如果每像元全为林地所覆盖，对应的NDVI为NDVI_f ^*；如果每像元全为灌丛所覆盖，对应的NDVI表示为NDVI_s ^*，可以得到(14)：

${\mathrm{NDVI}}_f^{*}{\mathrm{fc}}_{\mathrm{forest}}=\stackrel{\‾}{{\mathrm{NDVI}}_{f/b}}{\mathrm{fc}}_{f/b}-{\mathrm{NDVI}}_s^{*}({\mathrm{fc}}_{f/b}-{\mathrm{fc}}_{\mathrm{forest}})---\left(14\right)$

公式(14)变换后得到(15)：

${\mathrm{fc}}_{\mathrm{forest}}=\frac{\stackrel{\‾}{{\mathrm{NDVI}}_{f/b}}{\mathrm{fc}}_{f/b}-{\mathrm{NDVI}}_s^{*}{\mathrm{fc}}_{f/b}}{{\mathrm{NDVI}}_f^{*}-{\mathrm{NDVI}}_s^{*}}---\left(15\right)$

通过公式(15)计算得到林地的盖度，必须确定NDVI_f ^*和NDVI_s ^*值。参考研究区的植被类型图，对影像分类，在分类图像中选取一定数量的全为林地或灌丛覆盖“纯“像元，确定与这些像元对应的NDVI值，对这两种类型的纯像元的NDVI分别求平均值。得到的平均值为NDVI_f ^*和NDVI_s ^*值。

公式(15)中，NDVI_f/b可以利用NDVI对TVDI的敏感性之间的关系建立方程，即为：

$\stackrel{\‾}{{\mathrm{NDVI}}_{f/b}}={\mathrm{NDVI}}_{f/b}+\hat{T}*\mathrm{\α}---\left(16\right)$

灌丛的盖度可用公式(17)计算。

fc_shrub＝fc_f/b-fc_forest    (17)

参见图10，根据本发明基于植被-土壤水分响应关系的植被盖度估算方法，(1)计算地表温度：对影像进行辐射校正，然后利用6S模型进行大气校正，计算NDVI，计算比辐射率，然后利用辐射传输方程计算地表温度。(2)TVDI：根据地表温度，计算每个NDVI值对应的最大和最小温度，然后计算干边和湿边的方程，再计算TVDI。(3)NDVI对TVDI响应的敏感性：利用公式(11)-(13)计算。(4)线性分解模型的端元：在植被NDVI对土壤水分敏感性α和NDVI构成的散点图上，确定植被对TVDI敏感性较小，NDVI值最大的三角形端点为林地和灌丛混合类型的端元(α_f/s，NDVI_f/s)，植被对TVDI敏感性最大，NDVI值较大的三角形端点为草地的端元(α_g，NDVI_g)，对TVDI敏感性较小，NDVI值最小的三角形端点为裸地的端元(α_b，NDVI_b)。(5)裸地、草地、灌丛和林地混合类型盖度：在限制条件下，基于最小二乘法，对植被NDVI和对土壤水分的敏感性进行线性分解，得到像元内不同组分的覆盖度。(6)灌丛、林地盖度：根据公式(14)-(17)，对林地和灌丛的盖度进一步分解，得到灌丛和林地盖度。

2线性分解模型应用

2.1研究区

选择中国永定河流域一部分区域作为研究区(图3)。永定河是流经北京市的最大的河流，也是北京市的重要行洪河道之一。自桑干河上源至屈家店河长680km，流域面积47066km²。流域位于欧亚大陆东部中纬度地带，横跨我国南温带半湿润、半干旱气候区和北温带半干旱气候区。属大陆性气候，冬季较长、干燥寒冷、盛行西北风，夏季较短，春秋多风沙，冷暖变化显著。流域多年平均降水量为406mm。降水量年内分配不均，年降水量的77％集中在汛期。

2.2数据及预处理

选择3期TM影像(2006、2004、1997年)和1期ETM+(2000年)作为影像数据源，TM影像包括6个反射率波段，分辨率为30m，1个热波段，分辨率为120m，ETM+影像包括8个波段，热波段的分辨率为60m。每期影像的行列数分别为3365和2426。首先利用1∶100000地形图对TM和ETM+影像进行几何纠正，校正误差小于1个像元。利用(1)式进行辐射纠正。采用6S模型对影像进行大气纠正。参考研究区的植被类型图、DEM及其他辅助资料，对遥感影像进行分类，得到新的植被类型图。野外实测的土壤水分和各种植被类型盖度数据用于估算结果的验证。

2.3结果分析

2.3.1地表温度与TVDI

计算得到2006年地表温度空间变化如图4a所示，地表温度最低值为287.5K，对应于水面，最高值为320.4K，对应于裸地。在植被覆盖较少的地方，地表温度较高，大部分地区的地表温度高于305K。研究区的中部地区地表温度最高。北部和南部山区的地表温度相对较小，小于295K，山区的边缘地区的地面温度较山区高，大部分介于295～300K。

利用NDVI-Ts特征空间中的各像元相应最大和最小陆地表面温度，回归拟合获得2006年旱边和湿边方程，NDVI＝-38.64Ts+324.21，NDVI＝277.34Ts+7.68。NDVI-Ts特征空间中的旱边斜率小于0，表明随着植被覆盖度的增加，陆地表面温度最大值越小；相反湿边的斜率基本上大于，说明随植被覆盖度的增加，陆地表面温度最小值有增加趋势。根据旱边和湿边的方程，计算得到TVDI如图4b所示，TVDI与地表温度具有相似的空间变化格局。因为TVDI值越大，表明越干旱，值越小，表明越湿润，因此，研究区的中北部和南部山区的相对湿润，山区中间地带相对湿润。

2.3.2植被盖度

利用与2006年地表温度和TVDI相同的计算方法，得到1997、2000、2004年研究区的地表温度和TVDI。利用一元线性回归拟合模型计算NDVI对TVDI的敏感性。在2006年植被NDVI对土壤水分敏感性α和NDVI构成的散点图上，确定植被对TVDI敏感性较小，NDVI值最大的三角形端点为林地和灌丛混合类型的端元(α_f/s，NDVI_f/s)，其值为(0.001，0.791)，植被对TVDI敏感性最大，NDVI值较大的三角形端点为草地的端元(α_g，NDVI_g)，其值为(0.368，0.40)，对TVDI敏感性较小，NDVI值最小的三角形端点为裸地的端元(α_b，NDVI_b)，其值为(0.002，0.006)，各端元组分如图5所示。图5中，散点图是由2006年α和NDVI值构成的散点图，α表示每个像元NDVI对TVDI的敏感性。在限制条件下，基于最小二乘法，对2006年植被NDVI和对土壤水分的敏感性进行线性分解，得到像元内不同组分的覆盖度(参见图6)，从图6a中可以看出，在中部干旱地区，裸地盖度达到了0.8以上，而在湿润地区，裸地盖度小于0.2，在较湿润地区，裸地盖度小于0.4。裸地盖度表现出与TVDI相似的空间变化格局。从图6b可以得到，在TVDI大于0.7的地区，草地的盖度小于0.2；在中北部山区和南部山区，草地的盖度达到了0.4～0.6，有些地区甚至达到了0.8；整个研究区中，极少地区的草地盖度介于0.8～1.0之间。根据图6c，大部分地区林地和灌丛混合类型的盖度小于0.2或为0，山区林地和灌丛混合类型的盖度高于0.8，相当一部分地区盖度介于0.6～0.8。

根据公式(14)～(17)，进一步将林地和灌丛混合类型的进行分解。得到林地和灌丛的植被盖度，如图7所示，在灌丛的分布区，盖度主要介于0.4～0.6，部分地区的盖度达到了0.6～0.8，在部分湿润地区，灌丛的盖度高于0.8。在整个研究区，大部分地区林地的盖度为0，在林地分布区，大部分林地盖度介于0.4～0.6。在部分山区，林地的盖度介于0.8～1.0之间。从图6还可以得到，灌丛盖度为0的地区比林地盖度为0的地区分布面积小，表明灌丛比林地的分布更广。

3结果验证

3.1 TVDI与土壤水分关系分析

利用2006年野外实测的土壤水分值与TVDI进行相关性分析，验证TVDI作为干旱指数评价本地区土壤水分的有效性。整个永定河流域实测28个点土壤数据，每个点包含5层(0～10，10～20，20～30，30～40，40～50cm)土壤水分数据，选取本研究区内8个点不同层土壤水分与TVDI的散点图如图8所示，0～10cm和10～20cm土壤层水分与TVDI的相关性均较好，通过了0.05水平t-检验。其它层土壤水分与TVDI的无显著相关性。表明利用TVDI作为表征土壤水分变化趋势的指数，具有一定的合理性。可以利用TVDI在此区域进行土壤水分的相关研究。土壤水分越大，TVDI值越小，两者呈显著的负相关性。

4.2植被盖度验证分析

利用2006年研究区的实测灌丛和草地盖度评价基于线性分解模型估测结果。两种植被类型盖度实测值与估测如图9所示，灌丛的最大误差为0.19，最小误差为0，误差绝对值的均值为0.043。草地的最大误差为0.32，最小误差为0，误差绝对值的均值为0.137。可以看出，灌丛的估测精度高于草地估测精度。对整个估测结果来看，尽管某些点的估测误差较大，但总体估测结果能达到一定的精度，表明该线性分解模型可以运用到本区域进行植被盖度估测。

Claims (3)

1.一种基于植被-土壤水分响应关系来估算一区域内各个像元中的植被盖度的方法，所述植被盖度包括裸地盖度fc_b，草地盖度fc_g、林地和灌丛混合类型盖度fc_f/s，该方法包括步骤：

(1)通过以下步骤来确定所述区域的每一像元的归一化差值植被指数NDVI和干燥度指数TVDI值：

利用辐射校正参数和大气校正参数以及大气校正软件6S对表观反射率进行校正，得到TM第3、第4波段的光谱反射率；

利用所述TM第3、第4波段的光谱反射率并根据公式(3)确定每一像元的NDVI值，

$\mathrm{NDVI}=\frac{R_{\mathrm{NIR}}-R_{\mathrm{RED}}}{R_{\mathrm{NIR}}+R_{\mathrm{RED}}}---\left(3\right)$

式中，R_NIR为TM第4通道近红外波段的光谱反射率，R_RED为TM第3通道红光波段的光谱反射率；

利用中纬度夏季标准大气剖面，模拟得到与每一像元对应的大气参数：大气向下辐射亮度L_λatm↓为1.68Wm^-2μm^-1Sr^-1，大气向上辐射亮度L_λatm↑为1.74Wm^-2μm^-1Sr^-1，大气在热红外波段的透过率τ为0.77；

利用所述大气向下辐射亮度L_λatm↓、大气向上辐射亮度L_λatm↑以及大气在热红外波段的透过率τ并根据公式(4)-(6)确定每一像元的地表温度T_s，

ε＝1.0094+0.047Ln(NDVI)         (4)

$L_{\mathrm{\λ}}\left(T_s\right)=\frac{(L_{\mathrm{\λ}}-L_{\mathrm{\λatm}}\↑)-\mathrm{\τ}(1-\mathrm{\ϵ})L_{\mathrm{\λatm}}\↓}{\mathrm{\τ\ϵ}}---\left(5\right)$

$T_s=\frac{K_2}{\ln (\frac{K_1}{L_{\mathrm{\λ}}\left(T_s\right)}+1)}---\left(6\right)$

式中，ε为自然地表的比辐射率，L_λ(T_s)表示温度为T_s的黑体在热红外波段的辐射亮度，K₁＝607.76(Wm^-2μm^-1Sr^-1)，K₂＝1260.56K；

根据所确定的地表温度T_s并依据公式(7)计算每一像元的TVDI值，

$\mathrm{TVDI}=\frac{T_s-T_{S_{\min }}}{T_{S_{\max }}-T_{S_{\min }}}---\left(7\right)$

其中，Ts_min表示最小地表温度，对应的是湿边，Ts_min＝a₁+b₁NDVI，a₁、b₁是湿边拟合方程的系数；Ts_max为某一NDVI对应的最高温度，即干边，Ts_max＝a₂+b₂NDVI，a₂、b₂是干边拟合方程的系数；

(2)通过公式(11)-(13)确定所述区域的每个像元的NDVI值对TVDI敏感性α值，

$\mathrm{\α}=\mathrm{dNDVI}/d\hat{T}---\left(11\right)$

$\hat{T}=(\mathrm{TVDI}-\stackrel{\‾}{T})/\mathrm{\σ}---\left(12\right)$

$\mathrm{\σ}=\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{TVDI}}_i-\stackrel{\‾}{T})}^2/(n-1)}---\left(13\right)$

式中，

为所述区域的每一个像元的多年TVDI时相平均值，σ为TVDI的标准差；

其中，每一个像元的植被盖度满足公式(8)-(10)定义的线性分解模型，

fc_b+fc_g+fc_f/s＝1            (8)

NDVI_bfc_b+NDVI_gfc_g+NDVI_f/sfc_f/s＝NDVI    (9)

α_bfc_b+α_gfc_g+α_f/sfc_f/s＝α            (10)

式中，NDVI_b，NDVI_g，NDVI_f/s分别表示每像元内裸地、草地和森林与灌丛混合类型对应的NDVI值，α_b，α_g，α_f/s分别表示每像元内裸地、草地、森林和灌丛混合类型对应的NDVI值对TVDI敏感性；

(3)基于所述区域中的每一个像元的α值和NDVI值构建所述区域的三角形散点图，其中所述α值为自变量，NDVI值为因变量；

(4)基于三角形散点图确定三角形的每一个端点的端元的α值和NDVI值；

(5)利用最小二乘法，基于所确定的每一个三角形端点的端元的α值和NDVI值并且在约束条件0≤fc_b≤1，0≤fc_g≤1，0≤fc_f/s≤1下，对区域内每一像元进行线性分解以获得区域内每一像元的fc_b、fc_g、fc_f/s；

(6)通过以下步骤对所述林地和灌丛混合类型盖度fc_f/s进一步处理以得到每一像元内林地盖度fc_forest和林地盖度fc_shrub：

根据所述区域的植被类型图，从该植被类型图的分类图像中选取一定数量的全为林地覆盖的纯像元，并且确定与这些全为林地覆盖的纯像元对应的NDVI值并且求取这些NDVI值的平均值，作为全为林地覆盖的纯像元的

值；

根据所述区域的植被类型图，从该植被类型图的分类图像中选取一定数量的全为灌丛覆盖的纯像元，确定与这些全为灌丛覆盖的纯像元对应的NDVI值并且求取这些NDVI值的平均值，作为全为灌丛覆盖的纯像元的

值；

根据公式(15)-(17)确定每一像元的林地盖度fc_forest和灌丛盖度fc_shrub，

${\mathrm{fc}}_{\mathrm{forest}}=\frac{\stackrel{\‾}{{\mathrm{NDVI}}_{f/s}}{\mathrm{fc}}_{f/s}-{\mathrm{NDVI}}_s^{*}{\mathrm{fc}}_{f/s}}{{\mathrm{NDVI}}_f^{*}-{\mathrm{NDVI}}_s^{*}}---\left(15\right)$

$\stackrel{\‾}{{\mathrm{NDVI}}_{f/s}}={\mathrm{NDVI}}_{f/s}+\hat{T}*\mathrm{\α}---\left(16\right)$

fc_shrub＝fc_f/s-fc_forest     (17)。

2.如权利要求1的所述基于植被-土壤水分响应关系来估算一区域内各个像元中的植被盖度的方法，其中，在基于三角形散点图确定三角形的每一个端点的端元的α值和NDVI值的步骤中，植被对干燥度指数敏感性较小、NDVI值最大的三角形端点为具有α_f/s，NDVI_f/s的林地和灌丛混合类型的端元，植被对气候干燥度指数敏感性最大、NDVI值较大的三角形端点为具有α_g，NDVI_g的草地的端元，植被对干燥度指数敏感性较小，NDVI值最小的三角形端点为具有α_b，NDVI_b的裸地的端元，

所述利用最小二乘法，基于所确定的每一个三角形端点的端元的α值和NDVI值并且在约束条件0≤fc_b≤1，0≤fc_g≤1，0≤fc_f/s≤1下，对区域内每一像元进行线性分解以获得区域内每一像元的fc_b、fc_g、fc_f/s的步骤包括：

将所确定的每一像元的α值，NDVI值，以及所确定的α_f/s，NDVI_f/s，α_g，NDVI_g，α_b，NDVI_b的裸地的端元代入公式(8)-(10)中，调用MatLab软件中的最小二乘函数对公式(8)-(10)进行线性分解以获得每一像元的裸地盖度fc_b，草地盖度fc_g以及林地和灌丛混合类型盖度fc_f/s。

3.如权利要求1的所述基于植被-土壤水分响应关系来估算研究区内各个像元中的植被盖度的方法，其中所述辐射校正参数和大气校正参数如下确定：

利用公式(1)计算辐射校正参数L_λ，

$L_{\mathrm{\λ}}=\frac{L_{\max \mathrm{\λ}}-L_{\min \mathrm{\λ}}}{\mathrm{QCAL}-{\mathrm{QCAL}}_{\min }}\×(\mathrm{QCAL}-{\mathrm{QCAL}}_{\min })+L_{\min \mathrm{\λ}}---\left(1\right)$

式中，L_λ表示传感器所接收到的辐射强度，QCAL表示TM数据的像元灰度值DN值，QCAL_max为传感器所接收到的最大DN值，QCAL_min为传感器所接收到的最小DN值，Lmax_λ为传感器所接收到的最大光谱辐射值，Lmin_λ为传感器所接收到的最小光谱辐射值；

利用公式(2)计算大气校正参数，

${\mathrm{\ρ}}_{\mathrm{App}}=\frac{\mathrm{\π}{L}_{\mathrm{\λ}}{d}^2}{{\mathrm{ESun}}_{\mathrm{\λ}}{\cos \mathrm{\θ}}_4}---\left(2\right)$

其中L_λ为传感器在某个波段内接收到的亮度，d为日地距离，ESun_λ为太阳光谱在某波段内的平均辐射照度，θ为太阳入射天顶角。

Images