黄土高寒区地处干旱半干旱气候区, 降水量少且分配不均, 生态环境十分脆弱。退耕还林还草工程的实施, 使该区域的水土流失得到了极大改善, 但不考虑土壤水分条件, 盲目进行人工植被建设可能会导致土壤干层, 生态环境恶化^[1]。土壤水分作为该区域限制植被生长的关键因子、流域水量平衡和区域水文循环的重要因素^[2-3], 具有重大的研究价值。近年来, 有诸多关于该区域林地土壤水分特征、土壤水分分布及林地土壤水分利用等方面的研究, 例如：汪星等^[4]研究了山地密植枣林0—1000cm土层范围的土壤水分特征变化规律；贾小旭等^[5]研究了黄土高原区域尺度土壤水分的时空分异特征及其影响因素；常恩浩等^[6]研究了黄土区5种典型退耕植被群落的根系分布及其水分利用特征；董爽等^[7]研究了黄土区流域分水线深层土壤水分时空动态及干层分布特征；李蓝君等^[8]研究了黄土沟壑区小流域刺槐、侧柏、油松三种典型造林树种蒸散发对气候变化的响应。但这些研究大多采用经典统计方法来描述相关因子的变化, 并不能反映这些因子在整个研究尺度空间上如何具体变化, 而这恰好是地统计学方法的优势之处, 并且地统计学已经被证明是分析土壤特性空间分布特征及其变异规律的最为有效的方法之一^[9]。

本研究选择位于黄土高寒区的大通县安门滩小流域典型人工林地作为研究对象, 分别对5月、7月和9月的林地土壤水分和耗水量进行定位观测, 并采用地统计学方法, 研究小流域内典型人工林分的土壤水分、林地耗水及土壤水分盈亏量的时空分布特征, 旨在分析验证该区域典型林分土壤水分状况是否与林地耗水分布格局相匹配, 能否维持水分收支平衡, 从而为该区域的生态建设结构优化提供理论依据。

1 研究区概况

大通县位于青海省东部(东经100°51'—101°56'、北纬36°43'—37°23'), 地处祁连山地与黄土高原的过渡地带, 海拔2280—4622m。全县年平均温度为2.8℃, 无霜期97d, 年降水量450—820mm不等, 多年平均降雨量为508 mm, 属于大陆性高原气候。降水的地区差异较大, 从东南到西北降雨量逐渐增加。降水的季节分布也很不均匀, 多集中在7、8、9三个月内。生长季内(5—9月)降雨为400—600mm, 约占全年降水量的87%左右, 雨热基本同季。全县年平均蒸发量为1290mm, 湿润指数的变化范围为0.56—1.32。土壤的主要成土母质为黄土和黄土状物质, 另有少量坡积母质和残积母质等。

研究区内共有乔灌木树种119种、隶属25科48属, 其中裸子植物门有3科6属10种。乔木树种主要有青海云杉(Picea crassifolia)、祁连圆柏(Juniperus przewalskii Kom.)、青杨(Populus cathayana)、华北落叶松(Larix principis-rupprechtii)等。其中安门滩小流域的主要林分有18a生的青海云杉、祁连圆柏、油松(Pinus tabuliformis Carrière)、柠条(Caragana intermedia)、白桦(Betula platyphylla), 20a生的华北落叶松和30a生的青杨等。林下草本主要有早熟禾(Poa annua L.)、鼠掌老鹳草(Geranium sibiricum L.)、白莲蒿(Artemisia sacrorum Ledeb)、披碱草(Elymus dahuricus Turcz.)、刺儿菜(Cirsium setosum (Willd.) MB.)、风毛菊(Saussurea japonica (Thunb.) DC.)、乳白香青(Anaphalis lactea Maxim.)、飞廉(Carduus nutans Linn)和狗娃花(Heteropappus hispidus (Thunb.) Less.)等。

2 研究方法 2.1 观测方法 2.1.1 样点布设

通过无人机航拍, 获取精度为5m的DEM影像图(图 1), 结合野外实地勘测结果, 绘制植被类型图。并采用随机布设的方法共布设样点87个, 样点和植被的具体分布情况见图 2。

2.1.2 土壤水分测定与储水量、盈亏量计算

2017年8月, 用与CPN-503探头式中子土壤水分仪相配套的钻机在各样点较平坦处, 垂直地面打一直径略小50毫米、深度略大于100cm的钻孔, 在钻孔中插入外径为50毫米、内径为46毫米底部密封的1.1m长PVC管, 使其露出地面5—7cm, 并用塑料布将上口封住。采用中子土壤水分仪测定土壤体积含水量, 测定深度为100cm, 每20cm测定一次。监测时间为分别2018年5月、7月和9月的月初、月中和月末。测定前需在各样点对中子仪进行标定, 标定时, 0—20cm表层土壤, 用TDR测量中子管附近5cm范围内的表层土壤含水量, 20—100cm土层, 用土钻取样后采用烘干法测定土壤体积含水量, 从而获得标定方程。

土壤储水量按公式(1)计算。

(1)

式中, W_i为第i层土壤储水量(mm)；θ_v为第i层土壤体积含水量；h_i为第i层土壤的厚度, 此处为20cm。各月份的土壤储水量取每月三次测定的平均值。

阶段土壤水分盈亏量按公式(2)计算。

(2)

式中, ΔW为阶段土壤水分盈亏量(mm), 若取值为正, 说明该阶段内土壤水分有所盈余, 若取值为负, 说明该阶段内土壤水分有所亏损；W_a、W_b分别为阶段初、末1m土层内的土壤储水量(mm)。

2.1.3 林地耗水量测定与计算

在研究区典型区域设置若干Davis Vantage Pro2自动气象站测定同期降雨量与气温, 并在各样点林内沿等高线布设100cm×20cm×25cm雨量槽测定穿透雨、林外四周布设雨量筒测定林外降雨, 二者之差即林冠截留量。研究区生长季内降雨量与气温特征见图 3。

2.1.4 林地耗水量测定与计算

在各样点内采用瞬时剖面法^[10], 使用日产DIK 3130式张力计, 每隔10cm分层定时观测土壤土水势, 结合体积含水量, 确定相应范围内的水分特征关系曲线与非饱和土壤导水率^[11]。应用定位通量法^[12], 可计算出各样点地下部分的耗水量, 其数值与实测林冠截留之和, 为林地耗水量。在本研究点完成的相关林分调查表明, 试验林分的根系分布范围主要集中在0—1m的土层中, 因此可将定位点Z_r设定在根际层以下处, 本文中Z_r取90cm, ΔZ=20cm。

2.2 分析方法

地统计学也被称为地质统计学, 是以结构分析与变异函数理论为基础, 是在有限的区域内对区域化变量进行无偏最优估计的一种方法^[12]。半变异函数也叫半方差函数, 是地统计分析中的特有函数, 区域化变量Z(x)在点x和x+h处的值Z(x)与Z(x+h)差的方差的一半叫做区域化变量Z(x)的半变异函数^[13-14], 记为γ(h)。根据定义有：

(3)

式中：γ(h)为h的半方差函数值；h为两样本点的空间距离(m)；N(h)为间隔距离等于h的样本点的对数目；Z(x_i)为空间位置点x_i处样本的实测值；Z(x_i+h)为空间位置点x_i+h处样本的实测值。

本研究使用IBM Statistics SPSS 20.0统计软件分别对各样点的土壤储水量、林地耗水量和土壤水分盈亏量

进行描述性统计分析, 用K-S检验法(Kolmogorov-Smirnov test)进行正态分布检验, 显著性水平设定为α=0.05。利用GS+9.0软件和ESRI ArcGIS10.2软件进行半方差函数拟合、Kriging插值和空间分布图绘制等地统计学分析。根据已有等高线矢量数据, 使用ESRI ArcGIS10.2软件绘制出研究区的坡向图和坡度图, 详见图 4。

3 结果与分析 3.1 土壤储水量、林地耗水量和土壤水分盈亏量的统计特征和变异结构分析

土壤储水量、林地耗水量和土壤水分盈亏量的统计特征值列于表 1。由土壤储水量和林地耗水量的取值范围与均值可以看出, 土壤水分状况总体表现为9月最优, 其次为5月, 7月最差, 9月的土壤储水量较7月增加了27%左右, 这与赵荣玮^[15]、孔凌霄等^[16]的研究结果相同, 与图 3所呈现出的研究区生长季内降雨量与气温的变化情况是相呼应的。而林地耗水量则呈现出7月最大, 其次是9月, 5月最小的特点, 7月的林地耗水量为5月的3倍左右。研究期间, 植被类型和组成并未发现明显变化, 因而林地耗水变化主要是由于5月气温较低, 林木草本刚开始萌动, 7月气温上升, 植被生长繁盛, 土壤蒸发剧烈且植被耗水增加, 9月气温降低, 植被长势减弱, 土壤蒸发与植被耗水都减弱。由土壤水分盈亏量的取值范围和均值可以看出, 5—7月绝大多数林地的土壤水分呈亏损状态, 而7—9月所有林地的土壤水分都获得了一定程度的补充, 总体来看, 5—9月研究区多数林地的土壤水分有所盈余, 部分林地的土壤水分有所亏缺, 这与张进虎^[17]的研究结果一致。从变异系数来看, 土壤水分盈亏量变异最大, 土壤储水量次之, 两者属于中等变异(10%≤CV%≤100%), 林地耗水量变异最小。表 1中显示, 土壤储水量、林地耗水量和土壤水分盈亏量经K-S检验的P值均大于0.05, 即均服从正态分布, 因此都不必进行数据转换, 可直接进行地统计学分析。

经比对不同半变异函数模型产生的参数, 土壤储水量、林地耗水量和土壤水分盈亏量的最优理论变异函数模型均采用指数模型(表 1)。Co表示块金值, Co+C为基台值, Co/(Co+C)可反映区域化变量的空间自相关性大小。对于土壤储水量、林地耗水量和土壤水分盈亏量, 25%≤Co/(Co+C)≤75%, 表明三者均具有中等强度的空间自相关性。表 2为土壤储水量、林地耗水量和土壤水分盈亏量的交叉验证结果, 插值精度主要由预测平均误差ME、标准平均误差MSE和标准均方根误差RMSSE三个指标作为评价标准, ME和MSE越接近于0, RMSSE越接近于1表明插值精度越高。由表 2可知, 三者预测结果的ME和MSE均接近于0, RMSSE均接近于1, 插值结果的精度可以接受。此外, ME和MSE若为正值说明插值结果较实际值偏高, 若为负值说明插值结果较实际值偏低。可见, 各月土壤储水量和5—9月期间土壤水分盈亏量的插值结果较实际值偏低, 而各月林地耗水量和5—7月、7—9月期间土壤水分盈亏量的插值结果较实际值偏高。

3.2 土壤储水量、林地耗水量和土壤水分盈亏量的空间插值分析 3.2.1 土壤储水量的空间分布格局

林地土壤储水量是指在自然状况下, 林地某一特定土壤层次土壤的含水量, 包括重力水、毛管水和吸着水。土壤中水分含量状况极大地影响土壤物理性质, 并对整个生态系统的水热问题起决定作用^[18], 而流域内的气象、植被和地形等因素均会显著影响林地土壤储水量。由于研究区面积较小, 小流域内林地的降雨、温度等气象因子基本一致, 植被和地形因素起主导作用。研究区5月、7月和9月的土壤储水量分布情况如图 5所示。可以看出, 土壤储水量的空间分布呈条带状和斑块状格局, 且不同月份土壤储水量的空间分布格局存在较强的相似性。各月土壤储水量总体上呈南高北低、西高东低的空间分异规律, 且在西南-东北方向变异较东南-西北方向剧烈。结合研究区等高线图可知, 试验地土壤储水量的空间分异规律与其海拔高度相一致, 海拔越高, 土壤储水量越高, 这与王甜^[19]、段良霞等^[20]的研究结果相同。这可能是由于研究区为阶级状梯形分布, 坡度较缓, 天然降水很难顺坡流下, 且相对下坡位, 上坡位的林分较为矮小、稀疏, 土壤水分消耗较少。土壤储水量在西南-东北方向变异较东南-西北方向剧烈, 这很可能是东南-西北方向的地形为“山脊-鞍部-山脊”, 呈“U”型, 而西南-东北方向为“山脊-山脚”, 海拔变化区间较大所导致的。土壤储水量总体表现为西北部区域低于东南部区域, 结合图 4可知, 西北部为阳坡, 东南部为阴坡, 即阳坡的土壤储水量总体低于阴坡, 这与孔凌霄^[16]、史志华等^[21]的研究结果一致。

此外, 5月、7月和9月土壤储水量均在西南部和中部区域出现了高含量区, 土壤储水量分别达到224.82—253.39mm、208.94—240.35mm和248.54—307.08mm；东北部边缘区域出现了低含量区, 土壤储水量分别达到164.55—184.75mm、135.36—169.26mm和166.25—212.09mm。结合植被分布情况可知, 研究区西南部区域栽植的主要为山杏、白桦和青海云杉, 中部区域主要为华北落叶松, 东北部边缘区域主要为青杨。5月白桦、山杏林地土壤储水量最高, 其次是青海云杉、油松和华北落叶松林地, 斑叶稠李、祁连圆柏、柠条林地土壤储水量较低；7月山杏、青海云杉林地土壤储水量最高, 其次是白桦和油松林地, 斑叶稠李、柠条、祁连圆柏和华北落叶松林地土壤储水量较低；9月白桦林地土壤储水量最高, 其次是青海云杉、山杏和油松林地。各月均呈现出青杨林地土壤储水量最低的特点。

3.2.2 林地耗水量的空间分布格局

由图 6可以看出, 5月与9月林地耗水量的空间分布较为相似, 均呈明显的条带状格局, 而7月林地耗水量分布呈明显的斑块状格局, 这是7月林地耗水量的空间变异性较强所导致的。与土壤储水量的分布格局相同, 各月林地耗水量在西南-东北方向变异较东南-西北方向剧烈, 总体表现为西南部区域低于东北部区域。这是因为西南部区域所栽植的山杏和祁连圆柏生长较为矮小, 加之树种自身性质, 耗水量相对较小；而东北部边缘区域种植的青杨和中部区域种植的华北落叶松, 均生长高大茂盛, 耗水量相对较大。

5月、7月和9月林地耗水量均在东北部边缘和中部区域出现了高耗水区, 耗水量分别达到38.91—45.58mm、107.99—122.25mm和43.43—53.08mm；在西南部边缘区域均出现了低耗水区, 耗水量分别达到29.59—33.49mm、86.61—96.53mm和30.80—36.24mm。结合植被分布情况可知, 5月林地耗水量最高的是青杨和华北落叶松, 其次是斑叶稠李、青海云杉、白桦和油松, 祁连圆柏、柠条和山杏最低；7月林地耗水量最高的是青杨和斑叶稠李, 其次是白桦、华北落叶松、青海云杉、柠条和油松, 祁连圆柏和山杏最低；9月林地耗水量最高的是青杨和斑叶稠李, 其次是华北落叶松、柠条、青海云杉、白桦和油松, 山杏和祁连圆柏最低。其中, 作为阔叶树种的青杨和斑叶稠李, 林地耗水量始终较大, 而柠条由5月刚开始萌芽至长叶开花, 耗水量逐渐变大。

3.2.3 土壤水分盈亏量的空间分布格局

土壤水分盈亏量主要是降水量与同期蒸散量的差值, 它的取值及分布规律决定于降水量和影响植物蒸散量的太阳辐射、气温、湿度和风速等气象要素^[22]。研究区5—7月、7—9月和5—9月的土壤水分盈亏量分布情况如图 7所示, 不难看出, 这三个时期内的土壤水分盈亏量在空间上均呈条带状和斑块状分布格局, 三者间存在一定的相似性却又有所区别。相似性为：在5—7月、7—9月和5—9月这三个时期内, 土壤水分盈亏量的取值均呈现出东北部区域小于西南部区域的特点, 在降雨量基本相同的情况下, 这是东北部区域的林地耗水量大于西南部区域造成的。不同之处在于：5—7月研究区土壤水分盈亏量除东南部边缘小部分区域外, 总体取负值, 说明5—7月期间研究区土壤水分总体为亏损状态。土壤水分盈亏量表现出的趋势为南高北低, 且在东南部边缘区域出现高值, 在东北部边缘和中部区域出现低值。该期间华北落叶松和白桦林地的土壤水分亏损最为严重, 其次是青杨、山杏、柠条、祁连圆柏、油松和青海云杉林地, 斑叶稠李林地的土壤水分亏损最少。7—9月研究区土壤水分盈亏量都取正值, 说明7—9月期间研究区土壤水分总体为盈余状态。土壤水分盈亏量表现为西南部区域高于东北部区域, 且在西北部边缘和南部边缘区域出现高值, 在东北部边缘和中部区域出现低值。该期间土壤水分盈余量最多的是白桦和柠条林地, 其次是斑叶稠李、油松、祁连圆柏、青海云杉和山杏林地, 青杨和华北落叶松林地最少。5—9月研究区土壤水分盈亏量除东北部边缘和中部小部分区域外, 总体取正值, 说明5—9月期间研究区除东北部边缘和中部小部分区域林地的土壤水分有所亏缺外, 大部分林地的土壤水分有所补充。土壤水分盈亏量表现为西南部区域高于东北部区域, 且在西北部边缘和东南部边缘区域出现高值, 在东北部边缘和中部区域出现低值。在整个生长季期间, 土壤水分盈余量最多的是斑叶稠李和柠条林地, 其次是油松、青海云杉、祁连圆柏、白桦和山杏林地, 青杨和华北落叶松最少, 且东北部边缘处的青杨林地和中部区域的华北落叶松林地出现了土壤水分亏损现象。

在生长季期间, 斑叶稠李、柠条、油松、青海云杉、祁连圆柏、白桦和山杏等林地均能维持水分收支平衡, 只有东北部边缘处的青杨林地和中部区域的华北落叶松林地出现了土壤水分亏损现象。其中, 青杨林地在5月、7月和9月的耗水量均为最高, 天然降雨补充量小于林地蒸散量, 故而青杨林地在5—7月的土壤水分亏损量较多, 在7—9月土壤水分盈余量最少, 在整个生长季期间, 土壤水分有所亏损；而华北落叶松林地除5月外, 在7月和9月的耗水量均不突出, 但其在5—7月的土壤水分亏损最为严重, 在7—9月的土壤水分盈余量最少, 且在整个生长季期间, 土壤水分有所亏损, 这很可能是华北落叶松林地的部分土壤水分向深层土壤发生了渗漏而导致的。

4 结论

本研究运用地统计学方法对小流域典型林分土壤水分状况与林地耗水特征进行综合分析, 结果表明：

(1) 土壤水分状况总体表现为9月最优, 其次为5月, 7月最差, 而林地耗水量表现为7月最大, 其次是9月, 5月最小, 5—7月绝大多数林地的土壤水分呈亏损状态, 而7—9月所有林地的土壤水分都获得了一定程度的补充, 总体来看, 5—9月研究区多数林地的土壤水分有所盈余, 部分林地的土壤水分有所亏缺。土壤储水量、林地耗水量和土壤水分盈亏量均可采用指数模型作为最优理论变异函数模型；

(2) 5月、7月和9月土壤储水量总体上呈南高北低、西高东低的空间分异规律, 且西南-东北方向变异较东南-西北方向剧烈；各月林地耗水量在西南-东北方向变异较东南-西北方向剧烈, 总体表现为西南部区域低于东北部区域；在5—7月、7—9月和5—9月这三个时期内, 土壤水分盈亏量的取值均呈现出东北部区域小于西南部区域的特点。

可见, 当地土壤水分状况与林地耗水量分布格局并不完全匹配, 虽然绝大部分林分能够维持土壤水分收支平衡, 在生长季结束时有水分盈余, 但部分山脚处的青杨林地和中部区域的华北落叶松林地出现了土壤水分亏损的现象。为防止流域内水分环境恶化, 在之后黄土高寒区的植被建设过程中, 应适当调整林分配置, 在土壤水分条件不允许的情况下, 不宜继续栽植青杨等耗水量过大的树种, 从而保证林分正常生长, 实现水土保持措施可持续发展。