水资源保护：确定模型输入参数的方法

4.1.2.1　地理参数的确定

1.DEM的建立和处理

模型是基于地形变化对每个cell和reach来做模拟，它们的划分和参数的提取是整个模拟工作的基础，cell和reach数据集是AnnAGNPS模型输入数据集的一个非常重要的数据集。

首先将DEM数据调入接口程序界面，选择系统菜单“DEM Utilities”中的“Process DEM”来处理双阳河流域的DEM，Process DEM功能可根据给出的流域出口位置自动确定流网、分流域和流域出口。处理过程中会有一对话框提示输入初始化河流所需要的像元个数值，这个值可以是系统的缺省值也可以是用户根据实际需要确定的像元个数。初始化像元个数越大，计算出的河流相对较少，河网越稀，初始化像元个数越小，计算出的河流相对较多，河网越密。经过反复尝试，与实际的河网对比，确定初始化河流的像元个数为8000。完成处理过程后，视图中会添加：fillDem，FlowDir，FlowAcc，StreamGrd，LinkGrd，WatshdGrd和Outlets等图层，双阳河流域的数字化高程如图4.1所示。

栅格DEM建立的一般方法和过程为：①矢量化带高程值的地形图，并建立空间拓扑关系；②地图投影坐标系统设置，为了得到较精确的流域面积，采用高斯克里格投影，坐标系为北京54坐标，单位为m；③将矢量线图转换成TIN（Triangulated Irregular Network）数据结构，并进行空间插值；④将TIN数据结构转换成格网结构（LATTICE）；⑤对格网数据进行重采样，形成栅格DEM（GRID）。

本研究采用SRTM3 DEM原始数据。SRTM数据主要是由美国太空总署（NASA）和国防部国家测绘局（NIMA）联合测量的，SRTM的全称是Shuttle Radar Topography Mission，即航天飞机雷达地形测绘使命，2000年2月11日，美国“奋进”号上搭载的SRTM系统共计进行了222小时23分钟的数据采集工作，获取的雷达影像数据经过两年多的处理，制成了数字地形高程模型，该测量数据覆盖中国全境。

SRTM数据每经纬度方格提供一个文件，精度有1 arc-second和3 arc-seconds两种，称作SRTM1和SRTM3，或者称作30M和90M数据，SRTM1的文件里面包含3601×3601个采样点的高度数据，SRTM3的文件里面包含1201×1201个采样点的高度数据。目前能够免费获取中国境内的SRTM3文件，是90m的数据，每个90m的数据点是由9个30m的数据点算术平均得来的，在本文中采用SRTM3文件。

在ArcView-AGNPS系统中，经投影转换、图幅合并（Merge Grids）、网格重分和流域范围切边等几个步骤得到双阳河流域数字高程模型，如图4.1所示。

图4.1　双阳河流域数字高程图（DEM）

2.流域地面径流流向分析

水流方向分析是利用DEM进行流域地表其他水文特征分析的基础。由于DEM模式是由栅格像元组成，因此应计算每个栅格像元的水流流向。一个栅格像元（cell）的水流流向就是水体从其中流出的方向。水流流向的确定有多种算法，可分为单流向算法（SFD）、多流向算法（MFD）及其他算法［2］。AGNPS中的DEM Utilities水文分析模块采用单流向算法中的D8算法进行水流流向的计算。D8算法是生成水流方向最常用的算法，这种算法首先计算单元格与周围8个单元格间的坡度，按最陡坡度原则确定单元格的流向。通过栅格单元x的8个邻域栅格编码，水流方向便可以其中的某一值来确定。有效的水流方向定义为东北、东、东南、南、西南、西、西北和北，并分别用128、1、2、4、8、16、32和64这8个有效特征码表示。

3.集水区的划分与河流网络的分析

汇流能力强的栅格所在的区域就是集水区，集水能力特征值为零的则代表着区内的高地。据此可划分研究区内分水岭的空间分布和集水区的边界。根据栅格汇流能力的分析结果，将汇流能力超过一定阈值（一般大100）的栅格作为水道，将小于该阈值的栅格作为产流区，从而划分河流网络系统。河流网络栅格文件建立之后，可进一步确定河流的级序（Stream Order）和河流链等特征（Stream Link）。

激活图层“StreamGrid”，将图像放大到可以看清每个像元，以显示出流点的位置，否则难以正确确定出流点的位置。选择系统菜单“DEM Utilities”中的“Delineae Watershed”，选定一个流域出口，系统将根据高程的变化情况，生成流域的边界，经计算，该流域包含的面积是1763.71km2。完成这一步骤后，出流点出口的坐标也被记录下来。在“Data Prep”下拉菜单中选择“Assign Themes”选择程序要用到的一些图层，如果以前离开了Arcview，并已经存储了设计，现在要重新工作，就必须再次选定图层。放大到整个流域边界。用矩形选框工具圈定整个边界，激活图形。对剪切下来的这一部分DEM进行处理可以减少程序运行的时间和提高程序运行的效率。切下来的DEM命名为“DEM SUBSET”，后面的程序都是基于“DEM SUBSET”图层来计算的。

处理后的DEM如图4.1所示，生成的子流域与汇流点如图4.2所示，以及系统生成的流域边界如图4.3所示。

图4.2　由双阳河DEM生成的子流域与汇流点

图4.3　双阳河流域边界

4.1.2.2　基于AGNPS和DEM的双阳河流域地表水文模拟分析

地表水文模拟分析，具体的方法和步骤主要包括3个方面，即：①洼地（沉陷点）填充、平地处理及水流流向的分析；②水量计算及河流网络模拟分析；③流域划分。

1.洼地填充、平地处理与水流方向分析

洼地填充与平地处理是通过DEM提取流域特征点数据。在DEM生成过程中，由于数据采集的误差或地表存在凸凹，在进行水文分析时会造成产流进程的停止，因此在计算水流方向以前，必须对洼地和平地进行处理。DEM中的洼地分为凹陷型洼地和阻挡型洼地［2］；消除洼地的常用方法有平滑处理、填平抬高等方法，前者可消除孤立的、较浅的洼地，保留较大的洼地；后者可以消除所有洼地，但会产生大片平坦的地形。平地水流方向的确定也是DEM提取数字流域特征时会遇到的难题。处理方法一般为采用高程增量叠加算法设定平坦格网内的水流方向，最终使每个栅格元素（cell）都有一个明确的水流方向，从而能够产生合理的汇流水系。本例中，利用AGNPS中的DEM Utilities模块中的Neighborhood Statistics命令采用填平抬高的方法消除洼地。通过Sink函数确定DEM上的洼地，再用Fill函数将洼地填充，使洼地成为水流能够通过的平坦区域，再经过反复运算消除平地，最终可生成无洼地的数字高程模型（Filled DEM），保证了水流能畅通流至河口，如图4.4所示。在Filled DEM中由于无洼地区域存在，自然流水可以畅通无阻地流至区域地形的边缘。因此，经洼地填充之后双阳河河口就显示出来。再利用Process DEM命令便可进行水流流向分析如图4.5所示。

图4.4　洼地填充后的双阳河DEM

图4.5　生成的水流方向图

2.汇流量计算与汇流流网生成

汇入某个栅格元素的上游元素总数称为该栅格元素的汇流量。汇流量大的元素是主要水系和河流的元素，这些具有较高汇流量的栅格元素组成了水系。

在生成水系流网时，需事先设定汇流量的最小值，才能生成水网，水网密集程度与细节程度取决于给定的最小汇流量值，给定的数值越小，水系分布越密集。在水流流向分析完成后，AGNPS模型要求计算水流累积量。水流累计量由DEM Utilities模块计算。通过设定汇流量的最小阈值，如设100，1000，10000，（通过Delineate Watershed命令）计算水流累积量，可获得基于NO DATA背景的水系栅格图［4］。根据能体现双阳河主要水系，但不至于生成太多细小河流的原则，本文选定了阈值1000来模拟出水系栅格图如图4.6所示，该模拟结果与1980年出版的黑龙江省水系图中的双阳河流域的水系分布图吻合很好。反映了利用DEM模拟的河流网络的空间分布接近实况。

3.流域分水线及流域面积确定

流域是指流经其中的水流和其他物质从一个公共的出水口排水从而形成的一个集中的排水区域，是向一个给定点（汇入点）贡献汇流的区域，汇入点就是汇流集中的地方，它也是流域边界上的最低点。流域常用于定义自然资源管理区域，确定水系中发现的污染物的可能源头［5］。

图4.6　双阳河汇流量图（阈值1000）

流域分水线包围的区域面积就是流域面积。任何一个天然的河网都是由大小不等的、各种各样的水道所联合而成，而每个水道都有自己的汇水范围和流域面积，这样就形成了子流域，而较大的流域往往是由这些子流域联合组成。一个流域里的子流域个数取决于在汇流量计算时所确定的阈值，越高的阈值生成的子流域的数量就越少。本例中，在生成子流域后和各子流域的注入点后，如图4.6所示，通过AGNPS Data Prep工具里的Select Watershed Outlet命令，在水系上确定整体流域的注入点即可将所有水系最终汇入该点的所有子流域重新生成为一个流域，而流域边界即为全流域分水线如图4.7所示。

图4.7　双阳河水系汇水网络

4.双阳河全流域面积的确定

流域分水线确定后，在该流域的属性信息中可以看到流域面积（字段名为Shape_Area）由251728个栅格像元组成，由于每个栅格面积为90m×90m，计算后的流域面积为1963.19km2，此数据与黑龙江省省志（水利卷）统计的双阳河流域面积1861km2相差+5.5%。较好的模拟了双阳河流域的实际情况如图4.8所示。

图4.8　双阳河流域数字地表模拟（左）与TIN图（右）

执行“Find Outlet”寻找出流点，首先确保流域出流点位于流域下游出口处，出流点的坐标被记录在dednm.exe程序的控制文件dnmcnt.inp里。执行“create TopAGNPS Input Files”，这个菜单项将在1_TopAGNPS_DataSet文件夹中生成TopAGNPS格式的ASCII文件：Dednm.inp，Dnmcnt.inp。由于系统运行TOPAZ模型时，程序Dednm，Raspor和Rasfor是连续运行的，所以在运行该模型前，应建立Raspor.inp和Rasfor.inp文件，并复制到当前目录中，对于不同的区域，Raspor.inp和Rasfor.inp的内容基本上没什么变化。执行程序TOPAZ及Agflow，导入TOPAZ产生的*.arc文件，这一过程将产生栅格图层（netful，subwta，bound，ntgcod，和netw），然后将“subwta”和“bound”转换成图形文件，最后subwta（grid），bound（grid），netful（grid），netw（grid），ntgcod（grid），subwta.shp和bound.shp添加在Arcview图层中。此时，在“2_AgFlow_DataSets”文件夹中生成两个文件：AnnAGNPS_Cell.dat和AnnAGNPS_Reach.dat，子流域的生成个数是受流域模拟的河网所控制的，河网越密，子流域越多，而河网的生成则应与实际的河网相吻合。(www.xing528.com)

4.1.2.3　土壤参数的确定

双阳河流域的土壤大致分为黑土、黑钙土、草甸土、沼泽土4类。AnnAGNPS需要的土壤参数包括土壤的种类和理化性质，包括：土壤的机械组成、质地，剖面深度及层次、有机质含量，氮及磷的含量、pH值等，其土壤理化性质参数见表4.1，土壤机械组成见表4.2。

表4.1　土壤理化性质参数

表4.2　土壤机械组成

各土属的水力渗透类别由RULSE（修正的通用土壤流失方程手册）查询得到。K因子反映了土壤对侵蚀的敏感性及降水所产生的径流量与径流速率的大小。影响K因子的有多方面的因素，但一般说来，质地越粗或越细的土壤有较低K值，而质地适中的反而有较高的K值。

K值估算常用的方法是Wischmeier诺漠图，根据土壤质地、土壤有机质百分含量、土壤结构、土壤透水性等几个主要因子，查土壤可蚀性因子诺漠图。但并不适合我国大多数土壤。因此，本文采用RUSLE推荐的缺少资料时采用土壤颗粒的几何平均直径计算K值的方法。公式如下：

式中　Dg——土壤颗粒的几何平均直径，mm，计算公式如下：

式中　mi——不同组i下的土壤颗粒粒径的数学平均值，mm；

fi——在mi颗粒下的质量分数。

K值的计算结果见表4.3。

表4.3　土壤的可蚀性因子K

4.1.2.4　农作物参数的确定

流域内以耕地为主，为了便于研究，本研究把流域的作物盖化为两种：水田和旱田。其中，水田的化肥平均施用量为纯氮126kg/（hm2·a），纯磷64kg/（hm2·a）。旱田纯氮35kg/（hm2·a），纯磷40kg/（hm2·a）。作物参数主要参考模型的参考手册AH-703，对作物每隔15d的作物根部比重、表面覆盖比例及降雨阻隔高度等参数参考RUSLE手册取值。

水田不同时期的氮、磷吸收率，见表4.4。因流域面积较大，很难了解到具体的作物状况，所以只能用水田和旱田的参数，旱田的吸收系数参照小麦，见表4.5。其他参数取模型默认值或参考模型提供的参考资料确定。

表4.4　水稻生长的氮磷吸收系数

表4.5　小麦生长的氮磷吸收系数

4.1.2.5　作物覆盖因子C和土壤保持因子P的确定

作物覆盖因子C反映了地面植被对地面的保护作用。除与地表覆盖有关外，还与植被灌层类型、高度等有关。在完全没有植被保护的裸露地面，C值最大为1。而保护完好的森林区C值可能为0。土壤保持因子P反映了农田中水土保持措施对侵蚀量的影响。一般无任何水土保持措施的土地P值为1，其他情况取0～1之间的数值。通过对双阳河流域的遥感图进行处理分析，并确定作物的覆盖度，结合模型自带的参数手册，最后确定了双阳河流域的作物管理C因子和土壤保持因子P，见表4.6。

表4.6　作物管理因子C和土壤保持因子P

4.1.2.6　径流曲线数CN的确定

径流曲线数（Runoff Curve Number，记为CN）是美国水土保持局提出的估算降雨径流的方法SCS中的一个指标，综合反映了土壤、土地利用、农业耕作方式、水利条件等因素对流域储水量的影响。在理论上，CN取值介于0～100之间，但在实际条件下，CN值在30～100之间变化。根据土壤渗透率和产流能力大小，该模型将土壤分成4级土壤水文类型：

A型土壤（透水）：地表径流产生能力小，在完全湿润时应有很高的下渗率，这类土壤主要是由深厚的，易大量排水的砂土和砾石组成，渗透系数为7.6～11.4cm/h。

B型土壤：渗透系数为3.8～7.6cm/h，土壤质地由中等细到中等粗，土壤中的砂粒含量少于A型，这类土壤的产流能力属于中等。

C型土壤：渗透系数为1.3～3.8cm/h，土壤质地产流能力较强，土层较浅，存在有阻碍水向下运动的土层，土壤质地常常有中到细。

D型土壤：产流能力最强，渗透系数0～1.3cm/h，主要为含有高膨胀性的黏性土壤，以及在几乎不透水物质的浅层土壤。

流域内各种土地利用的径流曲线参数（CN）的取值主要参考模型自带的资料TR55，CN的取值见表4.7。

表4.7　不同土地利用类型的CN值

4.1.2.7　曼宁糙率和地表条件常数的确定

根据研究区的土地利用情况和模型用户手册上提供的数据选取了曼宁糙率和地表条件常数，见表4.8。

表4.8　曼宁糙率和地表条件常数取值

4.1.2.8　土壤数据和土地利用数据的生成

对于每一个栅格，不仅有其高程特征，也要有与它空间位置对应的土地利用数据和土壤数据，将这些数据在GIS下，生成*.shp的图层文件，导入ArvView与AnnAGNPS的集成界面，分别与Subwta.shp这一子流域图层相交，这样可以得到每个栅格的地形数据，土壤数据以及土地利用数据的属性信息。执行程序“Extract Cell and Reach Data”，系统将产生文件“ann_cell.dbf”和“ann_reach.dbf”，并将河段数据追加到subwta.shp的属性表中，同时还产生两个Excel文件“ann_cell.csv”和“ann_reach.csv”，这两个文件可以直接导入到AnnAGNPS模型的输入数据编辑器中，成为cell和reach数据，连同其他的数据集一起作为模型的输入文件“AnnAGNPS.inp”。