基于SVM的绿洲荒漠交错带土壤水分与地下水埋深反演

张钧泳 丁建丽 谭 娇

(1.新疆大学资源与环境科学学院， 乌鲁木齐 830046；2.新疆大学智慧城市与环境建模自治区普通高校重点实验室， 乌鲁木齐 830046；3.新疆财经大学计算机科学与工程学院， 乌鲁木齐 830032)

0 引言

地下水分布范围广泛，其可持续利用性强，在水资源中占有重要地位。浅层地下水作为地下水组分之一，是干旱植被生长环境的限制变量，地下水资源循环周期长，一旦赋存环境遭到破坏，再生能力减弱。随着人类社会发展，长期不合理、无规则开采使地下水水位持续下降，地下水系统遭受破坏[1]。大多数干旱内陆河流域面临严峻的大型灌溉农业用水[2-3]，地下水是水源补给的主要方式，故干旱区绿洲生态环境的调控作用显得尤为重要[4]。合理的浅层地下水位能满足干旱区绿洲植被的正常生长需求，这使得浅层地下水、植被与土壤之间的相互作用和植被群落分布特征与浅层地下水的相互响应研究尤为重要[5]。遥感技术具有宏观动态快速监测的特点，避免了传统检测方法中时间长和精力消耗多的缺点，在监测地下水位应用中可以快速获得地下水位，实现大面积监测[6]。随着卫星遥感技术不断进步，定量反演地表特征参数，结合相关实测数据，建立地下水相关模型，可快速获得地下水相关信息[7-11]。

目前，遥感手段在大尺度区域反演地下水埋深研究中已经成为主流趋势，但同时地下水遥感反演也是一个技术难题，地下水遥感监测主要是通过地下水与地表水[12]、地表温度[13]、土壤水分[14]和植被[15]等相关性进行研究。目前，对地下水进行遥感探测有以下4种类型：①热红外遥感监测法，早期该方法利用热红外波段判断地下水的存在，现今姜红等[12]分析了新疆开都河流域周围的绿洲，监测和分析该区域17年间的地下水位变化；付俊娥等[13]利用以MODIS遥感影像计算植被指数(NDVI)以及地表温度(LST)数据，通过SVM定量估算了河西走廊疏勒河流域的地下水位；尹涛等[14]通过1 d中不同时间段DNVI、LST进行数学变换，反演了在植被生长期的地下水埋深情况。②环境因素遥感信息法，根据与地下水相关的环境因子与地下水间的关系，来监测和反演区域地下水状况。田凯等[15]利用MODIS数据并计算NDVI以及条件植被覆盖率(VCR)，与野外实测的地下水埋深相结合，发现黄河中上游地区的植被生长所适合的地下水埋深。③遥感信息定量反演分析法，该方法案例较多，如毛德发等[16]在北京市大兴区消耗水用量以及本着地下水可持续利用的条件下，利用MODIS数据以及气象降水径流等数据，得出该区域综合节水的有效方法。④微波遥感监测分析法，RODELL等[17]利用GRACE卫星数据，对密西西比河主要的子流域地下水储量变化进行了估算。KOMAROV等[18]在水文因素、土壤介电常数和土壤对雷达波段的反射特性分析的基础上,建立土壤水分和地下水之间的实验方程。尹楠等[19]在不同极化模式下，研究了后向散射系数对垄行方向变化的敏感性及去除垄和地表粗糙度对土壤含水率反演的影响。综上所述，利用微波遥感、土壤水分、植被指数和热红外等基于遥感信息定量监测地下水位已经取得了较大的进展，其次，干旱植被指数在土壤水分和地下水的监测中有着至关重要的作用[20-21]。

本文以新疆渭干河-库车河绿洲(以下简称：渭-库绿洲)为研究区，选择Sentinel-1A影像和Landsat-OLI影像为数据源，通过微波和光学遥感数据以及土壤含水率、地下水埋深数据，结合植被以及土壤条件，利用支持向量机算法进行定量反演研究区土壤含水率以及地下水埋深，为监测绿洲的地下水埋深提供技术支持。

1 研究区概况和数据资料

1.1 研究区概况

渭-库绿洲位于亚欧大陆腹地。本文所界定的野外实地调查采样区间为81°26′～83°17′E，41°06′～41°56′N。从地貌条件看，北部地区海拔最高，其次为内部平原区，南部沙漠地区海拔最低，是一个典型的冲积扇平原。气候干燥少雨，气温呈北低南高的特点。绿洲内部的光照和热量资源十分充足，土壤类型丰富，主要农业经济作物为小麦、玉米、棉花、石榴和无花果等。植被类型多为盐生植被，呈片状分布，多分布于绿洲外部、绿洲东北部荒漠带及绿洲荒漠交错带[22]。土壤水分反演区域为渭-库绿洲全区，如图1a所示，地下水埋深反演区域为库车绿洲的东北角，如图1b所示。

图1 研究区位置图Fig.1 Location of study area

1.2 实测数据

2015年3月21日—4月30日对研究区开展了一次野外调查，调查的内容主要是了解春灌之前稀疏植被情况以及采集土壤样本，获取土壤含水率、温度、质地、容重、介电常数以及地下水埋深数据。根据实验室多年数据的累积，选择了38个具有代表性的采样点，均匀分布在绿洲、绿洲荒漠交错带以及荒漠带区域。其采用五点法进行土壤样本的采集，使用W.E.T型传感器、Hydra型土壤测试仪测量土壤含水率、温度以及介电常数等数据，土壤电导率(EC)和总溶解固体(TDS)含量由inolab cond 7310型台式电导率测试仪(德国WTW公司)测定，pH值由pH7310型台式酸度计(德国WTW公司)进行测定，土壤容重使用环刀法进行测定。使用美国Onset公司生产的HOBO型自动记录水位计(以下简称水位计)记录地下水埋深数据，该仪器可根据实验的要求设定水位计的记录周期，水位计产生的误差最大为1 cm，降低了实验误差。

1.3 遥感数据

Sentinel-1A是合成孔径雷达(Synthetic aperture radar，SAR)的一种，其主要的特点是空间分辨率高，工作状态不受云雨天气影响和波动，重访周期短，可全天时、全天候监测，具有单双极化方式，较强的干涉能力等优点[23]，在全球范围内可用于森林监测、农业监测、地震监测、陆地冰雪提取、洪水监测以及土壤水分、土壤湿度监测等[24-29]。本文Sentinel-1A数据来源于欧空局的哨兵科学数据中心(Sentinels Scientific Data Hub)，影像成像时间为2015年4月21日，为地距影像(Ground range detected，GRD)Level 1级产品，极化方式为VV，工作方式为干涉宽幅模式(Interferometric wide，IW)，分辨率为5 m×20 m，与太阳同步轨道，载波波段为C波段，工作频率为5.044 GHz，轨道高度为693 m，重访周期为12 d，入射角为20°～40°，幅宽240 km。

使用遥感数据为对地观测卫星Landsat OLI的影像，获取时间为2015年4月26日，空间分辨率为30 m，数据来源于http:∥glovis.usgs.gov/。

1.4 遥感影像预处理

利用欧空局(ESA)提供的Sentinel-1 Toolbox(S1TBX)软件对微波图像进行预处理：①辐射校正。利用S1TBX软件提供的辐射校正工具完成辐射校正，得到影像后向散射系数。②噪声处理。首先对图像进行热噪声移除；其次，为了降低影像中的斑点噪声，采用Refined Lee滤波进行处理，滤波器窗口设置为7像素×7像素。③几何校正。先将处理好的雷达图像转换成后向散射系数，然后进行几何校正。将处理之后的Sentinel-1A SAR数据输出为ENVI数据格式，并将其重采样到与Landsat 8数据具有相同的空间分辨率。

Landsat OLI需要进行辐射校正、大气校正、几何校正以及图像裁剪等预处理操作，具体操作参照文献[30]。

2 研究方法

2.1 不同温度-植被指数构建

SANDHOLT等[31]发现不同植被指数与温度指数均呈不规则简化三角形分布，在此基础上提出了温度植被干旱指数的概念。Ts-VI由植被指数和地表温度计算得到，其定义为

(1)

其中，对最大、最小地表温度散点进行线性拟合，得到TVDI特征空间的干湿边方程

Tsmax=a1+b1VI

(2)

Tsmin=a2+b2VI

(3)

式中TVDI——温度植被指数

Ts——地表任意像元的地表温度

VI——植被指数

Tsmin——相同植被指数值对应的最小地表温度

Tsmax——相同植被指数值对应的最大地表温度

a1、b1——干边线性拟合方程的系数

a2、b2——湿边线性拟合方程的系数

其中，Tsmin对应Ts-VI特征空间的湿边；Tsmax对应Ts-VI特征空间的干边。Ts-VI的阈值在0～1之间，当TVDI值逐渐向1靠近时，土壤干旱情况越来越严重；反之，当TVDI值越靠近0时，土壤湿度越高。因此Ts-VI与土壤湿度的相关性，在两种极端的情况下可以反映干、湿情况。

为消除大气散射的影响以及地表相邻点反射光折射的影响，利用ENVI图像处理软件通过数字高程(DEM)数据提取坡度、坡向数据。其操作步骤参照文献[32]，进而对LST数据进行地形校正。

2.2 土壤后向散射系数计算

为了得到比较准确的土壤表面后向散射系数，对有植被影响的后向散射系数进行处理，以消除植被所引入噪声带来的负面影响，从而获取区域表层土壤真实后向散射系数。根据研究区当时实际情况，选择适合本区域的Water-Cloud算法

(4)

(5)

γ2(θ)=exp(-2Bmvegsecθ)

(6)

γ2(θ)——农作物层的双层衰减因子

mveg——植被含水量，kg/m3

θ——雷达波入射角

其中，A和B分别为依赖于植被类型的参数，可以通过多次迭代衰减或最小二乘法获得。

估算土壤水分时应规避植被覆盖影响，则需要应用Water-Cloud模型消除植被层在土壤水分后向散射中的干扰。研究学者对于A、B参数的合理取值及地域性分析进行了研究，本文采用周鹏等[33]所标定的符合渭-库绿洲的最优解，即参数A=0.001 9，B=0.137。

在利用Water-Cloud模型进行有植被区的土壤水分定量估算时，其中一个至关重要的输入参数为植被含水量(VWC)。研究区植被多为棉花及低矮的植被，为了获取研究区的植被含水量，利用改进型归一化植被水分指数NDMI[34]，并根据经验方程[35]，则VWC计算式为

VWC=2.15NDMI+0.32

(7)

其中

NDMI=(ρNIR-ρMIR)/(ρNIR+ρMIR)

(8)

式中ρNIR——遥感影像近红外波段反射率

ρMIR——遥感影像中红外波段反射率

将式(7)代入式(5)与式(6)，运用Water-Cloud模型消除植被覆盖影响从而得到裸土后向散射系数，表示为

(9)

γ2(θ)=exp(-2×0.137mvegsecθ)

(10)

(11)

通过式(5)和式(11)可以计算出研究区有无规避植被影响的土壤后向散射系数，结果如图2所示。

图2 植被效应与去除植被效应后向散射系数比较Fig.2 Comparison of backscatter coefficients before and after removal of surface vegetation cover in April

2.3 SVM算法构建

已有研究证明，与其他核函数相比，径向基核函数(RBF)在土壤含水率预测模型中效果更优[36]。本文选择RBF核函数建立土壤含水率预测模型，其功能选择为回归功能。目前研究发现微波遥感数据对土壤水分变化的敏感度也非常强，而且微波数据具有实时观测能力，但仅考虑后向散射系数和土壤介电特性等因素的微波遥感反演土壤水分是远远不够的[37]，因其未考虑地表上的植被覆被等因素。因此，本文选择国内外水盐遥感研究最常用的多光谱数据(Landsat OLI)和合成孔径雷达数据(Sentinel-1A)，选择支持向量机(Support vector machine，SVM) 算法，此方法在土壤水分监测方面反演精度较高[38]。将干旱指数、土壤含水率、后向散射系数进行耦合，区域地下水埋深与影响因子之间的相互关系可以看作复杂的非线性函数问题。支持向量机非线性基本思想是：选择适合的核函数，通过空间变换(将输入向量z映射到一个高维的特征空间上)，实现其线性可分性，使得到的结果风险小。根据泛函数的相关理论，只要核函数K(xi，yi)满足Mercer条件，它就对应于某一变换空间中的内积，公式为

(11)

而对应的最优分类函数形式可写为

(12)

式中n——训练样本个数

xi、yi、yj——训练样本向量

b*——偏置项系数

2.4 支持向量机回归算法操作流程

基于Matlab平台，以土壤含水率预测为因变量，使用支持向量机回归算法进行模型的构建，其基本思想步骤如下：①建模集样本和验证集样本的划分，本文选择26个样点用于模型的训练，12个样本用于模型的验证。②SVM回归模型的创建，利用svmtrain函数进行样本训练，在训练时考虑数据的实际情况对数据进行归一化处理。选择高斯径向基核函数(RBF)，并通过交叉验证网格搜索法(Grid-search)，利用训练样本反复地调整惩罚因子c和参数g，获得最佳的参数组合。③利用svmpredict函数实现回归模型的预测。④精度评价通过svmpredict函数进行预测时产生的均方根误差RMSE和决定系数R2等评价标准，即可对建立的SVM回归模型的预测能力进行综合评价。

3 结果与分析

3.1 地下水埋深和土壤含水率的相关性分析

野外实际考察发现，农业生产中，当地于1月初进行冬灌和4月底进行春灌，以改善土地的水分条件。2015年1— 4月期间，没有灌溉用水、降雨量极少(图3)，且绿洲蒸发量极大，此时土壤中的水分来自地下水埋深的补给。

图3 月平均降雨量Fig.3 Monthly average precipitation

通过分析不同深度的土壤含水率与地下水埋深之间的相互关系来探求两者之间的线性关系，通过相关性分析得出(表1)，0～10 cm的土壤含水率与地下水埋深的决定系数达到了0.792 2，拟合效果相对较好，其次是40～60 cm的土壤含水率，其决定系数为0.722 3，10～20 cm和20～40 cm的土壤含水率与地下水埋深的决定系数分别为0.311 8和0.103 1，拟合效果相对较差。地下水埋深和不同深度的土壤含水率的拟合效果由高到低依次为0～10 cm、40～60 cm、10～20 cm、20～40 cm。

表1 地下水埋深与不同深度的土壤含水率的回归分析Tab.1 Regression analysis of groundwater depth and soil moisture content at different depths

注：** 表示在0.01水平上极显著，*表示在0.05水平上显著。

结果表明，0～10 cm的土壤含水率与地下水埋深的拟合效果最优，土壤表层的土壤含水率与地下水埋深之间的相互关系最为密切，利用土壤表层的含水率来反演地下水埋深是可行的。

3.2 土壤水分反演模型精度比较

从图4a～4d可看出，当只考虑TVDI作为SVM模型参数时，TVDIMSAVI、TVDINDVI、TVDIEVI、TVDIDVI均可以反映土壤表层的土壤湿度状况，可用于土壤表层的水分反演。建模集R2分别为0.74、0.66、0.71和0.55，RMSE分别为4.66%、4.96%、4.71%和4.83%，TVDIMSAVI在SVM建模的结果R2为0.74，RMSE为4.66%，验证集中的R2分别是0.70、0.69、0.74和0.59；RMSE分别为4.65%、4.91%、4.83%和4.56%。验证结果与建模结果基本相同，在4种不同植被指数构建特征空间得到的干旱指数中，TVDIMSAVI精度最优，其次为TVDIEVI，最差的为TVDIDVI。建模集和验证集不同植被指数的TVDI预测出来的土壤含水率相关性由大到小排序为：TVDIMSAVI、TVDIEVI、TVDINDVI、TVDIDVI。

表2 不同参数的SVM模型预测方案Tab.2 SVM model prediction scheme with different parameters

图4 基于不同参数的SVM建模精度和验证精度Fig.4 SVM modeling accuracies and verification accuracies with different parameters

表3 SVM回归模型的最佳参数组合Tab.3 SVM regression model of the best combination of parameters

3.3 土壤水分空间分布特征

图5 基于TVDI和参数反演的土壤含水率空间分布Fig.5 Spatial distribution of soil moisture content based on TVDI and parameter inversion

3.4 地下水埋深反演

图6 土壤含水率反演的地下水埋深Fig.6 Retrieved result of groundwater depth by soil moisture content

3.5 反演结果验证

在得到研究区地下水埋深分布图后，选取HOBO水位计的数据与反演结果进行验证，对水位计实测数据和反演结果进行比较，如图7所示。SVM模型反演的地下水埋深与实测值的平均误差为8.23%，优于李相等[40]得出的18.06%的结果。

图7 预测值与实测值结果比较Fig.7 Comparison of predicted and measured results

4 结论

(1)通过对不同深度土层的土壤含水率与地下水埋深进行相关性分析，得出与地下水埋深相关性最优的为0～10 cm的土壤含水率，40～60 cm土层的土壤含水率次之，10～20 cm和20～40 cm的土层相关性较低。