流域土壤重金属分布式运移模型研究

王圣伟 冀 豪 王 苗 娄天泷 张 畅

(西北师范大学计算机科学与工程学院， 兰州 730070)

0 引言

我国各流域尤其是中小流域土壤都面临着程度不一的水资源环境污染问题，对农业及经济发展造成了深远的影响。借助GIS技术合理估算流域土壤重金属运移量，能够对农业生产生态环境实现总体目标控制。国内外学者对于重金属污染运移的预测模型研究较为关注。RODE等[1]针对水环境的未来发展趋势，提出了一种综合水质模型，该模型的模拟结果表明，在河流较宽的部分，重金属积聚在通道边缘附近，并且由于侵蚀作用可能在未来会重新迁移。Al-MUR等[2]分析了吉达附近红海4个沉积物岩芯中重金属的分布，结果表明，在红海中重金属Mn、Cu和Pb具有较高的含量，核心沉积物中的重金属浓度在吉达中部附近增加，并且近年来已经逐步升高。CHO等[3]基于Clark单位线图及其空间分解方法的组合概念，并结合精细的空间变量流动力学，实现了空间分布的水文模拟，开发了一种集总概念和分布式特征混合水文模型。MEDIERO等[4]提出了一种基于分布式水文模型参数概率表示的标定方法，该方法用于曼萨纳雷斯河流域的案例研究，以校正后的模型作为实时洪水预报的决策支持工具。UNDUCHE等[5]开发了4种水文模型(WATFLOOD、HBV-EC、HSPF和HEC-HMS)，从复杂的完全分布模型到相对简单的集总模型，并对其有效性和准确性进行了校准和测试，以作为洪水预报的实用工具。这4种模型均能准确地模拟流量，可用于洪水预报，但是WATFLOOD和HEC-HMS模型在模拟高流量时表现不佳，HSPF和HBV-EC模型在大多数情况下对低流量预测过高，且这4种模型仅应用于草原流域的流量预测，对其他非草原流域的预测结果有待考证。

本研究以大夏河流域及其子流域的土壤和水体中重金属为切入点，对重金属的含量和迁移特征进行统计和分析，根据大夏河流域重金属的总体分布情况，选取5种易造成污染的重金属As、Cd、Cr、Cu和Pb进行分析，采用GIS图像处理技术和数据分析方法对流域遥感和水文数据进行处理，结合大夏河流域地貌形态特征和流域DEM数据进行分形计算，提出流域土壤重金属分布式运移模型，对大夏河流域土壤中重金属的运移量进行估算。

1 材料与方法

1.1 研究区概况

大夏河属于黄河上游的一级支流，地处青藏高原与黄土高原过渡带(东经102°2′～103°23′，北纬34°51′～35°48′)。大夏河经夏河县城东北流入，出土门关后进入临夏盆地，自西向东横贯临夏市、至折桥转北向东乡县泄湖峡，在康家湾喇嘛川塔张村注入刘家峡水库[6]。主河道总长203 km，总流域面积7 152 km2，多年平均流量34.3 m3/s，年均径流量11亿m3，多年含沙量3.49 kg/m3，年输沙量4.192×1010kg，林地面积3.27×106hm2，占流域总面积的21.2%[7]。

1.2 流域分级与划分

根据大夏河流域边界的经纬度，从地理空间数据云上高级检索，下载包含大夏河流域的DEM图像数据(共4块)，将下载的原始DEM图像导入ArcMAP软件，需要应用ArcMAP软件中数据管理工具—栅格—栅格数据集—镶嵌工具，把4块分散的栅格数据图合并到一起，然后根据大夏河的流域边界把大夏河流域从4块合并的栅格数据上裁剪，如图1所示。

图1 大夏河高程图Fig.1 Digital elevation model of Daxia river

本文首先对大夏河流域选取合适的阈值[8-9]进行河网水系提取，然后根据Strahler和Shreve这两种流域分级法[10]，并结合ArcMAP软件Spatial Analyst—水文分析—河网分级工具，对大夏河流域进行分级处理。而流域的划分首先需要确定各个子流域的出水口位置，然后根据水流流向[11-12]判断流入该出水口的所有上游区域，即为子流域，在ArcGIS软件中把大夏河流域划分成7个子流域，如图2所示。

图2 大夏河子流域区域划分Fig.2 Regional division of sub-basins of Daxia river basin

1.3 数据处理与分析

2016—2018年，使用手持式GPS记录位置，采集大夏河流域周边土样和水样，以全年存在稳定径流为依据，共26个采样点，采样位置具有代表性，见图3。在野外将样品放入自封袋中，编号后带回实验室，对土壤样品进行四分法取样后风干，研磨过200目尼龙筛，放入对应编号的样品容器内，对流域水样摇匀，也装入对应编号的样品容器内，大约40 mL[13]。样品以及标准样送往中国科学院兰州化学物理研究所进行检测，运用电感耦合等离子体发射光谱仪ICP-OES测定重金属的含量[14]，As、Cd、Cr、Cu、Pb 5种重金属的检出限分别为0.42、0.012、0.018、0.014、0.42 mg/L。

图3 大夏河流域采样点分布图Fig.3 Distribution map of sampling points in Daxia river basin

1.4 重金属运移模型

1.4.1SCS-CN运移模型

流域污染物的运移研究直接影响环境科学的进展，通常采用数学模型预测污染物的运移规律[15]。掌握流域重金属运移规律可以有效地控制与治理环境污染[16]。国外对于城市小流域进行地表径流污染特征的提取应用较为广泛[17]，国内相关资料较少。流域土壤的质地对重金属的迁移有一定的影响[18]。SCS-CN模型由美国农业部土壤保持局(USDA SCS)提出，由于对输入数据量要求不高、模拟精度高而被许多国家和地区广泛应用[19]。本文选取大夏河流域及其支流，以As、Cd、Cr、Cu、Pb 5种重金属为估算对象，结合不同土地利用类型，统计流域降雨径流的SCS-CN模型因子，计算不同土地利用类型下的地表径流量，估算大夏河流域的重金属运移能力[20]。

流域重金属地表径流运移量公式为

(1)

其中

S=25 400/CN-254

式中γi——第i类重金属的运移量，t/a

αi——第i类重金属的水溶出率

Aj——第j个水文响应单元(HRU)的面积，km2

P——研究区域年降雨量，mm/a

CN——径流曲线数

S——径流开始前的潜在最大滞留量，mm/a

根据大夏河流域多年降雨量，把大夏河流域按照子流域进行划分，统计各个子流域的年降雨量，并根据子流域划分图，确定子流域的流域面积，如表1所示。

表1 子流域统计结果Tab.1 Sub-basin statistical results

1.4.2分布式运移模型

定义河网分形维数与地表径流量的函数关系，可以得到大夏河各个子流域用分维值预测的径流量[21]。考虑到重金属运移量模型的单位换算，应用径流深这一概念消除单位不统一的影响，径流深指在某一时间段内通过流域指定断面的径流量除以该断面以上的流域面积的值[22]，公式为

(2)

式中R——径流深，mm

Q——年径流量，m3

F——该断面以上的流域面积，km2

根据分维值求得的年径流量，将流域重金属运移量公式定义为

(3)

式中Sj——第j个子流域的面积，km2

考虑到径流深的特点，将重金属运移模型公式修改为

(4)

式中Qj——第j个流域分维值模拟的径流量，m3

1.5 模型因子

1.5.1土地利用类型分类

在ArcGIS 10.4下，大夏河流域的土地利用类型经过遥感监督分类解析，主要选取了旱地、林地、草地、水域、居民用地和沼泽地，如图4所示。

图4 土地利用类型图Fig.4 Land use type map

1.5.2径流曲线数

径流曲线数CN取决于研究区域的土地利用类型和土壤类型，本文参照《美国国家工程手册》中的CN值[23]，结合流域多年降雨量和土壤湿润程度，确定CN值，如表2所示。参照美国水土保持局手册，选定CN为A类土壤组。

1.5.3网格法

应用网格法计算大夏河及子流域的分维数，网格法基本原理为使用不同边长的正方形网格覆盖被

表2 不同土地利用类型下的CNTab.2 Value of CN under different land use types

测流域等线体，当改变正方形的边长时，被测线体的网格数必然会出现相应的变化[24-25]，因此，拟合方程式为

lnN(a)=-Dlna+C

(5)

图5 相似流域图像Fig.5 Similar watershed image

式中a——边长N(a)——网格数

C——比例常数

D——流域的分维数

1.5.4地势起伏比的计算

径流是流域地貌形成的外营力之一，并且参与地壳中的地球化学过程，不仅影响植物的生长和湖泊、沼泽的形成，还影响土壤的发育程度。流域分维值又能够很好地反映流域地貌形态[26]，应用流域分维值拟合径流存在可行性。应用ArcGIS计算和大夏河流域形态相似的流域水文站的分维值，并计算该水文站的径流量。

由于研究区域所处的地势起伏度比较大，对地表径流量的影响也比较大，所以需要从ArcGIS中获取该研究区域最大与最小高程，然后计算地势起伏比，从而将分维值乘以地势起伏比进行归一化，公式为

(6)

式中Rm——地势起伏比

Hmax——研究区域的最大高程

Hmin——研究区域的最小高程

1.5.5相似流域的识别模型

1.5.5.1结构相似性

先分流域生成图像的像素分布函数，再通过分布函数的相似性，得到两幅图像的相似度。若用1表示图像中被点出的点，用0表示图像中未被点出的点，那么，该图像其实就是一个由1和0组成的0-1矩阵。该图像转换成0-1矩阵后可应用Matlab进行操作。因此，各个图像的结构相似问题就转换为比较矩阵的不同区域之间的相似性问题[27]。大夏河流域以及与大夏河流域比较相似的9个流域的图像如图5所示。

为了比较两幅图像的相似程度，可以定义任一图像的像素分布函数。设一幅图像F的 0-1矩阵为M。定义向量

(7)

式中vk——图像的竖直像素分布向量

mjk——矩阵M中j行k列的元素

由此，定义函数

(8)

式中k——矩阵的维数

f(k/n)——图像F的像素分布函数

根据式(8)的定义，可得到函数f(x)在x=k/n点的值。因为每幅图像的宽度为230像素，高度为230像素，所以1≤k≤230。

T=(u1,u2,…,u230)

(9)

式中T——向量F的特征向量

uk——函数f(x)在区间(0.01(k-1),0.01k)的平均值

假设两幅图像为F1和F2，它们对应的0-1 矩阵为M1、M2，对应的特征向量分别为T1、T2，由此定义

(10)

式中 〈·,·〉——两个向量的内积

α——两幅图像的结构相似度

如果α大于0.99，说明两幅图像结构相似。

1.5.5.2特征相似性

相似流域的识别不仅在结构上比较相似，在流域面积、河网长度、年平均气温、年平均降水量等某些特征指标方面也具有相似性[28]。表3对比了大夏河流域及与大夏河相似的9个流域的一些特征参数。

表3 相似流域的特征对比Tab.3 Characteristic comparison of similar basins

2 结果与分析

2.1 水溶出率

采用BCR连续分级提取法[29]，将水溶态和可交换态合并的酸溶态含量(质量比)作为水溶出率计算因子，在提取过程中可保持稳定，准确地描述污染物的运移能力。以5组土壤重金属样本作为测试提取对象，得到平均提取量和水溶出率如表4所示。

表4 BCR提取水溶态含量和水溶出率结果Tab.4 BCR extraction results of water soluble content state and water dissolution rate

2.2 重金属含量特征

由于流域位于甘肃省境内，因此采用按行政区域划分的甘肃省土壤背景值的几何平均值作为参考[30]，并对数据进行描述性统计分析，结果如表5所示。

表5 重金属含量描述性统计分析Tab.5 Descriptive statistical analysis of heavy metal content

从表5可以看出，5种重金属元素的平均含量从大到小依次为：Cr、As、Cu、Pb、Cd。

2.3 分维值模拟径流量结果

阈值对应的是区域的汇流量，也就是汇入该区域的栅格流量数，本文在ArcGIS 10.4中对该水系的矢量图取5 000的阈值，并对该矢量图进行要素转栅格、网格分析，在输入像元大小(正方形网格边长)里填入10～500不同的10组数值，得到不同边长所对应的栅格图，统计栅格图的网格数，对求得的数值取对数进行线性拟合，生成拟合方程为lnN(a)=-1.016 3lna+14.331，决定系数R2=0.999 8，拟合方程系数的斜率即为流域的分维值，故大夏河流域的分维值为1.016 3，小于1.6，判定大夏河流域地貌处于幼年期，此阶段流域水系发育不充分，地面较为完整。同理，分别计算7个子流域的分维值。

与大夏河流域比较相似的9个流域，其分维值基本与大夏河各子流域分维值相近，且某些特征指标也较为相近，所以求得相似流域分维值与径流量如表6所示。

对这9个水文站的归一化分维值和径流量做线性拟合，得到拟合曲线如图6所示。

由图6可知，相似流域分维值和径流量拟合方程为：y=-8.3×109x+6.1×109，决定系数R2=0.761 9，说明流域分维值和径流量相关性比较显著，该研究方法具有可行性。

表6 相似流域分维值和径流量Tab.6 Similar basin fractal dimension and runoff

图6 归一化的分维值和径流量拟合曲线Fig.6 Normalized fractal dimension and runoff fitting curve

根据拟合曲线求大夏河各个子流域的径流量，如表7所示。

从全国水雨情网站获取水文站的流量数据，以及从甘肃省水利厅官网下载《甘肃省水资源公报》，查询到折桥和夏河水文站的平均年径流量分别为7.265亿m3和3.817亿m3(2016年和2017年)，模拟的年径流量与实测数据基本吻合，说明用分维值模拟的径流数据可以真实反映流域径流量。

表7 大夏河流域模拟径流量Tab.7 Simulating runoff in Daxia river basin

2.4 流域相似度的计算

先对大夏河流域以及与大夏河流域比较相似的9个流域的图像进行旋转与裁剪，使每幅流域图像的宽度为230像素，高度为230像素。进而在Matlab R2014a中，把各个流域的图像导入并转换成0-1矩阵，再把数据导出，最后通过相似流域的识别模型进行计算。

根据相似流域的识别模型可以求得结构相似度α，如表8所示。

表8 流域相似度Tab.8 Watershed similarity

从表8可以看出，多个流域图块与参考图的结构相似度α均大于0.99(其中相似度为1即参考图所在的区域，是参考图与自身相比的结果)，说明大夏河流域以及与大夏河流域比较相似的9个流域具有结构相似性。

通过流域相似度的验证与分析，大夏河流域以及与大夏河流域相似的9个流域在结构上是相似的，并且根据表3可以看出，这些流域在不同尺度下的年平均气温、年平均降水量、流域面积与河网长度比值也具有相似性。

2.5 两种模型对比

根据SCS-CN模型估算出大夏河流域的重金属运移量，如表9所示。

表9 基于SCS-CN模型估算的重金属运移量Tab.9 Heavy metal migration estimated based on SCS-CN model t/a

SCS-CN模型的优点在于结构简单、输入参数少，缺点在于数据需要长期持续性观测。从表9可以看出，大夏河各个子流域的重金属运移量，As、Cd、Cr、Cu、Pb 5种重金属的总运移量为0.626 8、0.197 3、5.805 1、1.596 0、1.714 1 t/a。

根据分布式预测模型估算出大夏河流域的重金属运移量，如表10所示。

分布式运移模型优点在于从机理上去除了重金属河道滞留量和土壤入渗保留量，可以结合少量采样数据对远程遥感进行估算，缺点在于需要根据实际区域的土壤颗粒类型和水土保持因子来做进一步的剔除。从表10可以看出，大夏河各个子流域的重金属运移量，As、Cd、Cr、Cu、Pb 5种重金属的总运移量为0.671 0、0.209 9、6.281 6、1.746 5、1.837 7 t/a。

表10 基于分布式模型估算的重金属运移量Tab.10 Heavy metal migration based on distributed model estimation t/a

对比表9、10可以看出，应用分维值估算的径流量分布式模型求取的重金属运移量和应用SCS-CN模型求取的重金属运移量结果基本一致，说明分布式重金属运移模型相对合理。

3 结论

(1) Cu、Pb元素的实测平均含量分别为20.25、15.37 mg/kg，分别低于甘肃省Cu、Pb含量的地方背景值22.5、17.9 mg/kg，说明这两种元素没有超标。As、Cr、Cd元素的实测平均含量分别为20.43、94.48、1.43 mg/kg，分别高于甘肃省As、Cr、Cd含量的地方背景值11.7、69.3、0.110 6 mg/kg，说明这3种元素在大夏河流域存在不同程度的历史累积。

(2)采用网格法对大夏河流域及其子流域进行分析，其水系分维值为1.016 3。根据相似流域水文站的分维值和径流数据，计算分维值与径流量之间的线性关系，拟合曲线求出大夏河子流域的径流数据，对大夏河流域及与大夏河流域比较相似的9个流域的相似性进行了建模验证。

(3)基于SCS-CN和分布式运移两种模型的估算结果基本一致，但总体上应用分布式运移模型估算重金属的运移量偏高，比例较为合理，尤其对于远程遥感估算更具有实际应用前景。