微地形空间离散化水文响应单元确定方法

张明浩, 张 鑫, 王天宏, 徐仪玮, 张青峰

(1.西北农林科技大学 资源环境学院, 陕西 杨凌 712100; 2.杨凌职业技术学院, 陕西 杨凌 712100)

水蚀是土壤及其母质在降雨、径流等水力作用下进行的非常复杂的多重尺度地理界面过程[1]，对其高度定量化一直是地学科学家们的目标。从目前的科技水平来看，分布式参数模型因具有较高的空间分辨率依然是水蚀研究的主要方法。

微地形是指较小面积范围内，地表相对高程变化不大(通常不超过5～25 cm)的一种起伏地表[2]，它不仅是侵蚀发育的场所和直接表现形式，也是影响侵蚀发育和发展的重要因素[3]。以往侵蚀定量化研究多集中在流域、区域尺度水平[4-6]，格网多为cm级以上精度，所采用的DEM数据多具有各向同性特征。当研究尺度降至微地形尺度，DEM则表现出各向异性的空间异质性特征[7]。这也就使得以往的一些研究成果无法在微地形尺度上加以推广应用，因而难以指导生产和加深人们对水蚀发育机理的深入认识，有关微地形水蚀过程、特征和机理的研究仍显得非常薄弱。

开展微地形空间离散化(Spatial Discretization)研究，是将水蚀微地形分割成模型运行的基本地域单元[8]，根据水分运动的数学方程、边界条件以及初始条件建立相邻格网单元间的时空关系，确定各水蚀发育时期(溅蚀—片蚀—细沟侵蚀)的最佳水文响应单元(HRU，Hydrologial Response Units)，即运算时对研究区划分的最小水文单元，这也是实现水蚀分布式参数模型模拟的充要条件[9]。

元胞自动机(CA，CellularAutomata)是一种基于微观个体的相互作用，时间、空间、状态都离散的格网动力学模型[10]，可用来较好地模拟空间复杂系统的时空动态演变过程[11-12]。因其构建没有固定的数学公式，可根据不同研究需要构建不同类别的模型。

基于此，本文基于CA模型的微地形水沙汇集和传递关系分析产沙量，并与实测值进行比较，以确定最佳离散格网单元，以提出微地形空间离散化方法。

1 试验设计

降雨试验在黄土高原土壤侵蚀与旱地农业国家重点实验室进行。供试土壤为陕西省杨凌区坡耕地表层土(0—20 cm)。杨凌区(107°59′—108°08′E，34°14′—34°20′N)位于陕西关中平原中部、黄土高原南部，总面积135 km2，属温带半湿润大陆性季风气候，降水分布不均，多集中在夏秋两季，年均气温12.9℃，年均降水量和蒸发量分别为635 mm和1400 mm。

所采集的土样经风干后过5 mm筛，分5层填装在2 m×1 m×0.5 m的侵蚀槽中。每层10 cm，且土壤容重保持在1.30 g/cm3，含水率为10%，每层装土后耙磨整平表面，以保证土壤结构的均匀性和连续性。之后，由当地农民布设黄土区常见的人工掏挖(AD-Artificial Digging，从坡底逐渐向坡顶锄挖，形成空间上不称且凹凸相间的洼地)微地形表面(图1)。试验选择5个典型坡度(5°，10°，15°，20°和25°)，在90 mm/h雨强下进行分段(溅蚀阶段、片蚀阶段、细沟侵蚀阶段)降雨试验。试验中，降雨系统采用下喷式喷头，喷头距离地面18 m高，确保雨滴下落时到达地面前已达到终点速度，降水采用自来水。试验设置3组重复对照，每组试验降雨历时50 min。同时，在产生径流后每30 s收集径流，称重、静置、分离上清液，将泥沙样品置于105℃的烘箱内烘干后再次称重，得到产沙量数据，取3次试验的平均值。

利用三维激光扫描仪对分段降雨前后的微地形坡面进行扫描以得到地表点云数据，在ArcGIS中进行插值生成微地形数字高程模型(M-DEM)，其中M-DEM的原始分辨率为4 mm，具体方法参照文献[13]。

2 研究方法

2.1 CA模型

CA模型最基本的组成单元为元胞、元胞空间、邻居以及规则[14]，通过NetLogo平台编写程序代码得以实现[15]。基于本研究，CA模型的数学表达式可确定为：

CA={d,q,s,n′,j,Δt,D8}

(1)

式中:d为元胞边长(mm)；q为水流的交换量(ml)；s为水流携带泥沙的交换量(g)；n′为坡面糙度系数；j为平均坡度(°)；Δt为时间步长(s)；D8为坡面汇流规则。

在实际的降雨—径流过程中，水流由高处向低处运动。在CA模型中，水流方向指水流离开单位元胞时的指向，通过中心元胞与邻域元胞的最大距离权落差而确定[16]。本研究中，坡面汇流模型采用8方向的摩尔(Moore)型构建规则，即D8算法(中心元胞的8个邻域栅格编码)来模拟元胞间水流的流向。此外，在模型运行过程中，如果邻域元胞中存在两个及以上相同最低高程的元胞时，则流向随机选择。

水流交换量q可以借助水力学原理，采用曼宁公式得出[17]，即：

(2)

式中:h为两个元胞之间的水深差(mm);n为曼宁糙度系数。

在坡面水蚀过程中，当径流作用于土壤的力大于土壤阻力时，土粒会随着水流一起运动[18]，称之为泥沙的输移过程。水流携带泥沙的量s由泥沙运动的经验公式求得：

s=kquα

(3)

式中:u为水流的流速(mm/s)；k和α为水流的携沙系数,参考有关水力学的研究,结果取值分别为0.05,0.5[19]。

基于CA的微地形水沙汇集和传递关系严重影响着坡面的径流量和产沙量。本研究中，任意一个Δt内，元胞空间中各元胞自身的产沙量由本元胞空间降雨产生的流量与来自邻域元胞的流量的共同携沙量组成。CA模型假设水蚀模拟过程中水流从上游元胞流向下游元胞的同时对泥沙进行运输，且这一过程中无泥沙沉积，即上游元胞的泥沙均被运输到下游元胞中，则单位时间内微地形坡面出口断面累积产沙量就等于该单位时间内元胞径流量逐级随流向汇集到出口断面的累积携沙量。

在坡面水蚀过程中，影响产沙量的因素众多，其中较为重要的是坡面糙度。坡面糙度可通过区域内地表的表面积与其投影面积之比来计算，它也是反映地表形态的一个宏观指标[20-21]。研究表明，坡耕地的糙度在一定范围内，坡面糙度系数同曼宁糙度系数变化一致[22]。因此，本研究中曼宁糙度系数n通过计算微地形的坡面糙度获得，即：

(4)

式中:S曲面为栅格表面积；S投影为栅格的水平投影面积；γ为栅格单元的坡度。

利用ArcGIS对扫描获得的地表点云数据进行M-DEM、坡度、糙度系数等预处理(图1)，将处理后的一系列栅格文件编写到CA运行代码中，同时将元胞设置为1～10 mm(间隔为1 mm)十种大小进行CA模型模拟，并记录单位时间内的产沙量。

图1 试验区图

2.2 模拟结果验证

纳什效率系数(NSE)是用来评价模型质量的参数，常被用来验证水文模型模拟结果的好坏，用Ef表示。

(5)

均方根误差(RMSE)用来表示模拟值与试验值间的偏差，用Re表示。

(6)

式中:Re为均方根误差；X模型，i为CA模拟值；X实测，i为试验记录值；r为样本个数。在一组模拟数据中,Re越小表明模拟精度越高,反之则越低。

3 结果与分析

3.1 微地形坡面糙度的变化

对数据进行分析得出，整个试验过程中平均糙度系数整体上为1.178 9±0.06。当坡度相同时，3个侵蚀发育阶段的平均糙度系数并无较大差异；当坡度不同时，随着坡度的增加，平均糙度系数呈现增大的趋势，同时累计产沙量也呈现增大的趋势(表1)，说明粗糙地表会比光滑地表产生更多的泥沙，这种情况的发生往往归因于粗糙坡面使得径流发生集中，因而携带更多泥沙。

表1 不同水蚀阶段糙度系数与产沙量

3.2 不同格网/元胞大小对水沙传递的影响

运用CA模型对微地形不同坡度和侵蚀发育阶段进行模拟，得出150组模拟值与试验值对比分析结果如下(试验值与最佳模拟值产沙量如图2所示)：

注：单位时间为30s。

(1)不同坡度下各水蚀发育阶段所发生的时长不同，坡度越大时长越短；在相同坡度条件下，t细沟侵蚀>t片蚀>t溅蚀。随着水蚀现象的发生，累积产沙量与时间表现出公式为s=alnt±b(a，b为回归系数；s为产沙量(g)；t为降雨历时(s))的对数关系，且R2> 0.60，这与前人研究径流量与时间的关系[24]一致；

(2)从1～10 mm不同大小格网/元胞的模拟结果可以看出，并非格网/元胞越小模拟效果越好，过小时其模拟值与试验值产生较大偏差，拟合性总是呈现先增后降的趋势，模拟值与试验值贴合性最优时对应的格网/元胞大小即为最佳水文响应单元；

(3)对于不同水蚀发育阶段，初始阶段单位时间内水流携沙量相对来说是较大的，因为在该阶段，坡面土粒较为松散，土壤抗蚀能力弱，径流易对其进行搬运和堆积。随着水蚀的进行，到片蚀阶段时由溅蚀转化为薄层水流冲刷，一方面由于土壤表层颗粒减少，在一定程度上降低了土壤颗粒随径流流失的速度；另一方面溅蚀破坏了土壤表层结构，形成土壤结皮，阻碍了侵蚀强度的进一步增加，使水流携沙量减少。因此，片蚀阶段的产沙量逐渐变小，最终在细沟侵蚀阶段逐渐达到相对稳定的水平。

3.3 模型的验证及最佳HRU的判定

运用纳什效率系数(Ef)和均方根误差(Re)进行双重验证(图3)，可以看出：Ef值均处于较高的水平，CA模型模拟的结果较好。当Ef值达到峰值时，所对应的格网/元胞大小即为最佳水文响应单元，此时，所对应的Re值也为最佳；不同坡度下最佳Ef值所对应的格网/元胞大小和Re值见表2。

图3 不同格网/元胞大小对应的Ef及Re

表2 最佳HRU及其Ef，Re

4 讨论与结论

4.1 讨 论

研究微地形水蚀状况，对其进行空间离散化是不可缺少的基础性操作。在研究微地形水蚀时格网大小是否合理对模型模拟结果有着重要影响。M-DEM格网大小不同使得径流坡度、坡面糙度等均发生变化，从而影响水流在流域河网中的携沙量。研究表明：以M-DEM数据为基础建立水文模型时，并非M-DEM格网单元越小越好，只有当M-DEM中某一格网大小使得各个要素与流域内实际降雨过程中各要素值出现最佳拟合，模型的模拟效果才相对最好。因此在利用水文模型模拟坡面汇流时，选择最佳的水文响应单元至关重要，否则会使得模型产生较大的误差。

但不可否认的是，研究仍存在一些问题值得思考：(1)本研究内容仅考虑了单一雨强和耕作措施，对不同耕作措施、雨强下水蚀状况的模拟有待研究；(2)对于本模型中的一些假设，如坡面汇流采用D8算法，是否存在其他坡面汇流方式仍有待进一步探究；(3)CA模型的假设“水流从上游元胞流向下游元胞过程中无泥沙沉积”，在实际的降雨过程中是否存在泥沙在输移时沉降的现象，及其是否对模拟结果能够产生较大影响，等等。在以后可对此进行更系统、更深入的研究；(4)空间参数化是指按照模型的具体要求，对空间离散化地域单元进行说明和定值的方法。如何在离散化单元确定的基础上开展参数化的研究，进而构建微地形水蚀分布式参数模型，是今后重要研究的内容。

4.2 结 论

(1)CA模型可用来较好地进行微地形水蚀发育过程的定量化模拟，其模拟结果和实际试验并无较大差异。

(2)论文确定了一种基于CA元胞单元的微地形尺度下的空间离散化方法。通过此模型确定的最佳水文响应单元所对应的M-DEM可作为后续研究中的基础。

(3)空间离散化格网单元尺度的选择并非随意选择，必须要小到能够反映其地理因素的空间变化，又要大到可以较为准确地获取各种输入参数、运行模型的水平，因此需根据研究对象的具体情况来进行详细的判断。通过基于微地形空间离散化的研究，初步判定在90 mm/h雨强下的最佳水文响应单元为6 mm。