楊 旭， 陳建國

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

0 引言

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

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

圖1 河網(wǎng)拓撲編碼示意圖Fig.1 Schematic diagram of river network topology

1 河網(wǎng)拓撲編碼規(guī)則及基本算法

1.1 河網(wǎng)拓撲編碼規(guī)則

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

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

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

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

1.2 河網(wǎng)拓撲編碼基本算法

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

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

圖2 矢量河網(wǎng)拓撲編碼流程圖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)編碼平臺的構(gòu)建與實現(xiàn)

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

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

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

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

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

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

MAPGIS二次開發(fā)類庫，提供了一套強有力的VC++類，它屏蔽了基于MAPGIS API之上開發(fā)MAPGIS 實用程序的許多復雜性。在MapGIS SDK支持下,根據(jù)流域匯流關(guān)系的實現(xiàn)原理,使用VC++編程實現(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é)點編號

_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é)點高程結(jié)合水往低處流的特點判定流向

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

{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);}拓撲編碼寫入對應(yīng)河網(wǎng)屬性

if(Strahler_grade[i]==1)

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

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

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

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

else if(Strahler_grade[i]==5)

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

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

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

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

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

_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é)點編號完成對具有匯流關(guān)系河段編號的查找

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

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

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

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

_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)計河網(wǎng)各級分級數(shù)量，各級河網(wǎng)長度及各級河網(wǎng)總長度

for(i=1;i

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

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

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

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

……

運行程序，主要工作界面如圖4所示。

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

3 應(yīng)用實例

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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