基于MapGIS二次開發(fā)實(shí)現(xiàn)矢量河網(wǎng)拓?fù)潢P(guān)系的自動(dòng)編碼與統(tǒng)計(jì)

楊 旭， 陳建國

(1.中國地質(zhì)大學(xué) 地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室,武漢 710074;2．中國地質(zhì)大學(xué)(武漢) 資源學(xué)院,武漢 710074)

0 引言

河網(wǎng)編碼是水文信息管理工作的基礎(chǔ)，編碼能夠直觀準(zhǔn)確地反映河網(wǎng)中河段間水文拓?fù)潢P(guān)系[1]，是水文管理系統(tǒng)研究的關(guān)鍵。目前采用的河網(wǎng)編碼方法主要是通過DEM生成模擬河網(wǎng)，然后在模擬柵格河網(wǎng)上完成河網(wǎng)編碼[2-6]。然而由于DEM本身的精度以及生成模擬河網(wǎng)的單流向算法[7-8]不夠完善，使得實(shí)際河網(wǎng)與模擬河網(wǎng)存在較大的誤差。實(shí)際河網(wǎng)中會(huì)出現(xiàn)多河道匯集，河道分叉及匯流分叉的等復(fù)雜情況[9]，單流向算法模擬河網(wǎng)無法體現(xiàn)[10]。在對河網(wǎng)進(jìn)行編碼時(shí)，李鐵鍵等[11]提出了一種基于二叉樹理論,并以二元形式表示的河網(wǎng)編碼方法[12]。羅翔宇等[13]借助實(shí)測河網(wǎng)修正實(shí)測河道上部分柵格單元的高程,從修正后的DEM上提取出模擬河網(wǎng)后,再進(jìn)行河網(wǎng)的編碼。該方法雖提高了DEM數(shù)據(jù)的準(zhǔn)確性，但仍無法解決單流向算法所帶來的問題。

筆者提出直接在矢量河網(wǎng)上進(jìn)行編碼的思想，能夠很好地解決上述問題?；贛apgis具有強(qiáng)大的制圖編輯功能[13]，并且MapGIS 提供的完整二次開發(fā)函數(shù)庫，通過結(jié)合 MFC 類庫，可以建立適合的應(yīng)用軟件系統(tǒng)。通過基于節(jié)點(diǎn)大小平衡二叉樹的最小代價(jià)路徑搜索算法[14]，對原始DEM數(shù)據(jù)中的洼地進(jìn)行處理后的DEM與修正后的矢量河網(wǎng)進(jìn)行疊置分析[15]，并將河網(wǎng)屬性與水文屬性關(guān)聯(lián)，并導(dǎo)入Mapgis二次開發(fā)[16-17]程序平臺(tái)，實(shí)現(xiàn)對河網(wǎng)中河段的拓?fù)潢P(guān)系編碼并顯示，同時(shí)還可以實(shí)現(xiàn)對河網(wǎng)中河段基礎(chǔ)參數(shù)的統(tǒng)計(jì)，從而使得在MapGIS中實(shí)現(xiàn)水文分析成為可能。

圖1 河網(wǎng)拓?fù)渚幋a示意圖Fig.1 Schematic diagram of river network topology

1 河網(wǎng)拓?fù)渚幋a規(guī)則及基本算法

1.1 河網(wǎng)拓?fù)渚幋a規(guī)則

綜合Horton-Strahler分級[18-19]和Shreve-Smart分級[20]方法的優(yōu)缺點(diǎn)，開展對河網(wǎng)拓?fù)潢P(guān)系的管理方法研究，并建立河網(wǎng)的拓?fù)渚幋a體系。實(shí)際河網(wǎng)編碼時(shí)，需要綜合考慮會(huì)出現(xiàn)多河道匯集，河道分叉及匯流分叉的等復(fù)雜情況，因此在對河網(wǎng)進(jìn)行編碼時(shí)，遍歷各水文節(jié)點(diǎn)所連接河段，并采用繼承式編碼規(guī)則(圖1)，從流域各匯流源頭順流而下對各河段進(jìn)行編碼，使得編碼確保上流河段編碼始終包含下流編碼信息，從而實(shí)現(xiàn)上下游拓?fù)潢P(guān)系運(yùn)算。

在對n個(gè)河道匯流進(jìn)行分級時(shí)，規(guī)定下游河段Horton-Strahler分級增加1級，Shreve-Smart分級增加n級。Horton-Strahler分級表明該n條河道匯集的下游河道為匯流主干道，Shreve-Smart分級表明該下游河道由n條支流匯集；在對河道分叉現(xiàn)象進(jìn)行分級時(shí)，規(guī)定單河道分叉時(shí)，分叉河道以及分叉河道在下游匯流河道Horton-Strahler分級均不增加，但下游匯流河道Shreve-Smart分級增加，表明該河道是由分叉河道匯流而成；對于河道匯流分叉現(xiàn)象出現(xiàn)時(shí)，匯流分叉河道的Horton-Strahler分級與上級相同，但Shreve-Smart分級增加，表明該河道是由河流匯流分叉形成。

在完成對河網(wǎng)的分級后，對各河道進(jìn)行編碼，從而體現(xiàn)各河道編碼的唯一性。采用Horton-Strahler分級+河道順序碼+Shreve-Smart分級的方法對各河道進(jìn)行編碼，由編碼本身可以看出該河道的級別，以及該河道的匯流河道數(shù)(如河段編碼為206(5)時(shí)，表明該河段Horton-Strahler分級為2級；06為2級河段順序碼，(5)為Shreve-Smart分級，表明流向該河道有5個(gè)上游河道匯集，即代表匯流累計(jì)量)。在進(jìn)行拓?fù)渚幋a時(shí)，編碼規(guī)則為：該級Horton-Strahler分級+河道順序碼+Shreve-Smart分級->下一級Horton-Strahler分級+河道順序碼+Shreve-Smart分級。該編碼規(guī)則實(shí)現(xiàn)河網(wǎng)中河段的直接定位,而避免低效率的遍歷或搜索，編碼本身具備較顯著的拓?fù)湟饬x，并便于識(shí)讀。

水文節(jié)點(diǎn)對連接匯流河道起著重要的作用，水文節(jié)點(diǎn)與連接河段之間的對應(yīng)關(guān)系及子流域與水文節(jié)點(diǎn)之間的對應(yīng)關(guān)系構(gòu)成了水文網(wǎng)絡(luò)的拓?fù)潢P(guān)系；水文節(jié)點(diǎn)的上游和下游關(guān)系反映了流域的匯流關(guān)系。因此需要將水文節(jié)點(diǎn)類型(源頭節(jié)點(diǎn)；匯流節(jié)點(diǎn)；匯出節(jié)點(diǎn))，上下游連接河段編碼等屬性賦予水文節(jié)點(diǎn)。源頭節(jié)點(diǎn)位于外鏈河段的外側(cè)，內(nèi)鏈河段以及外鏈河段[21]內(nèi)側(cè)節(jié)點(diǎn)均為匯流節(jié)點(diǎn)，所有河段匯流至最終出水口即為匯出節(jié)點(diǎn)。

1.2 河網(wǎng)拓?fù)渚幋a基本算法

流域水文網(wǎng)絡(luò)拓?fù)潢P(guān)系的構(gòu)建是以河段以及水文節(jié)點(diǎn)之間的匯流關(guān)系為基礎(chǔ)的。在構(gòu)造水文拓?fù)浣Y(jié)構(gòu)時(shí)，確定流域尺度上的水流流向是首要環(huán)節(jié)，判斷哪些水文節(jié)點(diǎn)是水文源點(diǎn)，根據(jù)“水往低處流”的規(guī)律，每個(gè)水文邊線被賦予一個(gè)水流流向，水流順著水文邊線流入下游出水口，通過遍歷的方式記錄河段間的上下游關(guān)系，即為河段的拓?fù)潢P(guān)系。

具體算法是將預(yù)處理后的DEM數(shù)據(jù)和矢量河網(wǎng)進(jìn)行疊置分析，疊置結(jié)果綜合了DEM數(shù)據(jù)和矢量河網(wǎng)的要素屬性，從而獲取水文節(jié)點(diǎn)以及河網(wǎng)的基本屬性。將河網(wǎng)中河段起始節(jié)點(diǎn)(From_Node)和終止節(jié)點(diǎn)(End_Node)坐標(biāo)屬性與水文節(jié)點(diǎn)坐標(biāo)屬性相關(guān)聯(lián)，從而完成對河段兩端水文節(jié)點(diǎn)的編號及節(jié)點(diǎn)所對應(yīng)高程N(yùn)ode_Elev屬性，通過河段起止水文節(jié)點(diǎn)的高程屬性便可判斷流向，并采用繼承式編碼規(guī)則，從流域源頭節(jié)點(diǎn)開始順流而下，逐河段、逐級別按河段匯流流程進(jìn)行編碼。河流匯流關(guān)系是通過對各條河段之間的起始節(jié)點(diǎn)(From_Node)和終止節(jié)點(diǎn)(End-Node)屬性的判定來進(jìn)行。判定方法為:若河段X匯入河段Y中,則河段X的終點(diǎn)(End_Node)必然是河段Y的起點(diǎn)(From_Node),從而在河段X的匯出河段編碼屬性中記錄河段Y的唯一標(biāo)識(shí)碼,在河段Y的匯入河段編碼屬性中記錄河段X的唯一標(biāo)識(shí)碼，并通過河段起始節(jié)點(diǎn)的匯流關(guān)系確定水文節(jié)點(diǎn)類型以及水文節(jié)點(diǎn)所連接的上下游河段編碼。進(jìn)而通過逐層遞歸的方法實(shí)現(xiàn)對河網(wǎng)的拓?fù)渚幋a以及水文節(jié)點(diǎn)Node_JoLine，Node_Type(表1)等屬性的完善。河段以及水文節(jié)點(diǎn)完成拓?fù)渚幋a后，通過對河網(wǎng)屬性的讀取完成對河網(wǎng)基礎(chǔ)參數(shù)的統(tǒng)計(jì)(圖2)。在對各參數(shù)進(jìn)行操作運(yùn)算時(shí)，計(jì)算公式為式(1)～式(4)。

圖2 矢量河網(wǎng)拓?fù)渚幋a流程圖Fig.2 Vector river network topology coding flowchart

河道長度：

(1)

河道坡度：

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

(2)

河網(wǎng)密度：

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

(3)

河源密度：

j=1,2,…,Np

(4)

2 數(shù)字河網(wǎng)編碼平臺(tái)的構(gòu)建與實(shí)現(xiàn)

2.1 數(shù)字河網(wǎng)編碼平臺(tái)的構(gòu)建

MAPGIS二次開發(fā)類庫是建立在MAPGIS API之上的一個(gè)類庫層，用于支持基于MFC類庫的面向?qū)ο蟮腤indows程序設(shè)計(jì)，因此只需從類庫派生即可實(shí)現(xiàn)河網(wǎng)拓?fù)渚幋a及統(tǒng)計(jì)的各項(xiàng)功能。數(shù)字河網(wǎng)編碼平臺(tái)是以矢量河網(wǎng)與DEM進(jìn)行疊置分析后的屬性數(shù)據(jù)庫為后臺(tái)，VC++編寫的應(yīng)用程序?yàn)榍芭_(tái)，并結(jié)合MapGIS二次開發(fā)類庫實(shí)現(xiàn)對矢量河網(wǎng)的編碼及統(tǒng)計(jì)(圖3)。

MapGIS中對矢量河網(wǎng)(MapGIS線文件)的編碼過程主要是通過代碼對屬性結(jié)構(gòu)CATT_STRU運(yùn)算的形式進(jìn)行。屬性結(jié)構(gòu)中INFO_HEAD，F(xiàn)IELD_HEAD等記錄屬性值，因此通過_GetAtt獲取矢量河網(wǎng)屬性，并通過_LinkTbll與水文節(jié)點(diǎn)(MapGIS點(diǎn)文件)CATT_STRU屬性關(guān)聯(lián)，從而判斷各河段流向，進(jìn)而以流向?yàn)榛A(chǔ)通過遍歷河網(wǎng)判斷各河段Horton-Strahler及Shreve-Smart分級，結(jié)合兩者分級優(yōu)點(diǎn)對各河段進(jìn)行唯一編碼；根據(jù)河段編碼結(jié)合流向，可以判斷下游河段編碼，從而實(shí)現(xiàn)河網(wǎng)拓?fù)渚幋a。將河網(wǎng)編碼屬性CATT_STRU通過_LinkTbll與水文節(jié)點(diǎn)關(guān)聯(lián)，實(shí)現(xiàn)對水文節(jié)點(diǎn)類型以及鏈接河段的識(shí)別。根據(jù)對矢量河網(wǎng)CATT_STRU的遍歷，實(shí)現(xiàn)對河網(wǎng)參數(shù)的統(tǒng)計(jì)。

表1 河網(wǎng)主要屬性數(shù)據(jù)列

圖3 河網(wǎng)編碼平臺(tái)設(shè)計(jì)框圖Fig.3 Design block diagram of river network coding platform

2.2 數(shù)字河網(wǎng)編碼平臺(tái)主要功能實(shí)現(xiàn)

MAPGIS二次開發(fā)類庫，提供了一套強(qiáng)有力的VC++類，它屏蔽了基于MAPGIS API之上開發(fā)MAPGIS 實(shí)用程序的許多復(fù)雜性。在MapGIS SDK支持下,根據(jù)流域匯流關(guān)系的實(shí)現(xiàn)原理,使用VC++編程實(shí)現(xiàn)。主要代碼如下：

_GetLinNum(ai,&i,&n);讀取工作區(qū)河網(wǎng)河段數(shù)

for(i=1;i

{_GetAtt(ai,LIN,i,&Stru,&ptAtt);讀取工作區(qū)河網(wǎng)屬性

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

_GetFieldOnNumb(ptAtt,Stru,4,(char*)(&EndNode),sizeof(long),NULL);讀取河網(wǎng)起止節(jié)點(diǎn)編號

_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;

}通過水文節(jié)點(diǎn)高程結(jié)合水往低處流的特點(diǎn)判定流向

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

……

if(Grid[i]==5) {Strahler_Code[i]=500+n5;n5++;}基于Horton-Strahler分級進(jìn)行河段編號

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])通過對水文節(jié)點(diǎn)的遍歷查找本級河段匯入河段

{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,")");按照繼承式編碼規(guī)則編碼

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

_WriteAtt(ai,LIN,i,Stru,ptAtt);}拓?fù)渚幋a寫入對應(yīng)河網(wǎng)屬性

if(Strahler_grade[i]==1)

{_SetFldOnNumb(ptAtt,Stru,15,"源頭節(jié)點(diǎn)");如果河段Strahler分級為1級且節(jié)點(diǎn)為From_Node,則水文節(jié)點(diǎn)類型為源頭節(jié)點(diǎn)

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

_SetFldOnNumb(ptAtt,Stru,17,"匯流節(jié)點(diǎn)");如果河段Strahler分級為1級且節(jié)點(diǎn)為End_Node,則水文節(jié)點(diǎn)類型為匯流節(jié)點(diǎn)

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

else if(Strahler_grade[i]==5)

{_SetFldOnNumb(ptAtt,Stru,15,"匯流節(jié)點(diǎn)");如果河段Strahler分級為5級且節(jié)點(diǎn)為From_Node,則水文節(jié)點(diǎn)類型為源頭節(jié)點(diǎn)

_SetFldOnNumb(ptAtt,Stru,17,"匯出節(jié)點(diǎn)");如果河段Strahler分級為5級且節(jié)點(diǎn)為End_Node,則水文節(jié)點(diǎn)類型為匯出節(jié)點(diǎn)

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

else{_SetFldOnNumb(ptAtt,Stru,15,"匯流節(jié)點(diǎn)");

_SetFldOnNumb(ptAtt,Stru,17,"匯流節(jié)點(diǎn)");如果不滿足以上兩個(gè)條件，則其水文節(jié)點(diǎn)為匯流節(jié)點(diǎn)

_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]))根據(jù)水文節(jié)點(diǎn)編號完成對具有匯流關(guān)系河段編號的查找

{lstrcat(p[i],";");lstrcat(p[i],s[j]);}}統(tǒng)計(jì)水文節(jié)點(diǎn)所鏈接河段

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

_WriteAtt(ai,LIN,i,Stru,ptAtt);將水文節(jié)點(diǎn)所連接匯入?yún)R出河段寫入水文節(jié)點(diǎn)屬性

在對河流網(wǎng)絡(luò)進(jìn)行統(tǒng)計(jì)時(shí)，基于河網(wǎng)屬性表在Mapgis中計(jì)算流域總面積，從而完成對河網(wǎng)的分級，數(shù)量，長度范圍以及河網(wǎng)密度的統(tǒng)計(jì)，并以表格形式顯示。主要代碼如下：

_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];}

}統(tǒng)計(jì)河網(wǎng)各級分級數(shù)量，各級河網(wǎng)長度及各級河網(wǎng)總長度

for(i=1;i

Net_Lengsum+=length[i];計(jì)算河網(wǎng)總長度

Net_Density= Net_Lengsum *1.0/area*1000;計(jì)算河網(wǎng)密度

Head_Density=Strahler_grade1*1.0/area*1000000;計(jì)算河源密度

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];

}一級河網(wǎng)河網(wǎng)長度極值統(tǒng)計(jì)

……

運(yùn)行程序，主要工作界面如圖4所示。

圖4 數(shù)字河網(wǎng)編碼平臺(tái)界面Fig.4 Interface of digital river coding platform

3 應(yīng)用實(shí)例

通過對實(shí)驗(yàn)區(qū)流域的矢量河網(wǎng)進(jìn)行編碼，完成對矢量河網(wǎng)編碼實(shí)用性的檢驗(yàn)。由圖5可知，DEM

表2 河網(wǎng)屬性表

圖5 實(shí)驗(yàn)區(qū)DEM與矢量河網(wǎng)疊加圖Fig.5 DEM and vector river network overlay in experimental area

圖6 實(shí)驗(yàn)區(qū)矢量河網(wǎng)編碼結(jié)果圖Fig.6 Results of vector river network coding in experimental area

模擬河網(wǎng)其位置與實(shí)際河道位置有誤差，并且在主干河道處無法解決河道分叉的現(xiàn)象。修正后的矢量河網(wǎng)能準(zhǔn)確反映河道位置及分叉的現(xiàn)象，因此修正后的矢量河網(wǎng)更為準(zhǔn)確。

首先通過引入基于節(jié)點(diǎn)大小平衡二叉樹的最小代價(jià)路徑搜索算法對原始ASTER GDEM進(jìn)行填洼處理。然后將矢量河網(wǎng)與處理后的DEM進(jìn)行疊置分析，并將經(jīng)過疊置分析處理后的矢量河網(wǎng)線文件加載至Mapgis二次開發(fā)數(shù)字河網(wǎng)編碼平臺(tái)，通過河網(wǎng)分級編碼模塊自動(dòng)完成對河網(wǎng)流向的判斷以及分別對各河段進(jìn)行Horton-Strahler和Shreve-Smart分級編碼及拓?fù)渚幋a，并且程序會(huì)自動(dòng)識(shí)別河段的上游和下游河段，并完成對河段屬性的賦值(表2)。由表2可知，通過河網(wǎng)分級編碼程序自動(dòng)完成了對河道長度，河道坡度，河段起始水文節(jié)點(diǎn)的賦值。

將系統(tǒng)自動(dòng)處理完善的河網(wǎng)屬性表與水文節(jié)點(diǎn)進(jìn)行關(guān)聯(lián)，程序自動(dòng)完成對水文節(jié)點(diǎn)類型的判斷以及該水文節(jié)點(diǎn)所連接得到河段編號(表3)。由表3可知，通過河網(wǎng)分級編碼程序自動(dòng)完成了對水文節(jié)點(diǎn)屬性的完善，主要包括各水文節(jié)點(diǎn)的類型，水文節(jié)點(diǎn)所連接的上游或下游的河道編碼，從而使得水文節(jié)點(diǎn)與河網(wǎng)之間建立拓?fù)渎?lián)系。

表3 水文節(jié)點(diǎn)屬性表

圖7 實(shí)驗(yàn)區(qū)矢量河網(wǎng)參數(shù)統(tǒng)計(jì)圖Fig.7 Statistical diagram of vector river network parameters in experimental area

河網(wǎng)中各河段水系匯流編碼完成后，通過河網(wǎng)統(tǒng)計(jì)模塊，執(zhí)行河網(wǎng)統(tǒng)計(jì)命令，系統(tǒng)自動(dòng)完成河網(wǎng)中河網(wǎng)級別，各級別所對應(yīng)的河道數(shù)，河流的長度范圍，河流總長度，流域面積，河網(wǎng)密度，以及河源密度的的統(tǒng)計(jì)，并以表格形式彈出。運(yùn)行結(jié)果見圖7。

4 結(jié)論

筆者提出直接對矢量河網(wǎng)進(jìn)行編碼的思想，提高了河網(wǎng)拓?fù)渚幋a的實(shí)用性。結(jié)合Horton-Strahler分級以及Shreve-Smart分級的優(yōu)點(diǎn)為基礎(chǔ)，針對多河道匯集，河道分叉及匯流分叉等復(fù)雜河段的編碼提出了確實(shí)可行的方法。通過引入基于節(jié)點(diǎn)大小平衡二叉樹的最小代價(jià)路徑搜索算法，對原始DEM數(shù)據(jù)中的洼地進(jìn)行處理，并與矢量河網(wǎng)進(jìn)行疊置分析。該編碼采用繼承式編碼規(guī)則，從流域源頭節(jié)點(diǎn)開始順流而下，逐河段、逐級別按河段匯流流程進(jìn)行編碼，且上游子流域編碼始終包括其所流經(jīng)的下游子流域編碼。在此基礎(chǔ)上，設(shè)計(jì)矢量河網(wǎng)編碼算法，通過MapGIS SDK二次開發(fā)平臺(tái)編寫矢量河網(wǎng)自動(dòng)編碼系統(tǒng)，通過對實(shí)驗(yàn)區(qū)的矢量河網(wǎng)進(jìn)行編碼測試，從而驗(yàn)證了算法的可行性。