杜青松，李國(guó)玉，李金明，周 宇

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

0 引 言

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

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

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

1 研究區(qū)概況

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

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

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

研究區(qū)位于西天山中部，行政區(qū)域隸屬于新疆巴音郭楞蒙古自治州和靜縣。地理位置(見(jiàn)圖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ū)。現(xiàn)場(chǎng)觀測(cè)資料表明該區(qū)年最高溫度16 ℃，最低氣溫-30 ℃左右，日最大降水量為 146 mm，年降雨量超過(guò)1 000 mm，蒸發(fā)量425 mm，最大風(fēng)速 12 m/s，平均濕度約 43%，風(fēng)向以NEE為主。地表覆蓋以高山草甸為主，以草本為主，近流域范圍內(nèi)存在少量灌木，無(wú)高大喬木。海拔3 700 m向上為冰原帶，幾乎無(wú)植被生長(zhǎng)。草原是下游牧民放牧的根基，草生長(zhǎng)的好壞直接決定牧民的收入水平，草甸生長(zhǎng)需要的水、畜牧用水主要來(lái)源于高海拔地區(qū)的冰雪融水。

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

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

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

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

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

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

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

2.2 研究方法

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

S=NPix

(1)

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

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

圖2 方法流程圖Fig 2 Method flow

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

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

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

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

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

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

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

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

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

3 結(jié)果與分析

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

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

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

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

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

圖4描述的是集水面積與河流長(zhǎng)度的關(guān)系，分析得出，集水面積與河流長(zhǎng)度的關(guān)系和集水面積與河流條數(shù)的規(guī)律類似，總河流長(zhǎng)度和一級(jí)河流長(zhǎng)度受集水面積影響大，集水面積對(duì)其他各級(jí)河流長(zhǎng)度影響小，隨著集水面積的變化，1～4級(jí)河流長(zhǎng)度變化不明顯，趨于一個(gè)穩(wěn)定值。總的來(lái)說(shuō)，河流長(zhǎng)度隨著集水閾值的增大總體呈減小趨勢(shì)，但減小幅度慢慢變緩；當(dāng)集水面積超過(guò)5.85 km2時(shí)，4級(jí)河流長(zhǎng)度保持3.79 km不變，不再隨集水面積的增大而變化，此時(shí)與河道數(shù)保持不變相對(duì)應(yīng)。

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

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

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

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

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í)的數(shù)值方法確定最佳集水面積閾值。將式子y=-x+a與式(2)聯(lián)立，求解方程組，得出的唯一解x為最佳閾值。

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

□