基于公共DEM的平原湖泊流域特征提取研究及應(yīng)用

周 星，仲志余

(1. 長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究院，湖北 武漢 430010；2. 武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

流域特征(面積、河網(wǎng)結(jié)構(gòu)、主河流長(zhǎng)度和主河流比降等)在水利規(guī)劃設(shè)計(jì)、水文分析計(jì)算、水文模擬以及流域管理等方面有重要的作用。隨著3S技術(shù)的發(fā)展，基于數(shù)字高程模型(Digital Elevation model，DEM)提取流域水文特征已經(jīng)成為主流。國(guó)內(nèi)外學(xué)者提出了眾多基于DEM的河網(wǎng)提取方法，其中O Callaghan等最早提出的基于地形坡度計(jì)算流向(Deterministic-8 node算法，D8)和匯流累積量而確定河網(wǎng)的方法應(yīng)用最為廣泛[1-3]，美國(guó)ESRI(Environmental Systems Research Institute)公司的Arcgis軟件中的水文模塊、美國(guó)陸軍工程公司研制的GRASS(Geographic Resources Analysis Support System)軟件、River Tools軟件等都集成了該算法，從而實(shí)現(xiàn)了河網(wǎng)提取的自動(dòng)化及可視化。目前學(xué)者們對(duì)于流域特征提取研究多是研究DEM分辨率、河網(wǎng)閾值取值對(duì)河網(wǎng)提取結(jié)果的影響[4-8]，然而這些研究還缺乏對(duì)誤差的量化分析，同時(shí)對(duì)于水利規(guī)劃設(shè)計(jì)中重要的流域特征參數(shù)主河流長(zhǎng)度和主河流比降提取還缺乏相應(yīng)的研究。

在實(shí)際工作中，流域缺乏測(cè)繪數(shù)據(jù)或者只有部分測(cè)繪數(shù)據(jù)，對(duì)于這些流域的特征提取需要尋求解決的方法。作為全球尺度的公開(kāi)的基礎(chǔ)地形數(shù)據(jù)，SRTM(Shuttle Radar Topography Mission, 即航天飛機(jī)雷達(dá)地形測(cè)繪任務(wù)) DEM和ASTER (The Advanced Spaceborne Thermal Emission and Reflection Radiometer，即星載熱發(fā)射和反射輻射儀) GDEM的精度得到了廣泛的認(rèn)可[9-10]，并在許多大中型流域河網(wǎng)提取中得到成功的應(yīng)用[13-15]。此外隨著遙測(cè)技術(shù)的發(fā)展，越來(lái)越多的衛(wèi)星影像地圖公開(kāi)，諸如Google衛(wèi)星地圖、微軟Bing地圖、中國(guó)的天地圖等分辨率都達(dá)到了米級(jí)甚至分米級(jí)，通過(guò)這些影像可以輔助流域特征的提取。鑒于此，本研究在高清衛(wèi)星影像的輔助下，提出了基于公共DEM數(shù)據(jù)的流域特征提取的方法并將之應(yīng)用于典型的平原湖泊流域——花馬湖流域。

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域概況

花馬湖流域(圖1)位于鄂州市東部，東經(jīng)114°52′53.73″～115°4′30.68″，北緯30°12′8.63″～30°23′47.99″之間，流域集水面積291 km2，境內(nèi)地勢(shì)北高南低，北面東面面臨長(zhǎng)江，出口位于長(zhǎng)江?；R湖流域?yàn)榈湫偷钠皆戳饔颍饔騼?nèi)湖泊眾多，主要的湖泊有走馬湖、花馬湖和花家湖。流域內(nèi)地勢(shì)平坦，氣候?yàn)榈湫偷膩啛釒Т箨懶约撅L(fēng)氣候，根據(jù)流域鄰近的鄂州氣象站資料統(tǒng)計(jì)，多年平均氣溫17℃，多年平均降水量1 301.2 mm。此外流域開(kāi)發(fā)程度較高，尤其是東部及南部城市化明顯。

1.2 數(shù)據(jù)

采用SRTM DEM數(shù)據(jù)和ASTER GDEM數(shù)據(jù)。其中SRTM采用經(jīng)過(guò)多次優(yōu)化校正后的4.1版本的90 m分辨率數(shù)據(jù)；ASTER數(shù)據(jù)采用第二代版本的數(shù)據(jù)，空間分辨率為30 m。兩種數(shù)據(jù)均來(lái)源于中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心國(guó)際科學(xué)數(shù)據(jù)鏡像網(wǎng)站(http://www.gscloud.cn)。兩種數(shù)據(jù)的具體參數(shù)見(jiàn)表1。

2 研究方法

本研究采用Arcgis軟件作為流域河網(wǎng)提取平臺(tái)，具體采用Arc Hydro Tools作為主要工具欄。Arc Hydro Tools是基于ArcGIS 和Arc Hydro數(shù)據(jù)模型開(kāi)發(fā)的1套用于支持地表水資源應(yīng)用研究的工具集，相比于Arcgis自帶的Hydrology工具集功能更加全面，最新版本可以在ESRI 公司官方網(wǎng)站http://www.esri.com 免費(fèi)下載。對(duì)于完整的流域特征提取可分為矢量河網(wǎng)解譯、高程預(yù)處理以及流域特征提取計(jì)算3部分，每個(gè)部分通過(guò)單個(gè)或多個(gè)命令實(shí)現(xiàn)，并且可通過(guò)建?；蛘呔幊虒?shí)現(xiàn)集成，其總體流程見(jiàn)圖2。

表1 SRTM 、ASTER高程數(shù)據(jù)參數(shù)

2.1 矢量河網(wǎng)的獲取

采用Google無(wú)偏移衛(wèi)星影像識(shí)別流域內(nèi)主要河流及湖泊并矢量化，建立河網(wǎng)矢量化圖層River(圖1)。Google影像共分20個(gè)級(jí)別，級(jí)別越大影像空間分辨率越高，20級(jí)空間分辨率達(dá)到0.27 m，相當(dāng)于1∶800比例尺。20級(jí)影像下若沒(méi)有云層遮擋，能夠清晰識(shí)別大部分河流。對(duì)于湖泊將其中心線作為河流。在進(jìn)行矢量河網(wǎng)解譯前，需要確保衛(wèi)星影像與DEM的空間位置匹配。GoogleEarth原始影像采用坐標(biāo)系統(tǒng)為WGS84系統(tǒng)，在中國(guó)區(qū)域有系統(tǒng)性的偏差，在實(shí)際應(yīng)用中可在研究區(qū)域中選擇一些特征點(diǎn)對(duì)齊用GPS儀等準(zhǔn)確測(cè)量，然后采用Georeferencing工具進(jìn)行校正，消除系統(tǒng)性偏差。在本文研究中，選用的衛(wèi)星影像為商業(yè)軟件提供的Google衛(wèi)星影像，其在原始的影像基礎(chǔ)上對(duì)中國(guó)區(qū)域內(nèi)做了偏移校正，因此本研究中不再考慮衛(wèi)星影像的偏差。

2.2 矢量河網(wǎng)的獲取

DEM的預(yù)處理包括河網(wǎng)校正和填洼。諸多研究表明提取河網(wǎng)時(shí)易受洼地、人類(lèi)活動(dòng)及DEM測(cè)量精度的影響而產(chǎn)生不連續(xù)的河網(wǎng)、平行河網(wǎng)以及位置錯(cuò)誤的河網(wǎng)等，通過(guò)增加有效信息可減弱這些不利影響[12-13]。目前提高DEM提取河網(wǎng)精度的方法主要有洼地去除(sink)、河網(wǎng)校正(AGREE)和河流燒錄(Stream burning)等方法[16-17]。本研究采用河網(wǎng)校正和填洼方法DEM進(jìn)行預(yù)處理，預(yù)處理的主要步驟如下。

a) DEM由于精度不夠等原因不能夠正確反映河流的真實(shí)位置，尤其是在平原湖泊流域，地形起伏小，河流位置的水面高程與周?chē)嗖畈淮螅偌由螪EM分辨率無(wú)法識(shí)別堤防等小尺度建筑，有些河流的位置無(wú)法正確識(shí)別(圖3a)。河網(wǎng)校正使用DEM Reconditioning命令(AGREE算法)或者Stream burning來(lái)實(shí)現(xiàn)。AGREE算法的基本原理是通過(guò)人工輸入的矢量河網(wǎng)，準(zhǔn)確定位河流的位置，給予河流位置處DEM一個(gè)很深的強(qiáng)迫，并在河流位置一定的距離范圍內(nèi)線性插值，使得周?chē)乃髂軌蛄鞯胶恿魈?圖3b)；Stream burning的原理與AGREE類(lèi)似，不同的是其保持河流位置的DEM值不變，而增加其他位置的DEM。

b) 由于DEM數(shù)據(jù)制作的時(shí)候采樣效果和高程值設(shè)置為整數(shù)以及受測(cè)量時(shí)候噪聲的影響，當(dāng)DEM分辨率較大時(shí)往往會(huì)導(dǎo)致洼地的形成。顧名思義洼地就是高程低于周?chē)臇鸥?，所以水流都?huì)匯入其中，導(dǎo)致徑流最終斷流，需要進(jìn)行處理。本研究中，采用Fill Sinks命令實(shí)現(xiàn)。其基本原理見(jiàn)圖4，當(dāng)檢測(cè)到洼地時(shí)，將洼地處的高程抬高到與周?chē)裨懈叱套畹偷南裨恢?，并不斷地循環(huán)此步驟直至不再出現(xiàn)新的洼地。

2.3 水文分析

水文分析包括流向分析、匯流累積量計(jì)算、河流柵格定義、河流拓?fù)潢P(guān)系生成以及河流邊界生成等，主要步驟如下。

a) 流向分析的目的是明確每個(gè)單元的水流聚集點(diǎn)水流方向。使用Flow Direction命令(D8算法)，D8算法核心思想是通過(guò)計(jì)算柵格與周?chē)?個(gè)柵格的坡度，選取最大坡度的方向作為本柵格的流向。

b) 采用Flow Accumulation命令計(jì)算每個(gè)柵格的集水面積，該命令通過(guò)流向柵格圖搜索水流路徑，采用遞歸算法從流域出口柵格開(kāi)始遞歸搜索計(jì)算出每一柵格單元的上游匯水面積。

c) 得到匯水面積圖后，使用Stream Definition命令確定河流柵格圖，該命令的核心思想是設(shè)定1個(gè)最小匯水面積閾值(本研究設(shè)定為1 km2)，高于這個(gè)閾值的柵格將定義為河流。

至此，流域的河流已自動(dòng)提取出來(lái)，可進(jìn)一步做處理，如采用Stream Segmentation命令，建立柵格河流上下游拓?fù)潢P(guān)系；采用Catchment Grid Delineation命令，得到每條河段的邊界(即河段控制的子流域)；采用使用Drainage Line Processing和Catchment Polygon Processing命令將所得的柵格數(shù)據(jù)轉(zhuǎn)化為矢量數(shù)據(jù)以便進(jìn)行后續(xù)的流域邊界、河流長(zhǎng)度和坡度的提取。

2.4 流域特征提取計(jì)算

流域面積、河流長(zhǎng)度及比降是非常重要的流域參數(shù)，可根據(jù)通過(guò)地形數(shù)據(jù)提取。流域面積有著統(tǒng)一的定義即為出口斷面以上的集水區(qū)域。然而河流長(zhǎng)度的內(nèi)涵在不同領(lǐng)域時(shí)有所不同，以往的基于DEM流域特征提取的研究中很少關(guān)注這種內(nèi)涵的不同，幾乎所有的研究中河流內(nèi)涵都是以出現(xiàn)明顯水流(以匯水面積閾值來(lái)反應(yīng))為標(biāo)志的，這與傳統(tǒng)的認(rèn)知相符，顯然自動(dòng)提取的河流長(zhǎng)度是隨著閾值變化而變化，故較多研究集中在怎樣合理確定匯水面積閾值。然而在流域設(shè)計(jì)洪水計(jì)算中，河流長(zhǎng)度的內(nèi)涵為出口斷面沿主河流至分水嶺的長(zhǎng)度，這個(gè)長(zhǎng)度是確定的。傳統(tǒng)的研究中，在進(jìn)行水文分析操作后就結(jié)束了，為了研究河流的內(nèi)涵不同對(duì)流域參數(shù)的影響，需要對(duì)流域設(shè)計(jì)洪水領(lǐng)域中的河流長(zhǎng)度以及河流比降進(jìn)行提取，其中河流比降定義為主河流各高程轉(zhuǎn)折點(diǎn)分段比降的加權(quán)平均值，計(jì)算公式見(jiàn)式(1)。

(1)

式中H0——河流末節(jié)點(diǎn)高程；Hi——河段節(jié)點(diǎn)高程；li——河段長(zhǎng)度；L——總河流長(zhǎng)度。

流域設(shè)計(jì)洪水領(lǐng)域中的河流長(zhǎng)度及比降提取步驟如下。

a) 合理地進(jìn)行分區(qū)(子流域)，并確定子流域的出口即子流域匯入到主干的匯水口，使用Batch Point Generation命令建立子流域出口的文件。

b) 使用Batch Subwatershed Delineation命令，確定子流域面積，其核心思想是以分區(qū)匯水口為起點(diǎn)，通過(guò)流向柵格圖搜索水流路徑得到匯水口的上游面積。

c) 使用Longest Flow Path for Subwatersheds命令確定每個(gè)子流域內(nèi)從分水嶺至出口的最長(zhǎng)路徑即為子流域內(nèi)的主河流，其長(zhǎng)度即為主河流長(zhǎng)度，其思想是通過(guò)子流域區(qū)域內(nèi)流向柵格圖搜索最長(zhǎng)的匯流路徑。

d) 利用Stack Profile命令得到主河流的地形縱剖線，其思想是以主河流起點(diǎn)開(kāi)始，提取主河流特征點(diǎn)處的原始高程值及距離起點(diǎn)的距離，最后可按照公式(1)計(jì)算可得主河流的平均比降。

3 結(jié)果分析

3.1 河網(wǎng)校正對(duì)流域特征提取的影響

以SRTM為數(shù)據(jù)源，在是否進(jìn)行河網(wǎng)校正步驟下對(duì)花馬湖流域進(jìn)行流域特征提取，并分析河網(wǎng)校正對(duì)流域特征提取的改進(jìn)效果，結(jié)果見(jiàn)圖5。對(duì)比可知，不進(jìn)行河網(wǎng)校正提取的河網(wǎng)在河源及高程較高的區(qū)域與真實(shí)河網(wǎng)吻合較好，在地勢(shì)低的地區(qū)(平原地區(qū))誤差較大，出現(xiàn)較多的平行河流。輔助真實(shí)河流信息，對(duì)原始DEM進(jìn)行河網(wǎng)校正后，極大地提高了平原地區(qū)的河網(wǎng)精度，改善了平行河流的問(wèn)題。此外，采用河網(wǎng)校正后，流域出口的位置也校正到了正確的位置。

為了量化河網(wǎng)的偏差程度，提取矢量化的真實(shí)河網(wǎng)的節(jié)點(diǎn)，計(jì)算其離提取河網(wǎng)的距離(即偏移距離)，并計(jì)算偏移距離的統(tǒng)計(jì)特征值(表2)。采用AGREE校正后，河網(wǎng)的平均偏移距離從95 m減小到54 m，減小了43%；偏移距離的標(biāo)準(zhǔn)差也從130 m減小到97 m，減小了23%；采用Stream burning校正后，河網(wǎng)的平均偏移距離減小到51 m，減小了46%；偏移距離的標(biāo)準(zhǔn)差減小到92 m，減小了29%。

Stream burning校正方法提取的精度略高于AGREE方法，但是沒(méi)有顯著的區(qū)別；AGREE校正后提取的河網(wǎng)細(xì)節(jié)更加豐富，有更多的支流?？傮w而言，AGREE校正和Stream burning校正方法均明顯提升了河網(wǎng)提取精度，尤其適用于平原地區(qū)的河網(wǎng)提取，在實(shí)際應(yīng)用中可根據(jù)需要選擇應(yīng)用。

表2 校正算法對(duì)河網(wǎng)提取精度的影響

河網(wǎng)校正后，流域面積及河網(wǎng)密度有一定變化(表 3)。校正后河網(wǎng)密度有所下降，這是因?yàn)樾Ｕ舜蟛糠值钠叫泻恿鳌２恍Ｕ崛〉牧饔蛎娣e為289.51 km2，AGREE校正提取的流域面積為289.59 km2，Stream burning提取的面積為292 km2，三者的面積相差很小。整體上三者的流域邊界較為吻合，但在西北及東北局部有著較為明顯的區(qū)別，這是由于人類(lèi)活動(dòng)造成的。東部及北部河流經(jīng)歷過(guò)整治，沿江修建有長(zhǎng)江干堤，部分河流在流域邊界外，可對(duì)提取結(jié)果進(jìn)行人工修正，如東部邊界以長(zhǎng)江干堤為界限。

表3 不同數(shù)據(jù)源對(duì)流域特征提取的影響

3.2 不同數(shù)據(jù)源對(duì)流域特征提取的影響

分別選取SRTM和ASTER數(shù)據(jù)對(duì)花馬湖流域進(jìn)行AGREE校正算法的流域特征提取，分析不同高程數(shù)據(jù)源對(duì)于流域特征提取的影響(圖6、表3)。與SRTM數(shù)據(jù)類(lèi)似，ASTER數(shù)據(jù)提取的河網(wǎng)提取精度從山區(qū)至平原區(qū)逐漸降低。隨著DEM的水平分辨率增加，對(duì)地形起伏刻畫(huà)得更加精細(xì)，因此ASTER數(shù)據(jù)提取的河流長(zhǎng)度以及河網(wǎng)密度都略大于SRTM數(shù)據(jù)的結(jié)果，這與一些學(xué)者的研究結(jié)果相符[5,7]， 并且流域邊界更加光滑。但是ASTER數(shù)據(jù)提取的流域面積顯著大于SRTM數(shù)據(jù)提取的面積，也顯著大于水利部門(mén)審查通過(guò)的291 km2，甚至其流域的出口也偏移巨大。究其原因， 雖然GDEM的水平分辨率(30 m)高于SRTM(90 m)，但是其垂直精度較SRTM差(表1)，且易受到云層遮擋的影響出現(xiàn)空白，因此影響到了流域提取結(jié)果。

由此可見(jiàn)當(dāng)水平分辨率達(dá)到一定精度時(shí)，與水平分辨率相比，DEM的垂直精度對(duì)提取地形信息的真實(shí)性起決定作用。

表4 分片流域參數(shù)提取計(jì)算結(jié)果

3.3 流域參數(shù)提取

選擇流域上部走馬湖4條支流(新農(nóng)村港、鴨畈港、青山截留港和牌樓港)統(tǒng)計(jì)主河流長(zhǎng)度及比降參數(shù)，并與實(shí)測(cè)港道地形對(duì)比，結(jié)果見(jiàn)表4。從河流長(zhǎng)度上看，提取的主河流長(zhǎng)度均大于實(shí)測(cè)的港道，其中以牌樓港差異最大，通過(guò)DEM提取的牌樓港主河流長(zhǎng)度是其實(shí)際測(cè)量長(zhǎng)度的1.78倍。究其原因，實(shí)際測(cè)量的河流長(zhǎng)度是以是否有明顯水面作為判斷標(biāo)準(zhǔn)，而在設(shè)計(jì)洪水計(jì)算中河流長(zhǎng)度是河口至分水嶺的長(zhǎng)度，顯然其長(zhǎng)度會(huì)大于通常所講的河流長(zhǎng)度概念。采用自動(dòng)提取有一個(gè)好處就是避免了傳統(tǒng)的主觀判定主河流走向帶來(lái)的誤差。

通常而言，河流的上游比降大于下游比降，因此提取的主河流比降明顯大于測(cè)量的河流比降，新農(nóng)村港、牌樓港和青山截留港主河流比降分別為測(cè)量長(zhǎng)度比降的1.94、2.63和1.51倍。為了驗(yàn)證SRTM數(shù)據(jù)提取的主河流比降的精度，使用SRTM數(shù)據(jù)計(jì)算實(shí)測(cè)河段的比降，并與實(shí)際測(cè)量高程計(jì)算而得的比降對(duì)比。實(shí)際測(cè)量河流時(shí)，測(cè)量的為河床高程，而SRTM為地表高程即在河流處為水面高程，因此兩者比降不同(表4)。除了鴨畈港外，SRTM數(shù)據(jù)計(jì)算的河流比降略大于實(shí)測(cè)的比降，4條河段比降相差絕對(duì)值在-0.11‰～0.7‰?？傮w而言SRTM數(shù)據(jù)提取的河流比降精度較高，可滿足實(shí)際應(yīng)用，但是在人類(lèi)活動(dòng)較為明顯的青山截留港和鴨畈港誤差明顯，這是因?yàn)樵诖藚^(qū)域STRM數(shù)據(jù)時(shí)效性變差。

綜上所述，基于SRTM DEM數(shù)據(jù)的流域特征自動(dòng)提取方法能夠較為精確地提取主要的流域特征，尤其是對(duì)于主河流而言，避免了因主觀判定而導(dǎo)致錯(cuò)誤。河流比降的計(jì)算也較為精確，能夠滿足實(shí)際的應(yīng)用。

4 結(jié)論

本文針對(duì)無(wú)資料地區(qū)，提出了基于衛(wèi)星影像和公開(kāi)地形數(shù)據(jù)的平原湖泊流域特征自動(dòng)提取的方法，并在花馬湖流域進(jìn)行應(yīng)用，取得了比較好的效果，主要結(jié)論如下。

a) 在通過(guò)目視解譯高清衛(wèi)星影像獲取的真實(shí)河網(wǎng)信息輔助下，基于公共DEM數(shù)據(jù)能夠快速準(zhǔn)確地提取流域的河網(wǎng)水系結(jié)構(gòu)，并且在平原地區(qū)的精度有較大的提高。

b) 不同的數(shù)據(jù)源對(duì)于流域特征的提取影響較大。DEM的垂直精度對(duì)流域河網(wǎng)及邊界影響較大，相對(duì)于GDEM數(shù)據(jù)SRTM數(shù)據(jù)提取的結(jié)果更加可靠。

c) 基于SRTM DEM數(shù)據(jù)的流域特征自動(dòng)提取方法能夠較為精確地提取主要的流域特征，精度較高。

d) 人類(lèi)活動(dòng)尤其是城市化對(duì)于流域河網(wǎng)的結(jié)構(gòu)有著較為顯著的影響，同時(shí)在人類(lèi)活動(dòng)較為強(qiáng)烈的區(qū)域SRTM數(shù)據(jù)還存在時(shí)效性不足的特點(diǎn)。未來(lái)還需加強(qiáng)在城市區(qū)域流域特征提取方面的研究，此外如何結(jié)合多源觀測(cè)數(shù)據(jù)提高流域特征提取的精度也是進(jìn)一步研究的內(nèi)容。