典型喀斯特地区土壤水分遥感反演——以涟江流域为例

梁建方，周秋文，韦小茶，马龙生，颜红

(贵州师范大学地理与环境科学学院，贵州贵阳550001)

土壤是重要的自然资源，土壤水分作为陆地和大气能量交换过程中的重要因子，对陆表蒸散、水的运动和碳循环具有控制作用，对自然环境和人类社会的可持续发展有着重要影响[1-3]。获取区域乃至全球尺度的土壤水分空间分布信息，是研究全球范围内水循环、能量平衡和气象预报的先决条件，对改善区域及全球气候模式预报结果、进行农作物旱情监测及估产、自然和生态环境问题的研究起重要作用[4-5]。传统土壤水分测量方法只能得到有限测量点数据，仅能反应测量点周围有限区域的土壤水分情况。而通过星载遥感卫星传感器提供的监测手段，能获取大区域的土壤水分信息，同时具有时效快、动态性较强及时间空间分辨率高等优点，能为区域尺度土壤水分信息的获取提供良好手段[6-8]。

光学遥感监测土壤水分的方法主要有热惯量法、植被指数法等，但热惯量法仅适用于裸地及低值被覆盖区域，而植被指数法仅适用于植被覆盖较高区域[1,8]。Price和Carlson研究发现植被指数和地表温度之间存在负相关关系，且以植被指数和地表温度为横纵坐标生成的散点图呈三角形状[9-10]。Sandholt等[11]进一步研究植被指数与地表温度的关系，提出利用温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)来监测和估算土壤水分的方法。国内外学者基于该方法进行区域地表水分状况等开展大量研究，并取得了丰硕成果[12]。王秀君等[13]基于Landsat TM数据反演了河北省黄骅市土壤水分，结果表明TVDI与土壤湿度显著相关；曹雷等[14]利用TVDI在艾比湖地区进行试验，发现TVDI能有效反演区域土壤水分，且精度较高；陈斌等[15]以内蒙古锡林郭勒盟地区为研究对象，发现TVDI可以用来监测旱情变化及评价草原干旱状况；沙莎等[16]利用历史遥感数据建立的特征空间，发现TVDI能够较好地指示每年土壤相对湿度RSM的时空分布特征；吴孟泉、张喆等[17-18]研究均发现TVDI能够很好地用于山区干旱预警与监测。

喀斯特地区生态系统敏感度高, 灾变承受能力低, 环境容量小, 因而成为典型的生态脆弱区[19]。其土层浅薄、土壤贫瘠、成土速率低、水文过程变化迅速、地表土壤水受特殊地质地貌及水文条件等自然因素影响较大[20-24]。以贵州省为中心的中国西南喀斯特地区是世界上面积最大的喀斯特连续带，该地区地表破碎、地形复杂，由于森林砍伐、过度放牧、农业耕种和地下水过度开采等行为，地表植被的生长和人类生产生活深受土壤水分贫瘠影响[25-30]。郭巍等[26]采用TVDI对云南省2010年春季旱情进行监测，揭示了云南省入春以来的旱情动态变化及旱情分布特点；王海等[27]运用云南省境内的TVDI与昆明3个站点的土壤湿度做线性相关分析表明，TVDI对云南境内的干旱监测具有较好的适用性和可靠性；康为民等[28]利用同一时期TVDI与当地气象站点的土壤湿度数据进行定量验证，发现TVDI与土壤湿度显著相关，适用于大区域复杂地形的干旱检测与预警。

然而已有研究大都基于MODIS数据与Landsat TM/ETM+数据来进行干旱监测研究，而利用Landsat8 OLI_TIRS数据进行研究干旱监测研究较少。喀斯特地貌的典型特征是地表破碎、地形复杂，土壤水分分布情况随降雨及地表情况的不同而表现出不同区域特征。因此，本文基于Landsat8 OLI_TIRS数据，选取贵州省典型喀斯特地区的涟江流域为例，对流域的土壤水分进行遥感反演，并实现干旱时空分布情况的动态监测，同时结合该区野外实测土壤含水量数据进行验证分析。研究结果可为相关部门制定合理有效的抗旱、防旱措施提供科学依据，对喀斯特地区水土保持政策的制定、石漠化治理措施实施具有非常重要的指导作用，是喀斯特地区自然、经济、社会可持续发展的重要需求。

1 研究区概况

涟江流域(图1)位于贵州省西南部(106°14′～106°51′E，25°32′～26°26′N)，面积约3 006.67 km2，流域海拔高度在517～1 650 m之间，地势北高南低，中部平缓，其境内主要分布碳酸盐岩，喀斯特地貌典型发育。属亚热带季风气候，降水丰富，但降水量年际、季节性分布不均，而地表和地下水系构成的二元结构往往会导致降水渗漏，土壤储水能力低，加上不合理的人类活动造成植被不断遭受破坏，全区一年四季均存在轻重程度不相同的旱情。

图1 涟江流域地理位置及高程

2 数据来源与方法

2.1 数据来源及预处理

遥感数据来源于中国科学院计算机网络信息中心地理空间数据云平台(http://www.gscloud.cn)，获取2013年9月29日、2014年1月19日与2017年4月1日的Landsat8 OLI_TIRS数据，以及ASTER GDEM高程数据，空间分辨率为30 m。利用ENVI 5.1软件对遥感影像进行预处理，主要包括辐射定标、裁剪和大气校正等。选取涟江子流域(扬眉河流域)中具有代表性的林地、灌木林地、草地、耕地，于2017年4月2日使用TDR-300土壤水分仪测量得到24个样点表层(0～10 cm)土壤体积含水量，并用手持GPS获取采样点的经纬度信息。

2.2 温度植被干旱指数估算

Price和Carlson等研究发现以NDVI与Ts为横纵坐标生成散点图构成的空间关系为三角形关系。Sandholt等对NDVI-LST特征空间进行简化，提出温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)概念，计算公式如下：

(1)

式中TVDI——温度植被干旱指数；Ts——地表真实温度；Tsmax——某一NDVI对应最大地表温度，即干边；Tsmin——地表最低温度，即湿边。干边与湿边方程分别为：Tsmax=a·NDVI+b,Tsmin=c·NDVI+d，其中a、b、c、d为干、湿边拟合方程系数。

模型中地表真实温度Ts由辐射传导方程法获得，具体步骤如下。

a) 计算植被指数NDVI。NDVI由遥感影像中，近红外波段的反射值与红波段的反射值之差比上两者之和，其计算公式为：

(2)

式中NIR——近红外波段像元亮度值；R——红波段像元亮度值。

b) 计算地表比辐射率ε。本文采用Van经验公式[19]计算地表比辐射率：

ε=1.0094 + 0.047·ln·(NDVI)

(3)

式中，NDVI值在0.157～0.727之间，NDVI小于0的像元只要是水体，其体表比辐射率近似0.995；城镇很多建筑物像元NDVI位于0～0.157之间，其地表比辐射率可设置为0.923；而对于NDVI大于0.727像元可看作植被完全覆盖，地表比辐射率为0.986。

c) 计算相同温度下黑体的辐射亮度值B(TS)。卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑；地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。其辐射传输方程为：

Lλ=[ε·B(TS) + (1-ε)L↓]·τ+L↑

(4)

根据普朗克定律得到的黑体热辐射亮度B(TS)，其公式为：

B(TS) = [Lλ-L↑-τ·(1-ε)L↓]/(τ·ε)

(5)

式中Lλ——TIRS影像第十波段大气校正后的辐射定标值；ε——地表辐射率；TS——地表真实温度；τ——大气在热红外波段透过率。

d) 计算地表真实温度TS。在获取B(TS)后，根据普朗克公式的反函数，求得地表真实温度TS，其公式为：

TS=K2/ln(K1/B(TS)+ 1)

(6)

式中B(TS)——温度为TS黑体热辐射亮度；K1、K2——为传感器定标常数。

2.3 表层土壤水计算

利用野外表层土壤水分(0～10 cm)采样数据，建立土壤水分与温度植被干旱指数之间的线性相关模型：

W=a·TVDI+b

(7)

式中W——表层土壤水分；TVDI——温度植被干旱指数；a、b——线性系数。

3 结果与分析

3.1 干、湿边拟合

由于NDVI值小于0的像元主要是水体和云，湿度接近100%，NDVI值在0～0.157时的像元主要是城市建筑物与水泥地等。为避免水体、云及建筑物对特征空间的影响，在干湿边拟合时不考虑此类像元。因此，选取NDVI值大于0.157的像元，利用IDL编程实现相同NDVI值对应的TS的最大值和最小值提取，步长值为0.002，构成NDVI-TS特征空间的形状近似呈三角形(图2)。从3个季节的拟合结果可以看出，干边与湿边的决定系数R2分别在0.9和0.5以上，拟合精度较高。且NDVI与地表温度LST近似呈线性关系，随着NDVI值的增大，地表温度最大值LSTmax逐渐减小，最小值LSTmin逐渐增大。

a) 秋季

b) 冬季

c) 春季图2 Ts-NDVI特征空间

3.2 TVDI干旱等级划分

为从宏观了解涟江流域干旱情况，根据已拟合的干、湿边方程，计算研究区内各个像元TVDI，根据TVDI灰度值范围进行干旱等级划分指标[21]，将干旱情况划分为五级：湿润(0～0.2)、正常(0.2～0.4)、轻旱(0.4～0.6)、干旱(0.6～0.8)、重旱(0.8～1.0)。研究结果表明：研究区3个季节干旱情况在时间分布上总体一致，都为正常>轻旱>湿润>干旱>重旱。3个季节中干旱等级主要以湿润和正常为主，面积分别为1 876.911、1 983.149、1 905.917 km2，均占流域面积的62%以上(表1、图3)；全区3个季节均有轻旱分布，面积所占比例达30%以上，其中春季>秋季>冬季；干旱分布较少，具体分布情况为秋季>冬季>春季，研究区几乎无重旱。而研究区3个季节的旱情在空间分布上存在一定差异，全区3个季节轻旱分布较广，分布区域基本一致；秋季干旱主要分布在流域中部地势较为平坦的耕作区域，南部有较少分布；在冬季和春季，干旱则主要分布在流域中惠水县城至花溪大学城一带，该地区城镇化程度高，周围有大量裸地，可能是导致其干旱的原因。

3.3 土壤含水量反演

将GPS记录实测土壤水分采样点经纬度信息导入ArcGIS，提取采样点对应2017年4月1日TVDI值，与实测土壤体积含水量进行最小二乘法回归分析，得到两者关系(图4)及回归方程：y= -149.47TVDI+97.255(R2= 0.7572)，且拟合结果通过t检验，表明两者之间具有良好的线性关系。把该线性模型输入ENVI 5.1软件中进行反演，即可得到土壤水分分布(图5)。由图3、5可以看出，土壤水分与干旱情况分布基本一致，且土壤水分与TVDI呈现出明显负相关关系，土壤水分随着TVDI的增大而逐渐减小，土壤含水量越低，干旱情况则越严重，这说明了TVDI能够较好地反映研究区地表土壤水分状况。

表1 土壤干旱等级、面积及比例

a) 秋季

b) 冬季

c) 春季图3 不同季节土壤干旱等级空间分布

图4 土壤表层含水量和TVDI关系

a) 秋季图5 不同季节下土壤水分含量分布

b) 冬季

c) 春季续图5 不同季节下土壤水分含量分布

3.4 不同坡度土壤水分季节变化分析

基于DEM高程数据，利用ArcGIS空间分析得到研究区的坡度因子，根据1984年中国农业区划委员会颁发的《土地利用现状调查技术规程》把坡度分为五级，即Ⅰ(≤2°)、Ⅱ(2°～6°)、Ⅲ(6°～15°)、Ⅳ(15°～25°)、V(>25°)，通过分区统计，得到不同坡度等级下土壤水分柱状图(图6)。从图6可知，冬季和春季土壤水分均值随着坡度增大而增加，秋季则是先随着坡度增加，再到25°以上陡坡区域有所降低。这是由于秋季降水较少，而地表温度较高，且在坡度较陡、植被覆盖较差或地表裸露率高的区域，土壤水分相对较低。在第Ⅰ、Ⅱ、Ⅲ坡度等级下，3个季节的土壤水分均值中，冬季相对秋季和春季较高，并在第Ⅰ、Ⅱ坡度等级差异较为明显，分析其原因可能是该坡度等级下区域主要为耕地和建设用地，秋季和春季地表温度相对冬季较高，土壤水分则较少。

图6 不同季节坡度等级下土壤水分含量

3.5 不同坡向土壤水分季节变化分析

通过研究区DEM数据进行坡向提取，选取0°为正北方向，沿顺时针旋转，将坡向划分为9种类型，即FLAT(-1)、N(0°～22.5°，337.5°～360°)、NE(22.5°～67.5°)、E(67.5°～112.5°)、SE(112.5°～157.5°)、S(157.5°～202.5°)、SW(202.5°～247.5°)、W(247.5°～292.5°)、NW(292.5°～337.5°)。通过分区统计得到不同坡向等级土壤水分折线(图7)。从图7可以看出，坡向对3个季节的土壤水分影响基本一致，土壤水分均值在不同坡向上存在着差异性变化规律。土壤水分均值按坡向等级划分，从左到右先增加再减少，然后又先减少再增大，其在北坡和西北坡相对较高，出现峰值，而在无坡向和东南坡相对较低，出现低谷，可以看出研究区土壤水分含量阴坡比阳坡和无坡向多。相同坡向土壤水分均值相对冬季来看，春季与秋季差值相对较小，冬季则较大，产生此现象原因可能是冬季气温较低，地表蒸散发较小；秋季虽植被覆盖较好，但地表温度较高，蒸散发相对较大；春季与冬季相比，春季地表温度上升，蒸散发也随之增大。

图7 不同季节下坡向等级土壤水分含量

3.6 不同土地利用类型土壤水分季节变化分析

利用ENVI 5.1软件对2017年4月的Landsat8 OLI_TIRS数据进行假彩色合成，然后采取监督分类法对研究区土地利用类型进行分类，主要分为耕地、林地、灌木林地、草地、建设用地、水域和其他用地7个类型(图8)。选取前5个土地利用类型对3个季节土壤水分进行统计分析，得到不同土地利用类型土壤水分均值折线图(图9)。从图9可以看出，林地在3个季节中的土壤水分含量最高，灌木林地次之，其他3个地类则存在差异。而3个季节中不同土地利用类型的土壤水分分布在春季和秋季是一致的，主要表现为林地>灌木林地>草地>耕地>其他用地，冬季则为林地>灌木林地>耕地>草地>其他用地。5个土地利用类型中，林地、灌木林地和草地在3个季节土壤水分差异较小，其他用地和耕地差异则较大，以冬季>春季>秋季。其中，耕地在秋季地表温度相对春季和冬季较高，而秋季又是农作物丰收时期，地表植被覆盖减少，导致蒸散发增大，所以表层土壤水分最低；冬季则是一年中温度最低季节，地表蒸散发较少，所以土壤水分含量可能变高；春季是耕作季节，表层土壤被翻转，加之温度上升，地表蒸散发相对较大，所以土壤水分比冬季低。

图8 土地利用类型

图9 不同季节下土地类型土壤水分含量

4 结论

利用Landsat8 OLI_TIRS数据，结合温度植被干旱指数TVDI，对贵州省喀斯特涟江流域土壤水分反演评价，得出以下结论。

a) 研究区内3个季节旱情分布较广，以轻旱旱情为主，面积所占比例达30%以上，春季>秋季>冬季。干旱次之，几乎无重旱。该区旱情主要发生在流域北部、中部和南部地势较为平坦的耕作区域和城镇周围，TVDI适用于流域的土壤水分含量反演及干早监测评价。

b) 冬季和春季土壤水分均值随着坡度增大而增加，秋季则是先随着坡度增加，再到25°以上陡坡区域减少。相同坡度条件下，冬季土壤水分总体比秋季和春季多，春季与秋季差值相对较小，冬季则较大。

c) 土壤水分均值在不同坡向有较大差异，由无坡向、北坡按顺时针旋转至西北坡，土壤水分先增加再减少，然后又先减少再曾加。总体上同一坡向条件下土壤水分含量阴坡比阳坡和无坡向多，冬季比秋季和春季多，春季与秋季的差值相对较小，冬季则较大。

d) 在耕地、林地、灌木林地、草地、其他用地五个地类中，林地在不同季节土壤水分含量最高，灌木林地次之，其他3个地类则存在差异。相同地类在不同季节、林地、灌木林地和草地的土壤水分差异较小，其他地类和耕地的差异则较大，冬季>春季>秋季。

影响植被指数与地表温度空间特征关系的因素还有很多，而在贵州喀斯特山地地区，海拔落差、云雾、空气水分含量、降雨及研究区域大小和地表覆盖物类型都会对植被指数与地表温度空间特征关系产生影响。若在该经验模型中不能综合考虑这些因素，会引起一定偏差。