天然针阔混交林下穿透降雨特征及其模拟

周 威，康永祥，刘建军，夏国威，刘小林

（1.西北农林科技大学 林学院，陕西 杨凌712100；2.贺兰山森林生态系统定位研究站，宁夏 银川750021；3.甘肃小陇山森林生态系统定位研究站，甘肃 天水741022）

水在森林生态系统中的循环与分配整合了能量流动和养分循环等生态过程［1］，因而水文功能在森林生态系统的作用中是非常重要的。张焜［2］、王波等［3］、孔维健等［4］、刘小林等［5］分别在不同地域对天然针阔混交林的水文功能进行了研究，为探求天然针阔混交林的水文作用提供了有价值的参考。但由于研究对象的林分特征、气候条件、研究侧重点等存在差异，要全面认识天然针阔混交林的水文效应还任重而道远。

森林引起生态系统内降雨分配的时空差异，对林地水文过程产生显著影响［6］。国内外针对林冠层分配降水研究较多，但在秦岭西段小陇山的研究鲜见报道［7］。小陇山林区是兼有中国南北气候特点的典型天然次生林区，天然次生针阔混交林是一个重要森林类型。研究小陇山林区森林冠层下穿透降雨规律及与大气降雨之间的相关性，对揭示区域森林生态水文作用具有重要意义。

在森林生态系统中，不同冠层对降雨的截留效应响应不同，导致了林内穿透降雨的空间分布异质性［8］。地统计学是一种空间分析方法，用于描述区域化变量的空间分布特征，并提供了一种无偏最优的空间插值方法，可用来充分认识空间结构的特征，已广泛应用于矿产、土壤、生态、气象等因子分布图的绘制［9］。本文试图将地统计学的理论方法应用到穿透降雨的研究中，探讨穿透降雨空间分布的特征。

本文通过对小陇山天然次生针阔混交林下穿透降雨特征及其空间分布格局与降雨量关系的分析研究，旨在进一步明确针阔混交天然林生态系统的水分运转过程及功能的生态学机制，以求更深入地探究森林水文的过程及作用机理，以期为正确认识小陇山林区的森林水文效应及流域森林生态水文功能提供一定的观测数据和经验参数，对制定森林管理措施提供科学依据。

1 研究区概况

研究区位于甘肃省小陇山林业科学研究所沙坝实验基地，地处小陇山林区中心地带，行政区划属天水市秦州区娘娘坝镇，东接观音林场，南接高桥林场，西邻麻沿林场，北连李子园林场。属秦岭南坡石质山地，海拔在1 550～2 100m，年平均气温1.2～8.8℃，最高气温30.3℃，最低气温－22.4℃，≥10℃的年积温2 480.4℃，初霜期10月16日，终霜期5月4日，无霜期154d，年均降雨量757mm，平均年蒸发量1 012.2mm，平均相对湿度78%，属大陆性季风气候，为暖温带湿润区。土壤以山地褐土和山地棕壤为主，成土母质为岩石，土层厚度30—55cm，土壤质地为砂壤，表土层腐殖质含量丰富，矿物质养分含量中等，呈微酸性到中性。

样地设在沙坝实验基地乱石窖沟，地理坐标为105°53′72″—105°53′73″E，34°8′56″—34°8′57″N，海拔1 630～1 660m，样地坡度15°～25°，南坡中坡位。样地为以锐齿栎（Quercus aliena var.acuteserrata）为主的天然次生针阔混交林，几乎无人为干扰，植被繁密。样地内胸径大于5cm的乔木总计120株，平均胸径16.8cm，平均树高19.9m，平均郁闭度0.93，乔木林群落的树种组成为4锐齿栎，2华山松（Pinus armandii），1四照花（Cronus japonica var.chinen－sis），1灯台树（Bothrocaryum controversum），2软阔叶树种。样地内下木主要有白檀（Symplocos paniculata）、千金榆（Carpinus cordata）、鹅耳枥（Carpinus turczaninowii）、葛枣猕猴桃（Actinidia polygama）、灰栒子（Cotoneaster acutifolius）等，平均高度2.5m，盖度为60%。草本主要有菝葜（SmilaxL.）、羊胡子草（EriophorumL.）、悬钩子（Rubus L.）等，盖度15%。

2 研究方法

2.1 试验布设与数据收集

林外大气降雨量资料来源于小陇山森林生态系统定位研究站，位于甘肃省小陇山林业实验局林业科学研究所沙坝实验基地，地理坐标为105°54′E，34°07′N。在沙坝实验基地天然次生林中选取一块投影28.5m×27m的样地，借助测绳和罗盘仪沿等高线精确圈出样地。在样地内网格状机械布设100个雨量桶（直径20cm），保持雨量筒高度距地面60cm，按品字形水平布设，行与行、列与列之间投影距离均为3m，用于收集样地的穿透降雨量。每次降雨后及时测定各个雨量桶内的水量，计算雨量时根据收集面积将雨量桶内每次降雨量换算为mm单位。

2.2 数据分析

采用Microsoft Excel 2007和SPSS 18.0对数据进行一般处理分析，利用穿透降雨的测定数据，计算得到每次降雨事件每个雨量筒的穿透率（相同时间下林内穿透降雨量与林外大气降雨量的比值）。完成基本统计分析后，用SPSS 18.0对数据通过直方图和正态QQ图进行正态分布检验，确定各穿透率数据均符合正态分布，不需要进行数据转换。

然后利用地统计软件GS＋9.0对穿透率数据进行半方差函数分析和参数计算，所谓半方差函数，又叫半变异函数，就是两点间插值的方差的1／2，其表达式为：

式中：γ（h）——间距h的半方差函数值；N（h）——以h为间距的穿透雨观测点对的数目；Z（xi）——样点在xi处的穿透率；Z（xi＋h）——与xi相隔间距h处样点的穿透率。下同。半变异函数值随样点间距的增加而增大，并在一定的间距（称为变程，range）达到一个基本稳定的常数（称为基台，sill）［10］。

采用GS＋9.0中不同类型半变异模型进行拟合，因相比于圆形模型、指数模型、高斯模型和线性模型而言，球面模型拟合效果最好，故本文所用半变异函数模型均为球面模型［11］。

式中：C0——块金，即间距为0时的半方差函数值，是由实验误差和随机因素引起的变异；C——结构方差，指空间自相关部分引起的变异；（C0＋C）——基台，表示系统内的总变异，一般情况下，基台值越大表示总的空间异质性程度越高，反之越小；a——变程，指空间自相关距，表示观测值之间的距离大于该值时是相互独立的，小于该值时存在空间自相关性［12］。

采用ArcGIS 10.0这一GIS软件平台，利用地统计学中常用的数据插值空间变异分析方法——Kriging插值法，进行穿透降雨空间分布的预测模拟分析。

Kriging法是建立在半变异函数理论和结构分析的基础上，在有限的区域内对区域化变量的取值进行无偏最优估计的一种方法，可对周围的实测值进行加权以得出未测位置的预测，公式由数据的加权总和组成：

式中：Z（si）——第i个位置处的测量值；λi——第i个位置处的测量值的未知权重；s0——预测位置；N——测量值数量。权重λi取决于测量点、预测位置的距离和预测位置周围的测量值之间空间关系的拟合模型。

使用ArcGIS 10.0的地统计分析模块中普通Kriging法的球面模型对穿透率数据进行插值，预测出样地中300多万个点的穿透率值，绘制成栅格图像数据，直观地模拟样地的降雨分布。

通过计算得到各降雨事件的穿透率平均值、标准差（S）和变异系数（Cv），作为评价不同降雨量时穿透率的平均状况和总体不均匀度的指标。

3 结果与分析

3.1 各降雨事件穿透降雨特征

在本实验观测期2011年8月至2012年7月共收集了29次降雨事件数据，总降雨量为947.3mm，总的穿透雨量和穿透率分别为742.6mm和0.783 9。通过对这29次降雨事件实测数据的分析，发现当大气降雨量小于1.9mm时，林冠截留了几乎全部降水，没有穿透雨产生，当大气降雨量大于1.9mm时，样地中100个样点平均穿透降雨量与大气降雨量紧密相关，穿透降雨量随大气降雨量的增加而增大，呈显著的线性正相关关系（图1a）。拟合的方程如下：

式中：T——穿透降雨量（mm）；P——大气降雨量（mm）。

随着降雨量的增加，穿透率也呈增加的趋势。对穿透率与降雨量的关系进行回归分析，比较后得出对数函数具有较好的拟合性（图1b），其关系式为：

式中：Tr——穿透率；P——大气降雨量（mm）。

图1 穿透雨量、穿透率及穿透雨量变异系数与大气降雨量的关系

可见，穿透率随降雨量的变化基本分为3个阶段：在小雨雨量级（单场降雨量10mm以下），穿透率随着降雨量的增加而迅速增大，为快变期；在中雨雨量级（单场降雨量10～25mm），穿透率增加速率明显变缓，为渐变期；在大雨雨量级（单场降雨量25～60mm）和暴雨雨量级（单场降雨量60mm以上），穿透率逐渐趋于平稳，波动范围在0.793 9～0.907 5，为稳定期。表1的穿透率平均数也可以印证这一结论。

穿透降雨变异系数（降雨事件的穿透降雨标准差与穿透降雨量的比值）揭示穿透降雨的空间分布不均匀程度，其随着降雨量的增加而呈下降的趋势。对穿透降雨变异系数与降雨量的关系进行回归分析，经过比较得出Mitscherlich模型具有较好的拟合性（图1c），Mitscherlich模型是Richards模型的一个特例，最初是描述生物生长量的模型。其关系式为：

式中：Tv——穿透降雨变异系数；P——大气降雨量（mm）。

穿透降雨变异系数随降雨量的变化也可分为快变期、渐变期和稳定期3个阶段，但变化规律与穿透率随降雨量的变化相反。在小雨雨量级内，穿透降雨变异系数随着降雨量的增加而迅速减小，在中雨雨量级之后逐渐趋于稳定，波动范围在0.193 0～0.193 2。

表1 不同雨量级穿透率的不同指标特征

3.2 不同雨量级下穿透降雨分布格局

穿透率观测值的半方差函数分析（表2）表明：基台值（C0＋C）随雨量级的增大而减小，说明降雨量增大使得降雨穿透率在空间上趋向于均匀化；各穿透率半方差函数模型空间自相关范围（变程a）大致介于5.5～6.2m；模型残差值RSS是选择模型的主要依据，所有拟合的半变异函数模型残差值均极小，RSS越小则拟合效果越好，说明各模型拟合效果都很好；决定系数R2表示穿透率半方差模型的解释效率，即模型拟合精确度的数字表示，可见暴雨的模拟精度最高，小雨的模拟精度最低，模拟精度随雨量级的增大而升高，这可能与穿透降雨不确定度随降雨量的增加而减小有关。

表2 不同雨量级下穿透率半变异函数理论模型及其参数

表2中总体指各样点29次穿透降雨率平均值的集合，小雨、中雨、大雨、暴雨分别指各样点在不同雨量级下的穿透降雨率平均值的集合。

块金与基台的比值C0／（C0＋C）表示实验误差和随机因素引起的空间变异占系统总变异的比例。从表2可见，每组穿透率数据C0／（C0＋C）值都很小，说明实验误差和随机因素引起的空间变异很小，主要是空间自相关部分引起的变异。但C0／（C0＋C）值随雨量级的增大而增大，说明随着降雨量的增加，随机因素引起的空间异质性逐渐增加，由结构性的空间因素引起的空间异质性逐渐下降，即与树冠结构的相关性逐渐降低。

从空间自相关的角度来看，C0／（C0＋C）表示变量的空间相关性程度，如果C0／（C0＋C）＜25%，说明系统具有强烈的空间相关性，如果25%≤C0／（C0＋C）≤75%，说明系统具有中等的空间相关性，如果C0／（C0＋C）＞75%，则说明系统空间相关性很弱［13］。而表2中每组数据C0／（C0＋C）都显著小于25%，因此均表现出强烈的空间相关性。其中，小雨的空间自相关最显著。

用Kriging插值法绘制穿透降雨率在各雨量级下的空间分布图（图2），由图2中可直观地看出穿透降雨在不同雨量级的空间分布格局及其差异，穿透率的空间分布受雨量级的影响呈现出不同的特点。在小雨条件下，穿透率普遍偏小，其分布东北部及个别区域穿透率相对略高，其他各处较为平均，穿透率处于0.75以下（0.6～0.75）的占到样地总面积的80%以上。中雨条件下的穿透率略有增大，穿透率在整个样地中的分布比较平滑匀称，其分布更具有代表性。大雨条件下，穿透率整体较高，但从表2来看，其分布的不均匀性并不很强，变化相对比较平缓，85%以上的区域穿透率在0.8～0.95。相对来说，在暴雨条件下的穿透降雨分布呈现出一些差异，其分布的不均匀性更明显更突出。但在不同雨量级下，穿透降雨率分布形态也有不同程度的相似，如东北部均存在一个相对高穿透率区域，中部偏西及西南角的同一位置均存在相对低穿透率区域。因为除降雨量大小外，林冠下穿透降雨的空间异质性同时受多个其他因素综合影响，包括叶面积指数、冠层厚度、离树干距离、冠层郁闭度、冠形、雨强、风速、风向、坡度等。

图2 基于100个样点在不同雨量级下的穿透降雨率平均值预测的穿透降雨率在样地内分布

3.3 穿透率实测数据与预测数据的比较

由表1可知，实测与Kriging预测的穿透率的平均数差别很小，可以看成是一致的，但是实测数据与预测数据的标准差和变异系数差别很大，预测数据的标准差和变异系数明显小于实测数据。且根据表1实测数据可知，雨量级从小雨增大时，穿透降雨率的变异系数减小。由此可知，穿透降雨分布的不均匀性随雨量级的增加而下降。但预测数据不符合上述规律，在预测数据中，大雨的标准差和变异系数最小，其次是小雨以及中雨，暴雨的标准差和变异系数最大，且其远高于小雨、中雨和大雨，小雨、中雨和大雨的标准差和变异系数差别不大。由此推断，除大雨外，穿透降雨分布的不均匀性随雨量级的增加而增大。实测数据与预测数据产生矛盾。

虽然图3中预测值和实测值频数分布的基本趋势是一致的，但预测值的最大值比实测最大值小，而预测值的最小值比实测最小值大，这意味着预测值比实测值的极差缩小。因此，实测数据频数分布图相对矮宽，预测数据的频数分布图相对高窄，Kriging预测数据的离散度下降，分布更为集中，于是预测值的分布不能完全再现样本分布，即存在“平滑效应”。李超等［14］对平滑效应的产生原因做了详尽阐释，平滑效应的存在会导致预测值的空间连续程度增强，变异程度下降。

尽管Kriging预测值存在条件偏差，但由于在Kriging方程组中对无偏性进行了强制约束，使其仍然是全局最优的。对于平滑效应的认识也要全面、客观，要针对具体的研究目的而定。

虽然平滑效应导致了以上的问题，但正是Kriging预测值存在平滑效应才使得Kriging法能够用来绘制穿透率的等值线和趋势面，对预测结果的离散程度没有任何约束的插值方法是无法绘制等值线或等值面的。

因此，Kriging预测数据在反映总体分布规律时数据不一定精确，但可以表示相对大小的趋势，统计分析应以实测数据为主。同时，对插值结果的可靠性可以用Kriging方差来指示，Kriging方差越大说明结果的可靠性越差，否则结果就越可靠［15］。可见大雨条件下的穿透率插值结果最可靠，暴雨条件下的最不可靠。

图3 不同雨量级穿透率频数分布

4 讨论

本试验在观测期内，降水总量为947.3mm，总穿透雨量和穿透率分别为742.6mm和0.7839。穿透降雨量随大气降雨量的增加而增大，呈显著线性正相关关系。很多相似研究从不同尺度、不同林分、不同气候等角度均认为线形方程能较好地拟合穿透雨量和降雨量之间的关系［16－18］。

穿透率随着大气降雨量的增加而逐渐增大，最后趋于稳定值，对数函数可以较好地描述穿透率随降雨量的变化。这一结果与翟杰等［19］对紫金山麻栎的研究，鲍文等［20］对岷江上游油松的研究结果相一致，而刘章文等［21］对祁连山灌丛的研究认为指数函数拟合效果更好。实际上刘章文等所采用的是Mitscherlich模型，采用Mitscherlich模型拟合本实验数据的穿透率与降雨量的方程为：

可见Mitscherlich模型与对数函数拟合精度完全相同（R2＝0.439）。将 Mitscherlich模型式（7）与对数函数公式（5）中的Tr用T／P代替，将其分别经过转换得到穿透雨量与降雨量的关系公式：

以上两式与线性方程（4）拟合精度完全相同（R2＝0.994），表明其与线性方程均可描述穿透雨量与降雨量的关系。因而可知，本实验采用对数函数与Mitscherlich模型也可较好拟合穿透率与降雨量之间的关系。李振新等［22］认为穿透雨率同林外降雨量之间的关系可以用logistic曲线方程模拟的结论不适用于本研究，模拟精度相对稍差（R2＝0.428），但其穿透率存在一个最大值的特点与Mitscherlich模型相似。

穿透降雨变异系数与降雨量呈负相关关系，Mitscherlich模型可以较好地描述穿透降雨变异系数随降雨量的变化。说明穿透降雨空间分布的不均匀度随大气降雨量的增加有下降的趋势，最终逐渐趋于稳定值。这与战 伟 庆 等［23］、Rodrigo A 等［24］的 研 究观点一致，但所用拟合模型不同，他们分别认为指数模型和对数模型较为合适。战伟庆［23］认为造成这一现象的原因是，当林冠达到饱和持水量时，林外降雨几乎全部转化为穿透雨。而我们认为此时林外降雨并未全部转化为穿透雨，仍然有树体吸收和蒸发等过程发生，只是林外降雨转化为穿透雨的比例趋于稳定。

5 结论

穿透降雨率的半变异函数分析表明，天然针阔混交林冠下穿透降雨率的不均匀度随雨量级的增大而减弱，空间分布趋向于均匀化。且随着雨量级的增大，随机因素对穿透率的影响相对增加，而林冠结构等结构性因素对穿透降雨率产生的影响有逐渐变小的趋势，但结构性因素始终是最主要的影响因素。从穿透率的插值分布图可以直观地看到穿透率的空间分布在不同雨量级存在着明显的差别。

由于Kriging预测数据存在“平滑效应”，使得预测值的分布不能完全体现实际分布，故Kriging预测数据只用来反映总体分布规律，数据不精确，统计分析应以实测数据为主。

天然次生针阔混交林下穿透雨空间分布的差异对林地的土壤水分分布和养分循环利用及生物多样性必然会产生重要影响，有待于进一步研究。

［1］ Anna A，Anselm R.Trace metal fluxes in bulk deposition，through－fall and stem－flow at two evergreen oak stands in NE Spain subject to different exposure to the industrial environment［J］.Atmospheric Environment，2004（38）：171－180.

［2］ 张焜.重庆四面山4种类型天然林水文功能研究［D］.北京：北京林业大学，2012.

［3］ 王波，张洪江，杜士才，等.三峡库区天然次生林凋落物森林水文效应研究［J］.水土保持通报，2009，29（3）：83－87.

［4］ 孔维健，周本智，安艳飞，等.天然次生林和人工毛竹林水文生态特征比较［J］.水土保持研究，2010，17（1）：113－116.

［5］ 刘小林，郑子龙，蔺岩雄，等.甘肃小陇山林区主要林分类型土壤水分物理性质研究［J］.西北林学院学报，2013，28（1）：7－11.

［6］ Sun Ge，Mcnulty S G，Amaty A D M，et al.A comparison of the watershed hydrology of coastal forested wetlands and the mountainous uplands in the Southern US［J］.Journal of Hydrology，2002，263（1／4）：92－104.

［7］ 马熙渊，马立琦.小陇山林区主要林分对降水的再分配［J］.甘肃林业科技，2003，28（2）：1－4.

［8］ 谭俊磊，马明国，车涛，等.基于不同郁闭度的青海云杉冠层截留特征研究［J］.地球科学进展，2009，24（7）：825－833.

［9］ 时忠杰，王彦辉，熊伟，等.单株华北落叶松树冠穿透降雨的空间异质性［J］.生态学报，2006，26（9）：2877－2886.

［10］ 史利江.基于GIS和地统计学的土壤养分空间变异特征研究［D］.上海：上海师范大学，2006.

［11］ Esr.ArcGIS Resources［OL］.［2013－03－10］.http：∥resources.arcgis.com／en／help／previous－help／index.html.

［12］ 姜城.不同经营体制下土壤养分空间变异规律及管理技术的研究［D］.北京：中国农业科学院研究生院，2000.

［13］ 王丽霞，段文标，陈立新，等.红松阔叶混交林林隙大小对土壤水分空间异质性的影响［J］.应用生态学报，2013，24（1）：17－24.

［14］ 李超.土壤水分的空间变异特性分析和农田干旱评估［D］.北京：清华大学，2009.

［15］ Yamamoto J K.An alternative measure of the reliability of ordinary Kriging estimates［J］.Mathematical Geology，2000，32（4）：489－509.

［16］ 胡珊珊，于静洁，胡堃，等.华北石质山区油松林对降水再分配过程的影响［J］.生态学报，2010，30（7）：1751－1757.

［17］ 刘世海，余新晓.北京密云水库库区水源涵养林冠层水文特征研究［J］.林业科学，2005，41（1）：194－197.

［18］ 薛建辉，郝奇林，吴永波，等.3种亚高山森林群落林冠截留量及穿透雨量与降雨量的关系［J］.南京林业大学学报：自然科学版，2008，32（3）：9－13.

［19］ 翟杰，张志民，胡海波，等.紫金山麻栎林降水分配格局研究［J］.林业科技开发，2011，25（4）：63－67.

［20］ 鲍文，包维楷，何丙辉，等.岷江上游油松人工林对降水的截留分配效应［J］.北京林业大学学报，2004，26（5）：10－16.

［21］ 刘章文，陈仁升，宋耀选，等.祁连山典型灌丛降雨截留特征［J］.生态学报，2012，32（4）：1337－1346.

［22］ 李振新，郑华，欧阳志云，等.岷江冷杉针叶林下穿透雨空间分布特征［J］.生态学报，2004，24（5）：1015－1021.

［23］ 战伟庆，张志强，武军，等.华北油松人工林冠层穿透雨空间变异性研究［J］.中国水土保持科学，2006，4（3）：26－30.

［24］ Rodrigo A，Avila A.Influence of sampling size in the estimation of mean throughfall in two Mediterranean holm oak forests［J］.Journal of Hydrology，2001，243（3／4）：216－227.