基于MapGIS二次开发实现矢量河网拓扑关系的自动编码与统计

杨 旭， 陈建国

(1.中国地质大学 地质过程与矿产资源国家重点实验室,武汉 710074;2．中国地质大学(武汉) 资源学院,武汉 710074)

0 引言

河网编码是水文信息管理工作的基础，编码能够直观准确地反映河网中河段间水文拓扑关系[1]，是水文管理系统研究的关键。目前采用的河网编码方法主要是通过DEM生成模拟河网，然后在模拟栅格河网上完成河网编码[2-6]。然而由于DEM本身的精度以及生成模拟河网的单流向算法[7-8]不够完善，使得实际河网与模拟河网存在较大的误差。实际河网中会出现多河道汇集，河道分叉及汇流分叉的等复杂情况[9]，单流向算法模拟河网无法体现[10]。在对河网进行编码时，李铁键等[11]提出了一种基于二叉树理论,并以二元形式表示的河网编码方法[12]。罗翔宇等[13]借助实测河网修正实测河道上部分栅格单元的高程,从修正后的DEM上提取出模拟河网后,再进行河网的编码。该方法虽提高了DEM数据的准确性，但仍无法解决单流向算法所带来的问题。

笔者提出直接在矢量河网上进行编码的思想，能够很好地解决上述问题。基于Mapgis具有强大的制图编辑功能[13]，并且MapGIS 提供的完整二次开发函数库，通过结合 MFC 类库，可以建立适合的应用软件系统。通过基于节点大小平衡二叉树的最小代价路径搜索算法[14]，对原始DEM数据中的洼地进行处理后的DEM与修正后的矢量河网进行叠置分析[15]，并将河网属性与水文属性关联，并导入Mapgis二次开发[16-17]程序平台，实现对河网中河段的拓扑关系编码并显示，同时还可以实现对河网中河段基础参数的统计，从而使得在MapGIS中实现水文分析成为可能。

图1 河网拓扑编码示意图Fig.1 Schematic diagram of river network topology

1 河网拓扑编码规则及基本算法

1.1 河网拓扑编码规则

综合Horton-Strahler分级[18-19]和Shreve-Smart分级[20]方法的优缺点，开展对河网拓扑关系的管理方法研究，并建立河网的拓扑编码体系。实际河网编码时，需要综合考虑会出现多河道汇集，河道分叉及汇流分叉的等复杂情况，因此在对河网进行编码时，遍历各水文节点所连接河段，并采用继承式编码规则(图1)，从流域各汇流源头顺流而下对各河段进行编码，使得编码确保上流河段编码始终包含下流编码信息，从而实现上下游拓扑关系运算。

在对n个河道汇流进行分级时，规定下游河段Horton-Strahler分级增加1级，Shreve-Smart分级增加n级。Horton-Strahler分级表明该n条河道汇集的下游河道为汇流主干道，Shreve-Smart分级表明该下游河道由n条支流汇集；在对河道分叉现象进行分级时，规定单河道分叉时，分叉河道以及分叉河道在下游汇流河道Horton-Strahler分级均不增加，但下游汇流河道Shreve-Smart分级增加，表明该河道是由分叉河道汇流而成；对于河道汇流分叉现象出现时，汇流分叉河道的Horton-Strahler分级与上级相同，但Shreve-Smart分级增加，表明该河道是由河流汇流分叉形成。

在完成对河网的分级后，对各河道进行编码，从而体现各河道编码的唯一性。采用Horton-Strahler分级+河道顺序码+Shreve-Smart分级的方法对各河道进行编码，由编码本身可以看出该河道的级别，以及该河道的汇流河道数(如河段编码为206(5)时，表明该河段Horton-Strahler分级为2级；06为2级河段顺序码，(5)为Shreve-Smart分级，表明流向该河道有5个上游河道汇集，即代表汇流累计量)。在进行拓扑编码时，编码规则为：该级Horton-Strahler分级+河道顺序码+Shreve-Smart分级->下一级Horton-Strahler分级+河道顺序码+Shreve-Smart分级。该编码规则实现河网中河段的直接定位,而避免低效率的遍历或搜索，编码本身具备较显著的拓扑意义，并便于识读。

水文节点对连接汇流河道起着重要的作用，水文节点与连接河段之间的对应关系及子流域与水文节点之间的对应关系构成了水文网络的拓扑关系；水文节点的上游和下游关系反映了流域的汇流关系。因此需要将水文节点类型(源头节点；汇流节点；汇出节点)，上下游连接河段编码等属性赋予水文节点。源头节点位于外链河段的外侧，内链河段以及外链河段[21]内侧节点均为汇流节点，所有河段汇流至最终出水口即为汇出节点。

1.2 河网拓扑编码基本算法

流域水文网络拓扑关系的构建是以河段以及水文节点之间的汇流关系为基础的。在构造水文拓扑结构时，确定流域尺度上的水流流向是首要环节，判断哪些水文节点是水文源点，根据“水往低处流”的规律，每个水文边线被赋予一个水流流向，水流顺着水文边线流入下游出水口，通过遍历的方式记录河段间的上下游关系，即为河段的拓扑关系。

具体算法是将预处理后的DEM数据和矢量河网进行叠置分析，叠置结果综合了DEM数据和矢量河网的要素属性，从而获取水文节点以及河网的基本属性。将河网中河段起始节点(From_Node)和终止节点(End_Node)坐标属性与水文节点坐标属性相关联，从而完成对河段两端水文节点的编号及节点所对应高程Node_Elev属性，通过河段起止水文节点的高程属性便可判断流向，并采用继承式编码规则，从流域源头节点开始顺流而下，逐河段、逐级别按河段汇流流程进行编码。河流汇流关系是通过对各条河段之间的起始节点(From_Node)和终止节点(End-Node)属性的判定来进行。判定方法为:若河段X汇入河段Y中,则河段X的终点(End_Node)必然是河段Y的起点(From_Node),从而在河段X的汇出河段编码属性中记录河段Y的唯一标识码,在河段Y的汇入河段编码属性中记录河段X的唯一标识码，并通过河段起始节点的汇流关系确定水文节点类型以及水文节点所连接的上下游河段编码。进而通过逐层递归的方法实现对河网的拓扑编码以及水文节点Node_JoLine，Node_Type(表1)等属性的完善。河段以及水文节点完成拓扑编码后，通过对河网属性的读取完成对河网基础参数的统计(图2)。在对各参数进行操作运算时，计算公式为式(1)～式(4)。

图2 矢量河网拓扑编码流程图Fig.2 Vector river network topology coding flowchart

河道长度：

(1)

河道坡度：

i= 1,2,3,…,Ns(A)

(2)

河网密度：

i= 1,2,…,Nj= 1,2,…,Np

(3)

河源密度：

j=1,2,…,Np

(4)

2 数字河网编码平台的构建与实现

2.1 数字河网编码平台的构建

MAPGIS二次开发类库是建立在MAPGIS API之上的一个类库层，用于支持基于MFC类库的面向对象的Windows程序设计，因此只需从类库派生即可实现河网拓扑编码及统计的各项功能。数字河网编码平台是以矢量河网与DEM进行叠置分析后的属性数据库为后台，VC++编写的应用程序为前台，并结合MapGIS二次开发类库实现对矢量河网的编码及统计(图3)。

MapGIS中对矢量河网(MapGIS线文件)的编码过程主要是通过代码对属性结构CATT_STRU运算的形式进行。属性结构中INFO_HEAD，FIELD_HEAD等记录属性值，因此通过_GetAtt获取矢量河网属性，并通过_LinkTbll与水文节点(MapGIS点文件)CATT_STRU属性关联，从而判断各河段流向，进而以流向为基础通过遍历河网判断各河段Horton-Strahler及Shreve-Smart分级，结合两者分级优点对各河段进行唯一编码；根据河段编码结合流向，可以判断下游河段编码，从而实现河网拓扑编码。将河网编码属性CATT_STRU通过_LinkTbll与水文节点关联，实现对水文节点类型以及链接河段的识别。根据对矢量河网CATT_STRU的遍历，实现对河网参数的统计。

表1 河网主要属性数据列

图3 河网编码平台设计框图Fig.3 Design block diagram of river network coding platform

2.2 数字河网编码平台主要功能实现

MAPGIS二次开发类库，提供了一套强有力的VC++类，它屏蔽了基于MAPGIS API之上开发MAPGIS 实用程序的许多复杂性。在MapGIS SDK支持下,根据流域汇流关系的实现原理,使用VC++编程实现。主要代码如下：

_GetLinNum(ai,&i,&n);读取工作区河网河段数

for(i=1;i

{_GetAtt(ai,LIN,i,&Stru,&ptAtt);读取工作区河网属性

_GetFieldOnNumb(ptAtt,Stru,2,(char*)(&FromNode),sizeof(long),NULL);

_GetFieldOnNumb(ptAtt,Stru,4,(char*)(&EndNode),sizeof(long),NULL);读取河网起止节点编号

_GetFieldOnNumb(ptAtt,Stru,3,(char*)(&Node1_Elevation),sizeof(long),NULL);

_GetFieldOnNumb(ptAtt,Stru,5,(char*)(&Node2_Elevation),sizeof(long),NULL);读取河段起止高程

if(Node1_Elevation

{From_Node[i]=EndNode;

End_Node[i]=FromNode;}

if(Node1_Elevation>=Node2_Elevation)

{From_Node[i]=FromNode;

End_Node[i]=EndNode;

}通过水文节点高程结合水往低处流的特点判定流向

if(Grid[i]==1) {Strahler_Code[i]=1000+n1;n1++;}

……

if(Grid[i]==5) {Strahler_Code[i]=500+n5;n5++;}基于Horton-Strahler分级进行河段编号

lstrcpy(p[i],Strahler_Code[i]);lstrcat(p[i],"(");

lstrcat(p[i],shreve);lstrcat(p[i],")");}Horton-Strahler分级和Shreve-Smart分级完成各河段编码

for(int k=1;k

{if(From_Node[k]==End_Node[i])通过对水文节点的遍历查找本级河段汇入河段

{lstrcpy(p0,s[i]);lstrcat(p0,"(");lstrcat(p0,t[i]);

lstrcat(p0,")");lstrcat(p0,"->");lstrcat(p0,s[k]);

lstrcat(p0,"(");lstrcat(p0,t[k]);lstrcat(p0,")");按照继承式编码规则编码

...... _SetFldOnNumbFrom_NodeStr(ptAtt,Stru,11,(char*)p0);

_WriteAtt(ai,LIN,i,Stru,ptAtt);}拓扑编码写入对应河网属性

if(Strahler_grade[i]==1)

{_SetFldOnNumb(ptAtt,Stru,15,"源头节点");如果河段Strahler分级为1级且节点为From_Node,则水文节点类型为源头节点

_WriteAtt(ai,LIN,i,Stru,ptAtt);

_SetFldOnNumb(ptAtt,Stru,17,"汇流节点");如果河段Strahler分级为1级且节点为End_Node,则水文节点类型为汇流节点

_WriteAtt(ai,LIN,i,Stru,ptAtt);}

else if(Strahler_grade[i]==5)

{_SetFldOnNumb(ptAtt,Stru,15,"汇流节点");如果河段Strahler分级为5级且节点为From_Node,则水文节点类型为源头节点

_SetFldOnNumb(ptAtt,Stru,17,"汇出节点");如果河段Strahler分级为5级且节点为End_Node,则水文节点类型为汇出节点

_WriteAtt(ai,LIN,i,Stru,ptAtt);

else{_SetFldOnNumb(ptAtt,Stru,15,"汇流节点");

_SetFldOnNumb(ptAtt,Stru,17,"汇流节点");如果不满足以上两个条件，则其水文节点为汇流节点

_WriteAtt(ai,LIN,i,Stru,ptAtt);}

for(j=1;j

{if(End_node[i]==From_node[j]||(i!=j&&End_node[i]==End_node[j]))根据水文节点编号完成对具有汇流关系河段编号的查找

{lstrcat(p[i],";");lstrcat(p[i],s[j]);}}统计水文节点所链接河段

_SetFldOnNumb(ptAtt,Stru,18,(char*)p[i]);

_WriteAtt(ai,LIN,i,Stru,ptAtt);将水文节点所连接汇入汇出河段写入水文节点属性

在对河流网络进行统计时，基于河网属性表在Mapgis中计算流域总面积，从而完成对河网的分级，数量，长度范围以及河网密度的统计，并以表格形式显示。主要代码如下：

_GetFieldOnNumb(ptAtt,Stru,1,(char*)(&area),sizeof(double),NULL);获取流域面积

for(i=1;i

{if(Strahler_Grade[i]==1) {Strahler_grade1++;Net_Length[k1++]=length[i];

Net_Lentotal[0]+=length[i];}

……

if(Strahler_Grade[i]==5) {Strahler_grade5++;Net_Length[k5++]=length[i];

Net_Lentotal[4]+=length[i];}

}统计河网各级分级数量，各级河网长度及各级河网总长度

for(i=1;i

Net_Lengsum+=length[i];计算河网总长度

Net_Density= Net_Lengsum *1.0/area*1000;计算河网密度

Head_Density=Strahler_grade1*1.0/area*1000000;计算河源密度

MinLength[0]=MaxLength[0]=Net_Length[0];

for (i=1; Net_Length [i]!=0;i++)

{if(MinLength [0]> Net_Length [i])

MinLength [0]= Net_Length [i];

if (MaxLength [0]< Net_Length [i])

MaxLength [0]= Net_Length [i];

}一级河网河网长度极值统计

……

运行程序，主要工作界面如图4所示。

图4 数字河网编码平台界面Fig.4 Interface of digital river coding platform

3 应用实例

通过对实验区流域的矢量河网进行编码，完成对矢量河网编码实用性的检验。由图5可知，DEM

表2 河网属性表

图5 实验区DEM与矢量河网叠加图Fig.5 DEM and vector river network overlay in experimental area

图6 实验区矢量河网编码结果图Fig.6 Results of vector river network coding in experimental area

模拟河网其位置与实际河道位置有误差，并且在主干河道处无法解决河道分叉的现象。修正后的矢量河网能准确反映河道位置及分叉的现象，因此修正后的矢量河网更为准确。

首先通过引入基于节点大小平衡二叉树的最小代价路径搜索算法对原始ASTER GDEM进行填洼处理。然后将矢量河网与处理后的DEM进行叠置分析，并将经过叠置分析处理后的矢量河网线文件加载至Mapgis二次开发数字河网编码平台，通过河网分级编码模块自动完成对河网流向的判断以及分别对各河段进行Horton-Strahler和Shreve-Smart分级编码及拓扑编码，并且程序会自动识别河段的上游和下游河段，并完成对河段属性的赋值(表2)。由表2可知，通过河网分级编码程序自动完成了对河道长度，河道坡度，河段起始水文节点的赋值。

将系统自动处理完善的河网属性表与水文节点进行关联，程序自动完成对水文节点类型的判断以及该水文节点所连接得到河段编号(表3)。由表3可知，通过河网分级编码程序自动完成了对水文节点属性的完善，主要包括各水文节点的类型，水文节点所连接的上游或下游的河道编码，从而使得水文节点与河网之间建立拓扑联系。

表3 水文节点属性表

图7 实验区矢量河网参数统计图Fig.7 Statistical diagram of vector river network parameters in experimental area

河网中各河段水系汇流编码完成后，通过河网统计模块，执行河网统计命令，系统自动完成河网中河网级别，各级别所对应的河道数，河流的长度范围，河流总长度，流域面积，河网密度，以及河源密度的的统计，并以表格形式弹出。运行结果见图7。

4 结论

笔者提出直接对矢量河网进行编码的思想，提高了河网拓扑编码的实用性。结合Horton-Strahler分级以及Shreve-Smart分级的优点为基础，针对多河道汇集，河道分叉及汇流分叉等复杂河段的编码提出了确实可行的方法。通过引入基于节点大小平衡二叉树的最小代价路径搜索算法，对原始DEM数据中的洼地进行处理，并与矢量河网进行叠置分析。该编码采用继承式编码规则，从流域源头节点开始顺流而下，逐河段、逐级别按河段汇流流程进行编码，且上游子流域编码始终包括其所流经的下游子流域编码。在此基础上，设计矢量河网编码算法，通过MapGIS SDK二次开发平台编写矢量河网自动编码系统，通过对实验区的矢量河网进行编码测试，从而验证了算法的可行性。