朱曉潔, 褚忠信*, 于廣科,毛士博, 姜瑾斐, 朱龍海
(1. 中國海洋大學 海洋地球科學學院, 山東 青島 266100; 2. 海底科學與探測技術教育部重點實驗室, 山東 青島 266100)
山區河流的襲奪演化和分水嶺的遷移對流域地貌演變發育過程起著重要的作用,同時,河流水系演化也是地球表層循環系統中的重要組成部分[1].分水嶺遷移的結果是在一定的時間內使河道發生改道[2],進而引起河流發生襲奪.分水嶺遷移通常多發生在不對稱分水嶺區域,隨著時間的推移,下切侵蝕速率快、位置較低的一側河流率先切入并侵蝕破壞分水嶺,使分水線不斷地向侵蝕能力較弱、位置較高的一側河流移動,最終導致原先完整的分水嶺斷開并發生遷移[3-4],因此對斷開的分水線進行連接恢復是重建襲奪前水系的關鍵環節.對于河流襲奪前水系的重建,已有學者根據最陡路徑法[5]、流域面積插值法[6]等對其進行重建研究.
距今約8萬年前,在柴汶河上游流域發育一起較為年輕的河流襲奪事件[6-9],其襲奪地貌特征保存完整,兼具充足的地貌學及沉積學證據,在全球范圍內具有稀缺性與典型性.由于數字高程模型(DEM)囊括了豐富的地形、地貌及水文等空間分布信息[7],因此基于DEM數據對流域地形地貌特征進行提取已成為獲取流域信息的主要手段[8].本文對柴汶河河流襲奪前的水系進行重建分析,以期有助于進一步了解河流襲奪理論、認識水系變遷規律、豐富流域地貌認知及規劃河流的開發和治理.
此研究區域位于山東省淄博市沂源縣大張莊鎮西南邊緣,117°56′E~118°2′E,35°56′N~36°3′N之間,地處魯中腹地,西北毗鄰濟南市萊蕪區,西南與新泰市、蒙陰縣相連,東北地處沂源縣境內,屬4縣(市)交界位置[9].晚更新世時期,海拔較低的黃河支流——古柴汶河向上溯源侵蝕襲奪了海拔較高的淮河支流——古沂河[6],使古沂河流向由東北急劇向西南作了近90°的大拐彎,在形態上呈倒勾狀彎曲,因此古沂河上游匯入古柴汶河,形成現如今的柴汶河流域,如圖1所示.

圖1 柴汶河流域位置
柴汶河流域[10]源于沂源縣西南部牛欄峪一帶,東流至周科峪后西折向西南流,穿新泰市東周水庫,經新泰、寧陽縣境至大汶口入主流牟汶河,其干支流都是源短流急的山洪河流,洪水漲落迅猛,平時只有涓涓細流.流域內地勢東高西低,北高南低,北有泰山,東靠魯山和蒙山,西、南為丘陵和平原.柴汶河全長116 km,流域面積1 944 km2,是汶河的最大支流.
2.1 數據來源本文采用的實驗數據由地理空間數據云提供(http://www.gscloud.cn/),為ASTER GDEM-DEM數據,空間分辨率為30 m×30 m,空間參考為WGS_1984_UTM_Zone_50N.對所獲得的DEM數據利用ArcGIS 10.2軟件進行拼接裁剪處理,得到河流襲奪研究區的原始DEM數據,面積約為80 km2.同時,在原始DEM數據上疊加部分水系流域,得到柴汶河上游流域水系分布,如圖2所示.

圖2 柴汶河上游流域附近水系分布
2.2 最佳集水面積閾值確定將河流襲奪研究區的原始DEM數據加載到ArcGIS軟件中,對原始的DEM進行填洼[11]、流向提取、匯流累積量計算[12]以及河網水系的提取[13-15],再通過均值變點分析法[16]確定與實際水系接近的河網水系[17].均值變點分析法目前主要應用于提取地形起伏度的最佳單元[18-21],同時也逐漸應用于確定河流的最佳集水面積閾值[22].
2.3 分水線提取與恢復分水線通常為分水嶺最高點之間的連線,同時也為水流的起源點[23].因此,對DEM數據進行地表徑流模擬計算時,通過提取匯流累積量的零值可得到分水線[24].
斷開的分水線在DEM數據上主要表現為高程值由高變低,因此將分水線的高程值作為恢復襲奪前分水嶺的單一因子,即通過抬高斷開的分水線的高程值恢復襲奪前的分水嶺.已知斷開分水線兩端點的坐標及高程值(x1,y1,z1)、(x2,y2,z2),利用C++編寫程序對斷開的分水線之間的柵格依次重新賦予高程值,其原理如圖3所示.

圖3 分水線高程值重新賦值原理
首先,基于均值替換法[25],計算被連接的分水線兩端點(A、B)高程的平均值(即中點值)代替襲奪前的原中點(C)高程值;再基于二次曲線內插法相似原理,設Zi為被連接的分水線對應的各柵格單元高程值,可得到某一柵格單元(D、E)高程值的計算公式:
Z
(1)
Z
(2)
3.1 河網水系提取基于ArcGIS水文分析模塊,分別設置匯流累積量閾值自小到大為100、200、300、…、2 600共26組數據進行水系的提取,以100、600和1 600匯流累積量閾值為例,僅提取柴汶河上游流域河網水系,如圖4(a)~(c)所示.根據結果可知,隨著匯流累積量閾值的增大,提取出來的河網水系逐漸稀疏,河流總長度逐漸變短,水系分級逐漸不詳細,在后續研究分析時會對結果的準確度產生影響[26].因此,為提高結果準確度,采用均值變點分析法確定當集水面積閾值為0.54 km2時(即匯流累積量閾值為600時),根據提取出來的水系與1∶5 000的實際河網水系(如圖4(d)所示)進行對比發現較為吻合,最接近于實際河網水系.


3.2 分水線恢復河流襲奪與分水嶺遷移是相輔相成的2個過程,河流襲奪必然引起分水嶺的遷移;同樣,分水嶺遷移會引起水系重組,進而發生河流襲奪.分水嶺是分隔2個流域的界線,是地球內外營力共同作用的產物.當分水嶺兩側的水系侵蝕能力不同時,分水嶺會向侵蝕能力弱的一側遷移,當地貌演化至某一階段時,侵蝕能力強的水系會在某處率先侵蝕切穿分水嶺,并匯入侵蝕能力弱的水系,進而發生流域重組.
因此,斷開的分水嶺在DEM數據上最直觀的表現是:某處被切穿的分水嶺的高程明顯低于兩端的高程,在地勢上呈“凹”字形態,如圖5(a)所示.根據均值變點分析法確定的最佳集水面積閾值提取河網水系,并與分水線圖層進行疊加,同時結合野外地形實地考察分析,判斷確定需要重新被連接的分水線,如圖5(b)和(c)所示.

(a) 斷開的分水嶺示意圖 (b) 無人機航拍襲奪灣位置 (c) 需連接的分水線位置
圖 5 柴汶河上游襲奪區域分水線恢復連接位置
Fig. 5 The connection position of the dividing line of the captured area in the upper reaches of the Chaiwen River
根據分水線恢復原理,對被連接分水線途經柵格單元的高程值重新賦值.圖6中方格為柵格單元,數值為高程值.圖6(a)中填充的深色方格表示斷開的分水線兩端端點,斜直線表示需要被重新連接的分水線,斜直線途經的方格表示需要被重新賦值的柵格單元.由圖6(a)可知,在未恢復斷開的分水線之前,斷開的分水線兩端點高程值均高于需要被重新連接的分水線所途徑的高程值,呈“兩端高、中間低”的形態,證明本文判斷的需要被重新連接的分水線是正確的.圖6(b)中填充的深色方格表示已被重新賦值的柵格單元,方格內的數值為賦值后的新高程值.將分水線恢復前后的DEM高程值對比可知,原先斷開的分水線的高程值已由最初的低值重新賦值為高值,在恢復分水線前的高程值上抬高了約0.3~15.5 m,從斷開的分水線一端的高程值378 m逐漸抬高,最終逐漸接近另一端點高程值383 m.至此,柴汶河流域河流襲奪前的分水線已被成功恢復.同時,基于上述原理重新生成新的DEM,即柴汶河河流襲奪前的DEM數據.

(a) 分水線恢復前的DEM值 (b) 分水線恢復后的DEM值圖 6 柴汶河上游襲奪位置分水線恢復前后的DEM值
3.3 襲奪前水系流域恢復河流襲奪是引起水系結構改變,進而誘發流域盆地重組的重要原因之一.恢復襲奪前的水系流域有助于進一步了解河流襲奪產生的機理,以及河流襲奪對地貌產生的影響.流域的邊界線通常是相鄰兩個流域之間的分水線,因此根據恢復分水線成功后得到的河流襲奪前DEM數據,對其進行水系提取和流域劃分.將襲奪前的流域劃分為被奪河流域與襲奪河流域,并疊加新提取的水系圖層,從而得到柴汶河流域河流襲奪前的水系分布與流域劃分(如圖7(a)所示).為了便于對比與分析,同時提取并劃分了襲奪后的水系及流域(如圖7(b)所示).

(a) 襲奪前水系流域劃分 (b) 襲奪后水系流域劃分圖 7 柴汶河上游襲奪前后水系流域劃分
通過對襲奪前后的水系流域進行對比分析,并根據河流襲奪前水系流域分布圖中的水系形態可以推斷得出:在河流襲奪發生前,被奪河古沂河自西南向東北方向流經大張莊、田莊水庫后最終匯入沂河干流;河流襲奪發生后,古沂河在襲奪灣處發生了近90°的大轉彎后匯入了古柴汶河,襲奪的同時在襲奪灣生成具有高勢差的裂點,裂點一分為二,分別沿被奪河、倒流河向上遷移.倒淌河因失去源頭供給,且裂點遷移對其河道的調整,使流域內地勢發生倒轉,呈西南低東北高,遂則改向匯入古柴汶河,成為柴汶河流域的一部分.同時,河流襲奪也產生了斷頭河(如圖2所示),并與倒淌河之間形成了新的分水高地(即風口),將流域分為柴汶河流域和沂河流域.
4.1 河流襲奪與流域地貌形態演變關系流域地貌形態是由于地表受到侵蝕和堆積形成的,其形態的差異反映了流域被侵蝕的程度以及流域地貌演化的階段[27].河谷的形成與發育問題是地質地貌學長期關注的經典課題,有關河谷發育的學說已成為地貌學發展的支持基石和理論源泉[28].在山地和丘陵地區,由于地質構造復雜、巖層堅硬,河谷的形態主要由地質因素決定,表現為寬谷與峽谷沿河交替分布,河谷發展過程在基巖山地河谷的橫剖面發展過程中最為清楚.
以襲奪灣為界,分別選取古沂河與古柴汶河的上游部分河谷作橫剖面,剖面位置如圖7(b)所示.從橫剖面(如圖8所示)A-A’、B-B’和C-C’中可以看出,古沂河上游的河谷剖面為不對稱“V”字形.其河谷的凹岸發生了侵蝕,坡面較陡;凸岸逐漸堆積,坡面較緩,且出現了河灣和水小谷寬的現象.這是由于河流發生襲奪后,水流量減少,河流的下蝕能力減弱,側蝕能力加強,呈現出原河谷不相適應的現象.

圖8 柴汶河上游部分河谷橫剖面
在古柴汶河的上游,以下切侵蝕作用為主,河流深切入基巖,形成河身直,河床坡度陡,斷面狹窄的“V”型谷.剖面D-D’為“V”型谷的初級發育狀態——嶂谷,兩壁陡峭,常有基巖或礫石灘露出水面以上.嶂谷進一步發展形成峽谷(剖面E-E’),橫剖面呈明顯的“V”字形,出現谷中谷現象,這是因為在河流襲奪未發生前,該段河谷河流流量較小,下蝕較弱,以側蝕為主,從而形成較寬的河谷地貌;在襲奪發生后,被奪河的河水匯入,使河流流量加大,流速變快,河流溯源侵蝕和下蝕較為強烈,從而使河床加深變窄,形成落差較大的谷中谷地貌.F-F’橫剖面處的河谷形態“V”字形逐漸展寬,河谷曲線光滑,有從“V”型谷逐漸向“U”型谷演變的趨勢.
4.2 河流襲奪與河流縱剖面關系河流縱剖面是河流作用形成的,河流為了適應外在條件的變化,每一時刻任一河段不僅進行侵蝕,同時也發生堆積,因此河流縱剖面存在自動調整作用[29].同樣,以襲奪灣為界,分別提取古沂河和古柴汶河上游的河流縱剖面(如圖9所示).結果可知,山區相較于平原地區的河流縱剖面較陡,起伏較大,多跌水和瀑布,兩岸常有許多山嘴突出,使河床岸線犬牙交錯,呈凹凸不規則型.河流縱比降不僅反映出流經該河床的水流流速,還反映巖性的變化或構造運動,同時也能間接反映河流對地面侵蝕的能力.通過計算分別得出被奪河的河流縱比降為10.7‰,襲奪河的河流縱比降為15.5‰,古柴汶河上游的河流縱比降大于古沂河上游.當流域面積不大且其他條件相當時,河流縱比降越大,水流速度就越大,代表水動能就越大,河流對地面的下切侵蝕能力就越強[30].古柴汶河上游對地面的侵蝕能力要大于古沂河上游,因此能率先切穿分水嶺,襲奪古沂河上游.結合上述對其襲奪前水系流域的恢復研究可知,該研究結果是對海拔較低、侵蝕能力強的古柴汶河向上溯源侵蝕襲奪了海拔較高、侵蝕能力弱的古沂河的又一佐證.

圖9 柴汶河上游部分河流縱剖面
本文主要基于ArcGIS軟件為技術支撐,結合野外調查,利用DEM基礎數據對柴汶河襲奪河流域進行河網水系的提取以及河流襲奪前水系流域的恢復重建,取得了較好的結果.主要認識如下:
1) 在利用ArcGIS軟件提取河網水系的過程中,基于不同集水面積閾值提取河網水系時會得到不同的結果,因此會對后續的水文分析產生影響.本文采用均值變點分析法確定了柴汶河流域的最佳集水面積閾值為0.54 km2,即設置匯流累積量為600時,提取的河網水系最接近于實際水系,進一步印證了均值變點分析法在確定最佳集水面積閾值時具有一定的可靠性.
2) 在對襲奪前流域的重建過程中,主要基于二次曲線內插法來恢復其襲奪前的分水嶺,得到新的流域邊界,完成襲奪前水系的重建;同時尚未考慮其區域巖性、地質構造、構造運動、氣候變化等因素的影響.最終得到的結果說明這一想法應用在襲奪前河流的恢復重建方面是科學的,并且為襲奪前流域的恢復重建提供了重要進展.
3) 通過對柴汶河襲奪前水系流域的恢復重建以及流域地貌演化分析,河流襲奪對于被奪河與襲奪河的河谷形態發育產生了重要的影響,同時古柴汶河上游的河流縱比降大于古沂河上游,對地面侵蝕能力強,故能首先切穿分水嶺,襲奪古沂河,該研究結果對位于較低的河流通過向上溯源侵蝕襲奪較高的鄰近河流提供了又一例證.