基于NDVI的喀斯特水资源遥感定量监测及分析——以贵州省为例

贺中华，陈晓翔，梁 虹，黄法苏，赵 芳

（1.中山大学 地理科学与规划学院，广州510275；2.贵州师范大学 地理与环境科学学院，贵阳550001；3.贵州省水文水资源局，贵阳550002；4.贵州省贵阳市白云区职业技术学校，贵阳550014）

喀斯特无论在世界上还是我国，都是一类脆弱的生态环境，已引起国内外学术界的深切关注。喀斯特流域是具有特殊的双重含水介质，特殊的地表、地下双重分水岭，独特的地貌—水系结构的地域综合体［1］。喀斯特水资源是大气降水在喀斯特流域下垫面再分配的表现，根据喀斯特流域的特征，其水资源可分为地表水资源和地下水资源。影响喀斯特水资源的因素很多，除土地利用类型、岩组类型、地貌类型外，其植被类型也不容忽视。例如高大的乔木和低矮的灌木对减小降水对地表的冲击力不同，不同森林植被通过影响降水在地表的侧向流速、降水在流域下垫面滞流的时间及下渗率，进而影响喀斯特流域的水资源量。目前，衡量流域植被好坏的一个重要指标为植被指数。植被指数（VI）是对地表植被活动的简单、有效和经验的度量，在一定程度上反映流域下垫面的赋水信息。经过近20a的发展，植被指数已有几十种，其中归一化植指数（NDVI）被广泛地应用在全球与区域土地覆盖、植被分类和环境变化等研究领域［2－9］。对于喀斯特水资源的研究，课题组曾做过相关工作［10－11］，而基于 NDVI的喀斯特水资源的定量研究，无论在国内还是国外，未曾见有研究报道。本文在贵州省内选取20个具有连续5a观测水文数据和遥感资料的典型喀斯特流域，利用遥感技术，从TM影像中提取喀斯特流域NDVI，利用现代数学方法，探讨喀斯特流域水资源与NDVI的关系，建立水资源监测、预测模型，并利用5个研究样区进行模型检验。

1 数据选取及处理

1.1 水文数据选取

根据贵州省水文总站整编的《贵州省历年各月平均流量统计资料》以及贵州省水文水资源局整编的《贵州省水资源公报》，选其中都处于相同的气候带的20个水文断面，时间从2005—2010年，流域面积以中小流域为主，目的是为了保证流域下垫面的地质条件能尽可能相同或相近。各流域的9月平均径流深见表3。

1.2 遥感数据选取

数据选用TM影像的2005—2010年，成像时间分别为每年的9月，保证降雨对流域赋水影响较小，保证每个时段研究样区云量小于30%。

1.3 遥感影像处理及计算

1.3.1 遥感影像预处理

（1）大气校正。目前，大气校正的方法有很多，其中，大气辐射传输模型是大气校正中精度较高的方法。它是利用电磁波在大气中的辐射传输原理建立的模型来对遥感图像进行大气校正的。本研究采用的FLAASH模型是改进的MORTRAN模型，它不仅可以对高光谱数据进行大气校正，而且还可以对多光谱数据如 Landsat，SPOT，AVHRR，MERIS，IRS和ASTER等数据进行大气校正。

（2）几何校正。几何校正包括图像对地形图和图像对图像的配准。地形图是国家基础地理信息1∶25万数据，其坐标为地理坐标、采用克拉索夫斯基椭球。影像配准利用多项式中4项式，控制点选择40个左右，配准精度在5个像素内；精度图像对图像配准保证在0.3个像素，为了保证光谱信息，重采样时选取最近邻法。

1.3.2 表观反射率计算

（1）光谱辐射亮度的计算

如果没有定标参数Gain和Bias的资料，某一波段的L可以根据式（2）计算。

式中：QCAL——某一像元的 DN 值，即 QCAL＝DN；QCALmax——像 元 可 以 取 的 最 大 值 255；QCALmin——像元可以取的最小值。对于Landsat—7来说，式（2）［12－13］可改为式（3）（QCALmin＝1）。

（2）表观反射率的计算［13－15］。

式中：ρ——大气层顶（TOA）表观反射率（无量纲）；π——常量（球面度sr）；L——大气层顶进入卫星传感器；D——日地之间距离。根据表1，可以推算全年任何一天的日地距离；ESUN——大气层顶的平均太阳光谱辐照度，根据表2可查得［12］。

θ为太阳的天顶角，地面站提供的头文件给出的是太阳高度角，因此θ＝90°－β。另外，可以用式（5）直接求取cosθ［16］。

式中：φ——地理纬度；δ——太阳赤纬；h——太阳的时角。

1.3.3 NDVI指数的计算 根据植被的反射光谱特征，利用红波段、近红外波段的反射率和其他因子及其组合所获得的植被指数来提取植被信息，且这些波段常包含90%以上有关植被的信息。归一化植被指数是广泛使用的一种植被指数，由Rouse等人提出［17］。

式中：NIR——近红外通道反射率；R——红色通道的反射率。其中，－1≤NDVI≤1，负值表示地面覆盖为云、水、雪等，对可见光高反射；0表示有岩石或裸土等，NIR和R近似相等；正值表示有植被覆盖，且随覆盖度增大而增大。

表1 随时间变化日地距离（天文单位，D）

表2 Landsat－7和Landsat－5的大气层顶平均太阳光谱照度 W／（m2·μm）

本研究选用Landsat—7数据，首先，选用公式（2）计算光谱辐射亮度；其次，选用公式（4）、公式（5），并根据表1、表2计算表观反射率；再次，利用公式（6）计算光谱辐射亮度的LNDVI、表观反射率的ρNDVI，得表3。

表3 喀斯特流域研究样区水文数据及NDVI

2 遥感模型建立

2.1 模型原理

假定喀斯特流域某水文断面观测值Y和该流域植被指数X之间关系可用如下模型表示［18］：

其中，b0，b1，b2，b3是未知因素参数；ε—N（0，σ2）随机变量。为评价回归方程的精度，需对其进行显著性检验，通常用F检验。

因F服从自由度为（m，n－m－1）的F分布，对于指定的α，由F分布表可查Fα（m，n－m－1），如F＞Fα（m，n－m－1），则认为回归模型适合该组资料称它为显著的，否则称为不显著即不能使用。

2.2 模型建立

首先，根据表3，借助SPSS和Matlab统计软件，利用公式（8）计算喀斯特流域水资源与NDVI的相关关系，得到表4；其次，利用公式（7）建立喀斯特水资源监测、预测模型，模型系数如表5所示。图1，图2表示研究样区水资源径流深与NDVI的拟合效果。

表4 径流深与NDVI相关系数矩阵

（1）从表4可知，喀斯特水资源与其植被指数相关性都很高，尤其是地物光谱辐射亮度的归一化植被指数，高达0.857；其次，水资源与地物表观反射率的归一化植被指数的相关性也很高（0.652）；另外，反映喀斯特植被覆盖率的NDVI之间的相关性也很好（0.866）。

（2）从表5可知，分别由地物光谱辐射亮度的归一化植被指数、地物表观反射率的归一化植被指数来拟合喀斯特水资源，其拟合的效果很好，如图1—2所示，拟合度都很高，尤其是地物表观反射率的归一化植被指数来拟合喀斯特水资源，其值达0.971；利用公式（9）对拟合的效果进行F检验，其F的最大值达176.832，最小值为69.815，均大于给定的临界值5.29，说明由此建立的模型高度显著。

表5 模型系数表

（3）根据公式（7），利用表5，其喀斯特流域水资源监测、预测模型可表达为：

图1 喀斯特流域研究样区水资源径流深与LNDVI拟合效果

综上所述，在喀斯特地区，由于地表崎岖，地下洞隙纵横交错，水文动态变化剧烈，地表水渗漏严重，地下持水保水能力差；土层薄、肥力低、植被生长困难，水土流失严重，形成了独特的、脆弱的喀斯特自然环境，严重地制约喀斯特流域的持水、供水能力。喀斯特流域具有特殊的双重含水介质，形成独特的地表－地下水系结构，因此，喀斯特流域与正常流域特别是湿润地区常态流域相比，其流域水资源的形成机制、空间分布规律具有一定的特殊性。流域植被类型及覆盖率将直接影响降雨在喀斯特流域入渗及径流，即影响降雨在流域空间的再分配，因此，喀斯特流域植被指数是喀斯特流域赋水状况的重要性指标。

图2 喀斯特流域研究样区水资源径流深与ρNDVI拟合效果

3 遥感模型检验

为了评定监测、预测模型的精度，任选5个喀斯特流域作为研究样区，按上述方法对研究样区遥感影像进行处理，分别提取LNDVI、ρNDVI，如表6所示，分别代入模型（10），（11）进行计算，并与实测数值对比见表6。通过计算比较得出，模型（10），（11）相对误差值都比较小，说明用这两个模型对喀斯特流域水资源进行监测、预测，效果是比较理想的，且模型（11），效果更好，精度更高。从理论上分析，原始遥感影像的DN是未经过任何校正，包括辐射定标校正，只是进入传感器中的辐射能的一种数字转换形式，不能本质地反映地物的辐射特性。L和ρ都经过了辐射定标校正，但是，当ρ再经过大气校正后，它就是地物的反射率，能本质地反映地物的辐射特性，因此，由它构成的NDVI植被指数最接近地物的NDVI。

表6 模型检验结果

4 结论与分析

（1）喀斯特具有特殊的流域下垫面介质结构，其流域产、汇机制复杂，流域赋水影响因素多样，流域植被覆盖率起到重要的作用。

（2）利用地物光谱辐射亮度的归一化植被植指数（LNDVI）和地物表观反射率的归一化植被指数（ρNDVI）对喀斯特流域水资源进行监测、预测效果很好，尤其是利用ρNDVI进行监测、预测，精度更高。

（3）适合于喀斯特流域水资源监测、预测数学模型是：

通过方差分析和样区检验，得出很好的预测效果。

［1］ 杨明德，谭明，梁虹.喀斯特流域水文地貌系统［M］.北京：地质出版社，1998.

［2］ 何隆华，储开华，肖向明.Vegetation图像植被指数与实测水稻叶面积指数的关系［J］.遥感学报，2004，8（6）：672－676.

［3］ 王福民，黄敬峰，唐延林，等.采用不同光谱波段宽度的归一化植被指数估算水稻叶面积指数［J］.应用生态学报，2007，18（11）：2444－2450.

［4］ 杨曦，武建军，闫峰，等.基于地表温度—植被指数特征空间的区域土壤干湿状况［J］.生态学报，2009，29（3）：1205－1216.

［5］ 李玉霞，杨武年，童玲，等.基于光谱指数法的植被含水量遥感定量监测及分析［J］.光学学报，2009，29（5）：1403－1407.

［6］ Fabio Maselli1，Antonio Di Gregorio，Valerio Capecchi1 et al.Enrichment of land－cover polygons with eco－climatic information derived from MODIS NDVI imagery［J］.Journal of Biogeography，2009，36：639－650.

［7］ Karnieli A，Nurit A，Rachel T，et al.Use of NDVI and land surface temperature for drought assessment［J］.Journal of Climate，2010，23（3）：618－633.

［8］ Pu Ruiliang，Gong Peng，Tian Yong，et al.Using classification and NDVI differencing methods for monitoring sparse vegetation coverage：a case study of saltcedar in Nevada，USA［J］.International Journal of Remote Sensing，2008，29（14）：3987－4011.

［9］Thenkabail P S.Inter－sensor relationships between IKO－NOS and Landsat－7ETM＋NDVI data in three ecoregions of Africa［J］.International Journal of Remote Sensing，2004，25（2）：389－408.

［10］ 贺中华，梁虹，黄法苏，等.基于RS流域枯水资源的判读识别［J］.贵州师范大学学报：自然科学版，2004，22（2）：36－39.

［11］ 贺中华，杨胜天，梁虹，等.基于GIS和RS的喀斯特流域枯水资源影响因素识别：以贵州省为例［J］.中国岩溶，2004，23（1）：48－55.

［12］ Chander G，Markham B L.Revised Landsat－5TM radiometric calibration procedures and post－calibration dynamic ranges［J］.Geoscience and Remote Sensing，2003，41（1）：2674－2677.

［13］ Markham B L，Barker J L.Thematic Mapper bandpass solar exoatmospheric irradianees［J］.International Journal of Remote Sensing，1987，8（3），513－523.

［14］ Roderick M，Smith R，Lodwick G.Calibrating longterm AVHRR－derived NDVI imagery［J］.Remote Sensing of Environment，1996，58（1）：1－12.

［15］ Huete A.Overview of the radiometric and biophysical performance of the MODIS vegetation indices［J］.Remote Sensing of Environment，2002，83（1）：195－213.

［16］ LüS H.Physical Foundation of Remote Sensing［M］.Beijing：Commercial Press，1981：206－234.

［17］ Rouse J W，Haas R H，Schell J A，et al.Monitoring vegetation systems in the great plains with ERTS［J］.Third ERTS Symposium，NASA，1973.

［18］ 张超，杨秉赓.计量地理学基础［M］.北京：高等教育出版社，1985.