基于DEM的西天山中部河流提取研究

杜青松，李國玉，李金明，周 宇

(1.中國科學院西北生態(tài)環(huán)境資源研究院凍土工程國家重點實驗室，蘭州 730000；2.中國科學院 西北生態(tài)環(huán)境資源研究院 大興安嶺凍土工程與環(huán)境觀測試驗研究站，黑龍江 加格達奇 165000；3.中國科學院大學，北京 100049)

0 引 言

河流作為一種降水的自然特征，與降雨量、地形等因素息息相關。數(shù)字高程模型(Digital elevation model，DEM)是一種地表地形精確有效的數(shù)字表達，包含豐富的地形、地貌和水文信息[1-5]，無洼地DEM數(shù)據(jù)是提取自然界河流的基礎數(shù)據(jù)[3,6,7]?；贏rcGIS的水文分析擴展模塊，可以快速提取自然界的河流信息[3, 5]，這一方法是目前河網(wǎng)提取最行之有效的方法之一。自O′callagehan和Mark提出河網(wǎng)識別法[8-10]以來，學者們通過坡面流模擬，設定合適閾值快速提取自然界河網(wǎng)這一方法在國內(nèi)外已開展大量工作[11-13]并取得大量成果[1,14,15]。通過DEM數(shù)據(jù)提取河網(wǎng)信息不僅適用于自然山區(qū)[14]，對喀斯特地貌的小流域提取[15]也同樣適用。但在平坦地形地區(qū)，該方法提取的河網(wǎng)不理想，存在偽河流現(xiàn)象。

基于DEM數(shù)據(jù)提取河網(wǎng)的難點在于集水(閾值)面積的確定，對于同一流域利用不同集水面積閾值得到的河網(wǎng)不同[16]，合適的集水面積對準確河流地提取至關重要。集水閾值的確定方法大致可分為試錯法[17, 18]、數(shù)理統(tǒng)計法[3,14,19-21]、地形分維數(shù)法[22, 23]等。各方法在不同區(qū)域均得到應用，但都存在一定的人為主觀因素影響。其中，由于試錯法以定性的方式去描述提取的河網(wǎng)，需要大量的經(jīng)驗，主觀性較強，產(chǎn)生的誤差最大。

本文在綜合前人研究的基礎上，基于數(shù)理統(tǒng)計方法，以試錯的方式提取不同集水面積下的河網(wǎng)，擬合集水面積閾值-河網(wǎng)密度曲線，采用當河網(wǎng)密度變化與集水面積變化相等[24]時求解方程確定集水閾值，用該閾值提取研究區(qū)的河網(wǎng)信息，將提取的河網(wǎng)與實際河流對比作精度驗證，計算河流的霍頓參數(shù)并對該研究區(qū)域內(nèi)礦區(qū)通過河流對周圍環(huán)境造成的影響做簡要分析。

1 研究區(qū)概況

1.1 研究區(qū)邊界確定

在傳統(tǒng)方法中，對于研究區(qū)邊界的確定一般是以行政單位為邊界，或直接畫一個包含研究目標物的矩形框作為研究區(qū)的邊界。河流是一種自然要素，為了保證其本身的自然屬性，通過上訴的兩種方法均不能得到自然界完整的流域，因此為了保證研究區(qū)河網(wǎng)的完整性，本文對研究區(qū)的邊界獲取通過“盆域”工具獲取?！芭栌颉惫ぞ咄ㄟ^DEM數(shù)據(jù)和流向數(shù)據(jù)，可以提取河流完整的流域邊界，將通過“盆域”分析后得到的柵格矢量化可得到研究區(qū)的邊界。

1.2 研究區(qū)自然概況

研究區(qū)位于西天山中部，行政區(qū)域隸屬于新疆巴音郭楞蒙古自治州和靜縣。地理位置(見圖1)42°59′59.96″～43°34′14.13″N，84°51′22.21″～85°28′33.32″E。面積1 280.84 km2,海拔2 462～5 179 m，平均海拔3 406.79 m，屬于高寒高海拔地區(qū)?，F(xiàn)場觀測資料表明該區(qū)年最高溫度16 ℃，最低氣溫-30 ℃左右，日最大降水量為 146 mm，年降雨量超過1 000 mm，蒸發(fā)量425 mm，最大風速 12 m/s，平均濕度約 43%，風向以NEE為主。地表覆蓋以高山草甸為主，以草本為主，近流域范圍內(nèi)存在少量灌木，無高大喬木。海拔3 700 m向上為冰原帶，幾乎無植被生長。草原是下游牧民放牧的根基，草生長的好壞直接決定牧民的收入水平，草甸生長需要的水、畜牧用水主要來源于高海拔地區(qū)的冰雪融水。

圖1 研究區(qū)DEM分布圖Fig.1 The DEM data of study area

研究區(qū)被一條分水嶺分為兩個流域，分水嶺以北的瑪納斯河流域和以南的眾多細小河流流域，最終匯入伊犁河。分水嶺處有一個大型露天鐵礦(新疆諾爾湖鐵礦)，鐵礦由位于河流上游的礦區(qū)、選礦區(qū)和位于河流下游的生活區(qū)組成。在礦區(qū)的東北方向有一個自然湖泊諾爾湖，湖泊里的水由西向東注入瑪納斯河。河水主要來源于冰雪融水和位于地表下的石冰川和凍土融化，河流汛期為每年8月，10月后溫度降低，凍土凍結，活動層變淺，礦區(qū)以固態(tài)降水為主，河流進入枯水期。

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)源及數(shù)據(jù)預處理

利用SRTM v4 DEM數(shù)據(jù)(http:∥srtm.csi.cgiar.org/SELECTION/inputCoord.asp)，對下載的數(shù)據(jù)進行數(shù)據(jù)鑲嵌、坐標轉(zhuǎn)換、數(shù)據(jù)裁剪等預處理，得到與研究區(qū)大小一致的DEM原始影像。SRTM v4 DEM作為最常用的地表數(shù)字高程之一，被廣泛用于河網(wǎng)提取[25, 26]研究。

由于受空間分辨率和系統(tǒng)誤差的影響，DEM數(shù)據(jù)在生成過程中，會產(chǎn)生一些的錯誤的“洼地”，這會造成在河流流向提取時，DEM中洼地水流方向的不正確，進而影響河網(wǎng)提取的準確性。因此需要對原始DEM進行“填洼”處理以消除洼地的影響。但并不是所有的洼地都是由數(shù)據(jù)誤差造成的，有些洼地是真實地形的反應，如圖1中諾爾湖所在的位置，基于此理論，在對原始DEM作“填洼”處理時要設置合理的填洼閾值。

填洼閾值可通過實際計算得到，具體步驟分為：①對原始DEM應用“流向”工具計算流向數(shù)據(jù)；②通過“匯”工具流向數(shù)據(jù)為輸入數(shù)據(jù)計算DEM中的所有洼地；③利用“分水嶺”工具提取每個洼地的貢獻區(qū)域；④使用“分區(qū)統(tǒng)計”工具計算出每個洼地貢獻區(qū)域的最低高程和最低高程；⑤運行“柵格計算器”工具將最、最低高程值作差得到洼地填充閾值；⑥最后操縱“填洼”工具對原始DEM進行 “填洼”迭代運算，得到無洼地DEM。

2.2 研究方法

本文首先設置不同的集水面積閾值，并提取每一集水面積對應的河網(wǎng)信息，計算出其河網(wǎng)密度、河流長度、流域面積等，然后對每組數(shù)據(jù)作相關性分析，擬合合適的模型。最后借鑒前人研究，采用河網(wǎng)密度變化率等于集水值變化率[24]的數(shù)值方法確定該區(qū)域正確的集水閾值，該方法最大優(yōu)點是可以有效地避免人為主觀性的影響，求解的集水面積為唯一確定值。其原理是：隨著集水面積的增大，位于坡面的河道起點會向流域內(nèi)相鄰的地勢平坦區(qū)域靠近，位于坡面上的河道數(shù)減少，偽河道也會被刪除[1]。將河網(wǎng)密度與集水閾值的冪函數(shù)同切線方程聯(lián)立，求解方程組的唯一解作為河網(wǎng)提取的閾值，此時的閾值為最佳河流積水閾值[14, 24]。河流匯流面積通過公式計算：

S=NPix

(1)

式中：S為匯流面積，單位由柵格分辨率決定，小流域單位為m2，大流域單位為km2；N為流量柵格數(shù)；Pix為柵格分辨率。

通過ArcGIS10.6軟件的水文工具箱分別提取不同集水面積下的河網(wǎng)，在等面積投影坐標系下利用幾何計算器統(tǒng)計矢量化后的河網(wǎng)長度和流域面積等，再利用統(tǒng)計工具OriginLab計算河網(wǎng)密度，構建集水面積-河網(wǎng)密度曲線。將擬合的函數(shù)與切線方程聯(lián)立，求解方程組的唯一解作為河網(wǎng)的最佳集水面積。將這一集水面積下提取的河網(wǎng)疊加在Google Earth影像上，與實際河流對比以驗證其提取的準確性。方法流程圖(圖2)，具體操作步驟如下：

圖2 方法流程圖Fig 2 Method flow

(1)利用水文工具箱的“填洼”工具對原始DEM進行填洼處理，得到無洼地DEM；

(2)通過“流向”工具，采用D8算法提取河流方向；

(3)應用“流量”工具提取河流累計柵格河網(wǎng)柵格；

(4)在“柵格計算器”工具設置不同的匯流柵格數(shù)(條件函數(shù)等)，提取對應的柵格河網(wǎng)，同時計算該柵格數(shù)對應的集水面積；

(5)將(4)得到的柵格河流利用“河流鏈接”和“河流分級”工具進行河網(wǎng)分級處理，分級采用斯特拉勒(STRAHLER)分級法；

(6)輸入分級的河網(wǎng)數(shù)據(jù)，通過“河網(wǎng)矢量化”工具將柵格河網(wǎng)轉(zhuǎn)為矢量數(shù)據(jù)；

(7)使用“分水嶺”工具提取流域面積，同時將得到的結果矢量化；

(8)計算河網(wǎng)密度，擬合集水面積-河網(wǎng)密度曲線，聯(lián)立方程，求解集水閾值；

(9)提取最優(yōu)集水面積對應的河網(wǎng)，驗證其精度。

3 結果與分析

3.1 不同集水閾值的河網(wǎng)提取

ArcGIS水文分析中的流量定義與水文學中的流量含義不一樣，它不代表實際流量，該流量僅是流量分析結果的一個柵格累計計數(shù)。流量柵格中每一個像元記錄的是流向該點的柵格數(shù)量總和，流域集水面積與流量柵格數(shù)的關系為式1，通過該式可實現(xiàn)流量柵格數(shù)與集水面積的相互轉(zhuǎn)換。

以500為間隔，分別設置流量柵格數(shù)為500～8 500 共計17組數(shù)據(jù)，對應集水面積為0.45～7.65 km2之間的以0.45 m2為間隔的17個數(shù)據(jù)，同時提取各集水面積對應的河網(wǎng)，統(tǒng)計分析各集水閾值與河流數(shù)、河流長度、河流流域面積和河網(wǎng)密度之間的關系。

集水面積與河流數(shù)的關系如圖3所示，集水面積與總河流數(shù)和一級河流數(shù)相關性較好，與其他級別河流數(shù)相關性不明顯。隨著集水閾值的增加，各級河流數(shù)和總河流數(shù)均減少，一開始降幅大，在集水閾值為5.4 km2時，總河流數(shù)和一級河流條數(shù)趨于穩(wěn)定；在集水面積為1.35 km2時，最高河流等級由5級減少為4級，總河道數(shù)由2 039 減少到628條，減少69.2%；當集水面積超過5.85 km2后，最高級的4級河流條數(shù)穩(wěn)定為3條，不再隨集水面積的增大而變化。

圖3 集水閾值(面積)與河流數(shù)關系Fig.3 Relationship between catchment threshold (area) and number of rivers

圖4描述的是集水面積與河流長度的關系，分析得出，集水面積與河流長度的關系和集水面積與河流條數(shù)的規(guī)律類似，總河流長度和一級河流長度受集水面積影響大，集水面積對其他各級河流長度影響小，隨著集水面積的變化，1～4級河流長度變化不明顯，趨于一個穩(wěn)定值?？偟膩碚f，河流長度隨著集水閾值的增大總體呈減小趨勢，但減小幅度慢慢變緩；當集水面積超過5.85 km2時，4級河流長度保持3.79 km不變，不再隨集水面積的增大而變化，此時與河道數(shù)保持不變相對應。

圖4 集水閾值(面積)與河長數(shù)的關系Fig.4 Relationship between catchment threshold (area) and river length

從圖5可以看出，集水面積與流域面積的相關性不高，變化規(guī)律不明顯，除一級河流流域面積總體上隨著集水面積的增大而減小外，其他各級流域?qū)牧饔蛎娣e變化無規(guī)則，相關性極差，各級河流流域面積變化情況也不一致，有增有減。

圖5 集水閾值(面積)與流域面積的關系Fig.5 Relationship between catchment threshold (area) and watershed area

河網(wǎng)密度為單位面積上的河流長度，是河流的一個基本參數(shù)，與河流閾值變化息息相關。圖6顯示集水面積-河網(wǎng)密度曲線，可以得出，集水面積與河網(wǎng)密度的相關性極高(R2=0.999 78)，河網(wǎng)密度隨著集水面積的增大而減小，最終趨于穩(wěn)定。集水面積-河網(wǎng)密度呈現(xiàn)為冪函數(shù)圖像(圖6)，函數(shù)表達式為：

y=0.828 24x-0.488 93(R2=0.999 8)

(2)

式中：y為河網(wǎng)密度，km-1；x為集水面積，km2。

圖6 集水閾值(面積)-河網(wǎng)曲線Fig.6 Catchment threshold (area)-river network curve

3.2 集水面積閾值的確定

采用河網(wǎng)密度變化與集水面積變化相等時的數(shù)值方法確定最佳集水面積閾值。將式子y=-x+a與式(2)聯(lián)立，求解方程組，得出的唯一解x為最佳閾值。

(3)

求解方程組(3)，有a=1.659 4時，方程組的未知數(shù)x存在唯一解0.544 9，由此確定研究區(qū)最佳集水面積為0.544 9 km2，計算得出對應的流量累積柵格數(shù)為605。同時提取該集水面積閾值下的河網(wǎng)，結果經(jīng)矢量化后如圖7所示。

圖7 最佳集水面積下提取的河網(wǎng)Fig.7 River network extracted under optimal catchment area

以最優(yōu)集水面積0.544 9 km2提取的河網(wǎng)，河流由5種級別組成，1級河流為最基礎河流，5級為最高等級。提取的河網(wǎng)存在大量的基礎河流，這與研究區(qū)的高海拔、起伏大的地形地貌特征相符，即使是在相對平坦的區(qū)域也幾乎不存在平直的偽河流，這一結果與實際調(diào)研情況也相符合。

通過河流與礦區(qū)、生活區(qū)等的空間拓撲關系及河流的流向可知：礦區(qū)開采使用大量的炸藥、運輸車輛存在漏油、燃油加鉛等現(xiàn)象，采礦會對環(huán)境造成污染，是一大的污染源。經(jīng)過河流的作用會將這一點源污染擴展為面源污染。其中采礦區(qū)對諾爾湖沿線的水體造成污染嚴重并通過河流將污染物帶入瑪納斯河流域，生活區(qū)和下游草地主要受選礦區(qū)影響。河流流向?qū)Νh(huán)境污染研究采樣點的確定有指導意義，通過對河流沿線和污染源地區(qū)采樣對比研究，可大致了解研究區(qū)域的污染特征。

河流霍頓參數(shù)[22]是表征河流與地貌關系的基本參數(shù)，在水文模型模擬研究中具有重大意義，是提取是水文模型構建的基礎[19, 22]，研究區(qū)的河流基礎數(shù)據(jù)和地貌特征參數(shù)如表1所示。

研究區(qū)總河流共計1 661 條，基礎河流發(fā)育，占50%；一級河流面積占63%；河網(wǎng)密度為1.11 km-1，河源密度0.65 條/km2；河流的河數(shù)率、河長率、面積率分別為1.66、1.83和2.12。該區(qū)域由于海拔高、地形起伏大，有利于河流的形成；同時，存在大量冰雪、凍土和石冰川，受氣溫影響大，故存在大量由固態(tài)水融化而形成的細小河流。

3.3 提取河網(wǎng)的精度驗證

研究區(qū)范圍較小，缺乏能與之相匹配的地表基礎數(shù)據(jù)，如河流數(shù)據(jù)等。故很難采用傳統(tǒng)方法將提取的河網(wǎng)與“藍線河網(wǎng)”對比定量研究[21, 27, 28]得出河流提取的誤差。為驗證河網(wǎng)提取精度，將提取的河流與Google影像疊加(圖8)，對比其與Google影像中實際河流之間的誤差。通過對比發(fā)現(xiàn)提取的河網(wǎng)與實際河流基本一致，幾乎不存在偽河網(wǎng)或河網(wǎng)提取不足的情況。

表1 研究區(qū)河網(wǎng)特征參數(shù)Tab.1 Characteristic parameters of river network in study area

圖8 河網(wǎng)與Google影像對比Fig.8 River network in Google image

將提取的矢量河網(wǎng)轉(zhuǎn)為kml格式可加載到Google Earth Pro中(圖9)更好地從多視角、多尺度、多維度(2D、3D)對比，判斷河網(wǎng)與Google地形圖是否一致從而判定河網(wǎng)提取精度。結果表明，提取的河網(wǎng)與Google地形契合度很高，河流分布符合地形規(guī)律，提取的河網(wǎng)與真實河網(wǎng)具有一致性。

圖9 河網(wǎng)與谷歌地形對比Fig.9 Comparison of river network and Google terrain

經(jīng)過將提取的河網(wǎng)與Google 影像和Google地形圖對比研究，得出在最優(yōu)閾值(0.544 9 km2)提取的河網(wǎng)與實際河網(wǎng)誤差不大，具有很高的一致性。提取的河網(wǎng)可以作為研究區(qū)流域信息提取的基礎數(shù)據(jù)，進而可以計算更多的河流特征參數(shù)，從而實現(xiàn)對流域生態(tài)環(huán)境的深入研究。

4 結 論

本文基于DEM數(shù)據(jù)提取河網(wǎng)信息，采用數(shù)值計算方法解決了研究區(qū)基礎河流數(shù)據(jù)缺乏的問題。

(1)確定了該研究區(qū)最優(yōu)流域集水面積為0.544 9 km2，獲得的河網(wǎng)數(shù)據(jù)真實可靠，提取的河網(wǎng)數(shù)據(jù)可作為輸入數(shù)據(jù)進行下一步的研究，完善了該研究區(qū)域的基礎數(shù)據(jù)庫；

(2)在河流的眾多特征參數(shù)中，河網(wǎng)密度對河流集水面積變化最敏感；各級流域面積變化與集水面積變化不相關；

(3)礦區(qū)開采直接影響諾爾湖的水質(zhì)情況，從而可能會對瑪納斯河造成污染；

(4)選礦區(qū)污染物可能會隨河流擴散對下游的工人生活區(qū)及草原造成影響。

□