喀斯特地区地表侵蚀产沙对前期有效降雨的响应

潘佑静 唐丽霞 杨 智 杜 波 罗 金

(1. 贵州省瓮安县猴场水库管理所，贵州 黔南 550400；2. 贵州大学林学院，贵州 贵阳 550025；3. 贵州省水土保持技术咨询研究中心，贵州 贵阳 550002；4. 贵州省清镇水务局，贵州 贵阳 551400；5. 贵州省瓮安县天文水库管理所，贵州 瓮安 550400)

喀斯特地区是典型的生态环境脆弱区，下垫面的高度异质性加之人为破坏使得该地区的水土流失极为严重，土壤侵蚀产沙的成因非常复杂[1]。长期的岩溶作用构成地表和地下双层空间结构[2]，降雨是造成该地区地表土壤水蚀的主要动力因素，如何确定有效的侵蚀性降雨及其对地表侵蚀产沙的影响尤为重要[3]。目前对侵蚀性降雨的研究主要集中在雨强、降雨历时及降雨量等方面[3-4]，对前期有效降雨研究较少考虑。前期有效降雨通过影响坡面土壤的含水量，改变土壤的粘性、抗剪强度、入渗初期的水势梯度等物理特性，影响土壤水文过程及土壤侵蚀。国内外对前期有效降雨的研究主要集中在有效前期降雨值的确定[5-7]，泥石流形成中的作用[8]，以及常态地貌下与流域产流产沙的关系研究[9]，相对常态地貌，前期有效降雨对喀斯特地区地表产流产沙影响更为显著，然而这方面现有的研究却很少。喀斯特地区碳酸盐岩风化发育的石灰土土层薄、土壤总量少、储水能力低，流域以蓄满产流为基本产流机制[10]，在非超渗产流机制主导下，前期土壤水分对入渗能力有着直接的影响，降雨发生时的土壤水分越高，土壤入渗能力就越弱，从而更有利于地表径流的产生[11]，产沙量也随之呈现增加的变化趋势[12-13]。

本研究以贵州省喀斯特地区的关岭县享乐小流域为研究区，通过实测的日降雨量，采用产沙累积频度-降雨量分形关系计算流域前期有效降雨量[14-15]，划分流域发生或不发生侵蚀产沙的降雨条件，从而获取喀斯特区地表侵蚀产沙的降雨条件阈值，以期为贵州省喀斯特地区防治坡面水土流失、坡耕地整治、生态恢复与重建提供参考依据。

1 研究区概况

研究区为贵州省关岭县的享乐小流域，地处东经105°32′59″～105°34′34″，北纬25°48′22″～25°50′12″，流域面积为3.93 km2。享乐小流域属长江流域乌江水系，地貌类型是以溶峰林、峰丛-洼地、谷地为主的岩溶低中山丘陵，流域内有地下水出露点与落水洞，流域岩溶面积占83.84%，水土流失以水力侵蚀为主。地势为西高东低，由西北向东南倾斜，西部最高达1 775 m，中部次之为1 470 m，东部最低点1 160 m，相对高差为615 m。流域地处亚热带湿润季风气候区，年平均气温19.2 ℃；月平均最低气温1月最低，为6.5 ℃； 8月最高，为24.1 ℃。年平均降雨量为1 236 mm，无霜期298 d。土壤主要以水稻土、黄壤、石灰土为主；现有植被类型主要为次生乔木林、疏幼林，森林覆盖率为22.42%，树种以杉木 (Cunninghamialanceolata)、柳杉 (Cryptomeriafortunei)、光皮桦 (Betulaluminifera) 为主。

2 研究方法

2.1 数据观测

于流域小型气象站和卡口站内进行观测，采用HOBO小型自动气象站记录流域历次降雨过程数据，包括降雨量，降雨开始及停止时间，降雨历时，共获取流域2009—2013年5 a间的降雨事件数据。卡口站采用矩形薄壁量水堰，径流数据采用SWZ 型周记式水位计进行记录，同时采用自动水位取样器取样，将水样用滤纸过滤，获取河道泥沙悬移质、推移质数据，计算出每次降雨事件中流域产沙量，2009—2013年共有562 d发生了降雨，其中70 d的降雨造成了侵蚀产沙，详细数据见表1。

表1 2009—2013年降雨、产沙天数统计Table 1 The statistics of rain fall and sediment during 2009-2013

2.2 分形分布拟合

分形分布广泛存在于自然界，所有规模的冲动频率都符合分形分布，其函数关系为幂律，其形式见公式 (1)。在双对数坐标下，幂律分布表现为1条斜率为负数的直线，这一线性关系是判断随机变量是否满足幂律的依据。

p(x)=Cx-D

(1)

式中：p(x)是具有特征线度为x的物体的数目；C为比例常数；D为标度指数也称分形维数或分维，是描述分形对象的空间分布格局或内在的等级层次等最基本的特征量。

采用分形分布理论拟合流域2009—2013年562次降雨事件70次产沙事件中的降雨产沙累积频度 (Y) 和降雨量 (L) 间的幂律关系。

2.3 前期有效降雨量计算

引发流域产沙的累积降雨量阈值为产沙累积降雨阈值 (Lin)，对于触发产沙的降雨，当流域累积降雨量达到Lin值，就有某一产沙累积频度的可能在该流域内引发产沙，将该产沙累积频度定义为Ycf。即当产沙累积频度为Ycf时 ，产沙前n天的累积降雨阈值减去前n-1天的累积降雨阈值即为第n天的降雨经衰减后对n-1天产沙发生影响的有效降雨量，产沙发生前多天降雨衰减后的累计有效降雨量即为流域产沙发生时的前期有效降雨量，据此可得流域产沙发生时的前期有效降雨量，方法如下：

Pean=Linn-Linn-1

(2)

(3)

式中：PEA为前期有效降雨量；Pean为第n天的降雨经衰减后对n-1天产沙发生影响的有效降雨量，Linn为产沙前n天的累积降雨阈值；Linn-1前n-1天的累积降雨阈值，m为前期降雨天数。

3 结果与分析

3.1 累积产沙频度与累积降雨量分形分析

2009—2013年期间流域内最长降雨持续历时为11 d，考虑导致产沙的降雨历时时段内至其前10天的降雨情况。以侵蚀产沙当天作为前第1天，将每次产沙事件划分前1、2、4、6、8、11天，6个时间段的产沙累积频度和降雨量绘制成底数为10的双对数坐标图 (图1)。

图1不同时段内累积产沙频度-累积降雨分布图
Fig.1 Distribution cumulative frequency of sediment yield-cumulative rainfall in different periods

从图1可以看出，6个时段的产沙累积频度随降雨量的变化均出现2个不同的无标度区间，分别按照公式 (1) 拟合得到的2条不同斜率 (标度指数) 的直线。在6个时段的产沙累积频度与降雨量分形拟合中决定系数 (R2) 最大值为0.999 4，最小值为0.782 2，均值为0.924 7，整体拟合效果较好，即表明产沙累积频度与降雨量之间遵循分形分布。

对于关岭享乐小流域的产沙发生当天至其前2、4、6、8、11天，6个时间段，不同标度区间交点处的Y值为58%～78%，平均为68.88%，整体在70%左右，当降雨量超过交点处的L值后，产沙的累积频度仅为22%～42% (平均31.12%)。该点的降雨量为引发产沙的累积降雨阈值的上边界，即在该流域引发70%左右产沙的累积降雨阈值 (Lin) 分别为31.6 mm (1 d)、37.8 mm (2 d)、41.4 mm (4 d)、51.7 mm (6 d)、69.9 mm (8 d)、77.8 mm (11 d)。

3.2 累积降雨阈值与时段分形分析

将流域各时段70%产沙可能时的累积降雨阈值 (Lin) 作为累积降雨时段天数 (n) 的函数标绘为横坐标以2为底数，纵坐标以10为底数的双对数坐标图 (图2)。

图2累积降雨阈值与累积降雨时段的分形关系
Fig.2 A power-exponential relation between cumulative rainfall threshold and cumulative rainfall time

可以发现两者之间的函数遵循幂律分布：

Lin=29.147n0.3741≤n≤11

(4)

3.3 前期有效降雨分析

将侵蚀产沙当天作为前1天，对于前n天由公式 (4) 有

Linn=29.147n0.374(1≤n≤m)

(5)

对于前n-1天有

Linn-1=29.147(n-1)0.374

(6)

联合式 (5) 有

(7)

联合式 (2)、式 (3)、式 (6)、式 (7) 得到对于前n-1天有效的降雨量Pean：

(8)

(9)

对于Lin有

(10)

式中：Pi为相对产沙当天前第i天的降雨量，当i取1时，P1为产沙当天的降雨量。

前期累积降雨对流域产沙发生的影响即为前期有效降雨量。对于流域侵蚀产沙发生时 (前第1天) 的有效前期降雨量PEA，由式 (9) 有

(11)

使用有效前期降雨量模式 (11)，计算出关岭享乐小流域2009—2013年562个降雨事件的累积前10 d前期有效降雨量。其中日降雨量最大值为102.9 mm，最小值为0.1 mm，当日前期有效降雨量为分别13 mm、44.89 mm；前期有效降雨最大为117 mm，并有多天为0 mm。不同降雨日降雨情况异质性强，降雨量、前期有效降雨量变异系数分别为1.69、1.38，降雨量的变异性大于前期有效降雨量，流域日降雨量与前期有效降雨量变化不同步，发生产沙的日降雨量最小值为13 mm，该天前期有效前期降雨量为29.33 mm。研究期内共有106 d降雨量大于13 mm，其中38 d未发生产沙，未产沙降雨日最大降雨量为47.4 mm，其前期有效降雨为16.65 mm。可知，日降雨量或前期有效降雨量单一因素不能预测流域产沙事件的发生。

3.4 侵蚀降雨产沙预测

用日降雨量与前期有效降雨量共同确定流域产沙事件的发生，将小流域2009—2012年中共452个降雨日的流域日降雨量作为前期有效降雨量的函数标绘在双对数坐标图上 (图3)。连接产沙时的最小日降雨量与最小前期有效降雨点，获得引发该流域70%以上的产沙可能的降雨阈值下边界线 (图4中的直线)，线性方程为：

P1=20.879-0.282 9PEA

(12)

图3引发产沙的最小降雨阈值
Fig.3 Minimum rainfall threshold triggers sediment yield

将阈值下边界线作为流域产沙降雨条件的判别线，当前期有效降雨量与当日降雨量点 (PEA,P1) 落在阈值线以上时，视为流域发生侵蚀性产沙，当点 (PEA,P1) 落在阈值边界线以下时流域不发生侵蚀产沙。代表降雨日的452个点中有128个落于阈值边界线上，324个落于阈值边界线下，与流域实际59 d产沙、303 d未发生产沙相比，准确率为86.95%。则有当前期有效降雨量与当日降雨量点 (PEA,P1) 落在阈值边界线以上时，该流域发生侵蚀性产沙可能性较高，当点 (PEA,P1) 落在阈值线以下时发生侵蚀产沙的可能性较低。

将关岭享乐小流域2013年发生的降雨前10天前期有效降雨量与日降雨量标绘在双对数坐标图上，并添加日降雨阈值下边界线 (图4)。最小降雨阈值以上和以下的 (PEA,P1) 点分别为9个与101个，即2013年该流域应有9次侵蚀产沙降雨与101次非侵蚀性降雨。对比发现该流域2013年实际共有110 d发生降雨，其中侵蚀产沙型降雨11次，非侵蚀型降雨99次，与预测值相差4次。在70%侵蚀产沙可能情况下，预测准确106次，预测准确率达到96.36%。对于享乐小流域，将日降雨阈值下边界线与产沙发生有关的区域平均降雨条件联合能很好的预测流域的侵蚀产沙的发生。

图42013年降雨-产沙情况
Fig.4 The relationship of rainfall and sediment in 2013

4 结论与讨论

1) 关岭享乐小流域地表侵蚀产沙事件的累积频度和降雨量、产沙事件的累积降雨阈值和降雨时段之间均遵循分形分布的幂指数关系。利用分形分布理论结合流域降雨、产沙情况可建立适用于喀斯特地区的前期有效降雨模型。

2) 基于累积产沙频度-降雨量分形关系建立的前期有效降雨模型考虑了降雨衰减系数并非经验系数，而是由流域内引发产沙的累积降雨阈值随观测时段天数变化的标度指数所决定，能去除衰减系数的取值随意性，相比传统模型可较为符合实际的计算出产沙发生时的流域前期有效降雨量。

3) 将日降雨量与有效前期降雨量结合可以确定喀斯特流域发生侵蚀性降雨产沙的降雨阈值。关岭享乐小流域基于前期有效降雨量PEA、当日降雨量P1的流域地表侵蚀产沙发生阈值线性方程为P1=20.879-0.282 9PEA，当前期有效降雨量与当日降雨量点 (PEA,P1) 落在阈值线以下不发生流域地表侵蚀产沙，落在阈值范围以上发生侵蚀产沙。该阈值线性方程对关岭享乐小流域2013年产沙事件的预测准确率达到96.36%。

4) 对喀斯特流域，通过降雨阈值根据日降雨预报即可得到有效、较为准确的流域产沙预测结果。因此，可将前期有效降雨量作为是否发生产沙的一个重要参数，对喀斯特流域侵蚀产沙发生的预测具有重要意义。

目前前期有效降雨计算方法多是是基于蒸发量、降雨量等基本指标构建的简单模型[4-7]，但是各自所取的衰减系数因各地区太热辐射、气温等小生境气候因子不同而具有较大差异，多采取区域经验值，具有随意性。采用累积产沙频度-降雨量分形关系计算该流域前期有效降雨量，能去除衰减系数的取值随意性。但因各流域的产沙情况与日降雨情况不同，其标度指数也不同，对于不同流域需各自计算其标度指数，构建前期有效降雨模型，相比传统模型计算量大且复杂，对该方法的应用有一定的难度。通过建立多个喀斯特流域的累积产沙频度-降雨量前期有效降雨模型计算更为精确、适用于喀斯特地区的经验衰减系数可进一步的研究。

基于分形分布计算自然灾害预测被较多的应用于常态地貌的滑坡、泥石流等地质灾害预测中，对于喀斯特地区鲜见。喀斯特地区由于复杂的地形地貌和地下、地表二元空间结构影响，各种环境因子之间存在复杂的相互关系，釆用常规分析方法很难准确界定各因子的影响程度，难以预测土壤侵蚀的发生，更难以有效预测滑坡、崩塌、泥石流、岩溶塌陷等自然灾害的发生。分形分布是普遍存在于自然界的，所有规模的冲动频率都符合分形分布，本研究通过利用分形分布计算出流域产沙阈值，对流域侵蚀产沙做出很好的参考，基于分形分布计算自然灾害发生前环境因子状态，也可以建立灾害发生时的影响因子自然临界阈值，对今后在喀斯特地区自然灾害的预测提供参考。

[参 考 文 献]

[1] 凡非得, 王克林, 熊鹰, 等. 西南喀斯特区域水土流失敏感性评价及其空间分异特征[J]. 生态学报, 2011, 31(21): 6353-6362.

[2] 曹建华, 袁道先, 章程, 等. 受地质条件制约的中国西南岩溶生态系统[J]. 地球与环境, 2004, 32(1): 1-8.

[3] 李瑞, 李勇, 刘云芳. 贵州喀斯特地区降雨与坡面土壤侵蚀关系研究[J]. 水土保持研究, 2012, 19(3): 7-11.

[4] 刘正堂, 戴全厚, 倪九派, 等. 喀斯特地区裸坡面土壤侵蚀的人工模拟降雨试验研究[J]. 水土保持学报, 2013, 27(5): 12-16.

[5] Glade T, Crazier M, Smith P. Applying probability determination to refine landslide-triggering rainfall thresholds using an empirical “Antecedent Daily Rainfall Model” [J]. Pure and Applied Geophysics, 2000, 157(6): 1059-1079.

[6] Zêzere J L, Trigo R M, Trigo I F. Shallow and deep landslides induced by rainfall in the Lisbon region (Portugal): assessment of relationships with the North Atlantic Oscillation[J]. Natural Hazards and Earth System Science, 2005, 5(3): 331-344.

[7] 姜志文, 刘树春. 柴河水库流域前期影响雨量的计算方法[J]. 现代农业科技, 2010(9): 258, 260.

[8] 崔鹏, 杨坤, 陈杰. 前期降雨对泥石流形成的贡献: 以蒋家沟泥石流形成为例[J]. 中国水土保持科学, 2003, 1(1): 11-15.

[9] 李昌志, 刘兴年, 曹叔尤, 等. 前期降雨与不同沙源条件小流域产沙关系的对比研究[J]. 水土保持学报, 2001, 15(2): 36-39.

[10] 陈洪松, 王克林. 西南喀斯特山区土壤水分研究[J]. 农业现代化研究, 2008, 29(6): 734-738.

[11] 沈紫燕, 王辉, 平李娜, 等. 前期土壤含水量对粘性红壤产流产沙及溶质运移的影响[J]. 水土保持学报, 2014, 28(1): 58-62.

[12] 王海斌. 不同水土保持措施下径流小区降雨与产流产沙关系研究[J]. 水土保持研究, 2011, 18(5): 63-66.

[13] 纪启芳, 张兴奇, 张科利, 等. 贵州省喀斯特地区坡面产流产沙特征[J]. 水土保持研究, 2012, 19(4): 1-5.

[14] 李长江, 麻土华, 孙乐玲, 等. 降雨型滑坡预报中计算前期有效降雨量的一种新方法[J]. 山地学报, 2011, 29(1): 81-86.

[15] 李长江, 麻土华, 朱兴盛. 降雨型滑坡预报的理论、方法及应用[M]. 北京: 地质出版社, 2008.