柴汶河河流袭夺前水系的重建

朱晓洁, 褚忠信*, 于广科,毛士博, 姜瑾斐, 朱龙海

(1. 中国海洋大学 海洋地球科学学院, 山东 青岛 266100; 2. 海底科学与探测技术教育部重点实验室, 山东 青岛 266100)

山区河流的袭夺演化和分水岭的迁移对流域地貌演变发育过程起着重要的作用,同时,河流水系演化也是地球表层循环系统中的重要组成部分[1].分水岭迁移的结果是在一定的时间内使河道发生改道[2],进而引起河流发生袭夺.分水岭迁移通常多发生在不对称分水岭区域,随着时间的推移,下切侵蚀速率快、位置较低的一侧河流率先切入并侵蚀破坏分水岭,使分水线不断地向侵蚀能力较弱、位置较高的一侧河流移动,最终导致原先完整的分水岭断开并发生迁移[3-4],因此对断开的分水线进行连接恢复是重建袭夺前水系的关键环节.对于河流袭夺前水系的重建,已有学者根据最陡路径法[5]、流域面积插值法[6]等对其进行重建研究.

距今约8万年前,在柴汶河上游流域发育一起较为年轻的河流袭夺事件[6-9],其袭夺地貌特征保存完整，兼具充足的地貌学及沉积学证据,在全球范围内具有稀缺性与典型性.由于数字高程模型(DEM)囊括了丰富的地形、地貌及水文等空间分布信息[7],因此基于DEM数据对流域地形地貌特征进行提取已成为获取流域信息的主要手段[8].本文对柴汶河河流袭夺前的水系进行重建分析,以期有助于进一步了解河流袭夺理论、认识水系变迁规律、丰富流域地貌认知及规划河流的开发和治理.

1 研究区域概况

此研究区域位于山东省淄博市沂源县大张庄镇西南边缘,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 数据与方法

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 结果与分析

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 讨论

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 柴汶河上游部分河流纵剖面

5 结论

本文主要基于ArcGIS软件为技术支撑,结合野外调查,利用DEM基础数据对柴汶河袭夺河流域进行河网水系的提取以及河流袭夺前水系流域的恢复重建,取得了较好的结果.主要认识如下:

1) 在利用ArcGIS软件提取河网水系的过程中,基于不同集水面积阈值提取河网水系时会得到不同的结果,因此会对后续的水文分析产生影响.本文采用均值变点分析法确定了柴汶河流域的最佳集水面积阈值为0.54 km2,即设置汇流累积量为600时,提取的河网水系最接近于实际水系,进一步印证了均值变点分析法在确定最佳集水面积阈值时具有一定的可靠性.

2) 在对袭夺前流域的重建过程中,主要基于二次曲线内插法来恢复其袭夺前的分水岭,得到新的流域边界,完成袭夺前水系的重建;同时尚未考虑其区域岩性、地质构造、构造运动、气候变化等因素的影响.最终得到的结果说明这一想法应用在袭夺前河流的恢复重建方面是科学的,并且为袭夺前流域的恢复重建提供了重要进展.

3) 通过对柴汶河袭夺前水系流域的恢复重建以及流域地貌演化分析,河流袭夺对于被夺河与袭夺河的河谷形态发育产生了重要的影响,同时古柴汶河上游的河流纵比降大于古沂河上游,对地面侵蚀能力强,故能首先切穿分水岭,袭夺古沂河,该研究结果对位于较低的河流通过向上溯源侵蚀袭夺较高的邻近河流提供了又一例证.