陳 浩， 董廷旭， 李 勇， 文星躍

（1.綿陽師范學院 資源環(huán)境工程學院，四川 綿陽621006； 2.成都理工大學 油氣藏地質(zhì)及開發(fā)工程國家重點實驗室，四川 成都610059；3.西華師范大學 國土資源學院，四川 南充637002）

水系對于構(gòu)造、氣候等外部變量因素的改變非常敏感，水系演化通常記錄了造山帶系統(tǒng)新近時期的細微變動，因而，地質(zhì)構(gòu)造活躍區(qū)域水系發(fā)育模式及其對構(gòu)造活動、氣候條件的響應(yīng)機制一直以來是構(gòu)造地貌與河流地貌研究的熱點［1］.大致呈南北向展布的岷江上游流域盆地，位于第四紀以來活動較為強烈的岷山構(gòu)造帶西緣和龍門山構(gòu)造帶北側(cè)，雖然一些學者對岷江水系特征進行了較為深入的分析，但研究重心主要集中于岷江下切過程、河流階地樣式與區(qū)域隆升的耦合關(guān)系［2-10］，有關(guān)岷江上游流域地貌發(fā)育期次的研究目前還未見報道.近年來，地理信息系統(tǒng)技術(shù)的不斷發(fā)展和數(shù)字高程模型（DEM）數(shù)據(jù)的應(yīng)用使地貌過程研究由定性向定量化階段轉(zhuǎn)變［11］.基于DEM數(shù)據(jù)提取符合實際狀況的水系網(wǎng)絡(luò)是岷江上游流域地貌發(fā)育階段研究的前提，最佳集水面積閾值的確定是水系網(wǎng)絡(luò)提取過程中的關(guān)鍵環(huán)節(jié)，目前，基于DEM確定最佳集水面積閾值的方法主要包括流域?qū)挾确ā⒑泳W(wǎng)密度法、分形法和二階導數(shù)法等［12-17］，雖然這些方法可以較準確地提取最佳河網(wǎng)閾值，但是仍然存在一定的主觀性及不便捷性.常直楊等［18］通過均值變點法獲取了青藏高原東緣雜谷腦河流域河網(wǎng)密度值變化的拐點，并計算了該流域最佳集水面積閾值，證明均值變點法可以較為準確地確定河網(wǎng)的提取閾值，所提取河網(wǎng)與實際河網(wǎng)相比誤差小、可信度高.本次研究利用均值變點法確定最佳集水面積閾值，在此基礎(chǔ)上提取岷江上游流域地表水系網(wǎng)絡(luò)，分析岷江干流東西兩側(cè)各亞流域盆地的發(fā)育階段，從而為深入研究青藏高原東緣隆升歷史提供地貌依據(jù).

1 研究區(qū)概況

發(fā)源于岷山南麓郎架嶺的岷江是長江水量最大的一級支流，其干流河道自北向南穿行于崇山峻嶺之間，進入都江堰灌區(qū)后，都江堰引水樞紐魚嘴將干流河道分為內(nèi)外兩江，兩江在彭山江口鎮(zhèn)復(fù)合后南流至樂山市，右納大渡河和青衣江，至宜賓匯入長江.岷江源頭至都江堰魚嘴段為岷江上游，全長約340 km，流域面積約為21 580 km2.岷江上游地處北西西向的松潘—甘孜地槽褶皺帶、近東西向的西秦嶺地槽褶皺帶和北東向的龍門山斷裂帶的交匯部位，焚風效應(yīng)顯著，使岷江上游干流沿線形成典型干熱河谷，河谷兩岸植被稀疏，表土裸露.在岷江上游地區(qū)，主要發(fā)育2個活動斷層體系，北部為岷山構(gòu)造帶西緣的岷江斷裂，南部為龍門山?jīng)_斷帶的汶川—茂縣斷裂.岷江上游干流河道沿岷江斷裂和汶川—茂縣斷裂南流，橫切龍門山構(gòu)造帶，向東注入成都盆地.

2 研究方法

對岷江上游流域地貌發(fā)育階段進行研究，運用均值變點法確定水系網(wǎng)絡(luò)提取的最佳集水面積閾值，以此為基礎(chǔ)生成岷江上游水系網(wǎng)絡(luò)；以各支流和岷江干流河道的交匯點作為各亞流域盆地的匯水點，利用起伏比法對所提取岷江干流東西兩側(cè)各亞流域盆地的發(fā)育階段進行劃分，進而探討岷江上游流域地貌發(fā)育階段的特征.

2.1 均值變點法目前，ArcGIS水文分析模塊是最常用的流域水系網(wǎng)絡(luò)提取工具.其基本原理為：在對DEM圖進行填洼處理生成無洼地圖件的基礎(chǔ)上，運用D8算法確定水流離開每一個柵格后的流向，再根據(jù)水流方向數(shù)據(jù)統(tǒng)計區(qū)域匯流累積量，基于匯流累積量設(shè)置不同集水面積閾值，以所設(shè)置閾值的數(shù)值為基準，經(jīng)過柵格運算和矢量化處理，生成最終水系網(wǎng)絡(luò).值得注意的是，集水面積閾值大小與所提取河網(wǎng)疏密程度呈反比，即隨著閾值由小變大，河網(wǎng)結(jié)構(gòu)由密集趨于稀疏，究其原因，主要是隨著集水面積閾值的增大，匯流累積柵格中超過集水面積閾值的柵格變少，河道起點位置向地勢平坦處“退縮”［18］.因而，集水面積閾值的確定是水系網(wǎng)絡(luò)準確提取的關(guān)鍵因素.

均值變點分析法是確定集水面積閾值較為成熟的方法.如前所述，隨著集水面積閾值的增大，提取的水系將變得稀疏，河網(wǎng)密度會逐漸降低.河網(wǎng)密度曲線存在一個由陡變緩的拐點，拐點所對應(yīng)的閾值即為研究區(qū)域最佳集水面積閾值.為確定這一拐點［18-21］，設(shè)樣本數(shù)為N，樣本以Xi點為界分為2段，分別計算每段樣本的算術(shù)平均值ˉxi1和ˉxi2及全體樣本均值ˉx，有

其中，Si為樣本分段后的統(tǒng)計量，S為原始樣本統(tǒng)計量，N為樣本數(shù)，下腳標t1=1，2，…，i-1；t2=i，i+1，…，N.S和Si的最大差值對應(yīng)的點即為河網(wǎng)密度曲線由陡變緩的拐點，其對應(yīng)的閥值即為最佳集水面積閾值.

2.2 起伏比法起伏比法能夠?qū)α饔虻孛舶l(fā)育階段進行定量描述，在地貌學研究中得到了廣泛的應(yīng)用［22-23］.

將都江堰市定義為整個岷江上游流域盆地的匯水點，在岷江干流東西兩側(cè)系統(tǒng)提取研究區(qū)域所有亞流域盆地.在此基礎(chǔ)上，通過起伏比法分析各亞流域盆地發(fā)育階段的期次.

式中，HI為亞流域盆地的起伏比，Hmean為亞流域盆地平均高程（m），Hmax為亞流域盆地最大高程（m），Hmin為亞流域盆地最小高程（m）.按照戴維斯地貌旋回理論，流域地貌演化時間越長，侵蝕程度越高，起伏比越低.當HI＜0.4時，地貌處于均衡發(fā)育階段，為地貌發(fā)育的老年期；反之，當HI＞0.6時，地貌發(fā)育階段為幼年期；HI值介于兩者之間時，則為地貌發(fā)育的壯年期階段［11］.為了更詳細地劃分研究區(qū)各亞流域盆地的地貌發(fā)育期次，將起伏比為0.4～0.5的壯年期細化為晚壯年期，0.5～0.6的壯年期細化為早壯年期.

3 結(jié)果與討論

3.1 水系網(wǎng)絡(luò)生成本次研究采用的集水面積閾值為柵格的數(shù)目，試驗設(shè)置的柵格數(shù)目分別為1 000、2 000、3 000、4 000、…17 000、18 000、19 000、20 000.從圖1可以看出，隨著集水面積閾值（柵格數(shù)目）由小變大，岷江上游流域河網(wǎng)趨于稀疏，其密度值逐漸減小，集水面積閾值與河網(wǎng)密度擬合方程的判定系數(shù)R2=0.952 8，擬合效果較好.進一步分析可知，隨著閾值增加，河網(wǎng)密度在最初階段迅速減小，當達到某一閾值時，河網(wǎng)密度的減速變緩而趨于相對平穩(wěn)，減速變緩所對應(yīng)閾值即為研究區(qū)域最佳集水面積閾值.

由圖1可知，河網(wǎng)密度曲線必然存在一個由陡變緩的拐點，且該點唯一.由（1）式計算出S和Si的差值，最大差值對應(yīng)的點即為河網(wǎng)密度曲線由陡變緩的拐點，其對應(yīng)的最佳集流閾值柵格數(shù)目為4 000（圖2），據(jù)此通過ArcGIS水文分析模塊所提取的水系網(wǎng)絡(luò)如圖3.

圖1 河網(wǎng)密度與集水面積閾值（柵格數(shù)目）的擬合曲線Fig.1 The fitting curve of drainage density and threshold value of accumulation area（grid number）

圖2 原始樣本統(tǒng)計量S和分段樣本統(tǒng)計量Si差值變化曲線Fig.2 Change curve of the difference value between original samples statistics S and segmented samples statistics Si

3.2 流域發(fā)育階段劃分將岷江上游干、支流交匯點作為各亞流域盆地匯水點，在研究區(qū)域系統(tǒng)提取136個亞流域盆地，其中，1＃～66＃盆地位于岷江干流東側(cè)，67＃～136＃盆地位于岷江干流西側(cè).利用起伏比法計算各亞流域盆地發(fā)育期次，計算結(jié)果見表1和圖4.

表1 岷江上游各發(fā)育期次亞流域盆地面積、周長占比Tab.1 The proportion of area and perimeter of Minjiang drainge sub-basins in different development periods

從圖3和圖4可以較為直觀地看出，岷江干流河道兩側(cè)亞流域盆地呈不對稱分布特征，干流西岸支流數(shù)量多、流程長，東岸支流數(shù)量少、流程短，面積較大的亞流域盆地主要分布于岷江干流河道西岸的川西高原之上.究其原因，主要是川西高原廣闊、平坦的地表為各支流水系提供了足夠的發(fā)展演化空間.另外，由圖4可知，岷江干流河道主要沿呈近南北走向的岷江斷裂和呈北東—南西走向的汶川—茂縣斷裂發(fā)育，河谷的位置與斷裂的位置基本重疊，表明岷江上游干流水系的流向主要受斷裂帶及與之大致平行的巖層破碎帶所控制.斷裂的逆沖及右旋走滑活動易造成巖石破碎，從而在巖石表面形成抗風化軟弱帶，河流易于沿薄弱的抗風化軟弱帶下蝕，使河流的平面展布方向與活動斷裂的走向一致.

圖3 岷江上游流域水系網(wǎng)絡(luò)（集流閥值4 000）Fig.3 Drainage networks in the upper reaches basin of Minjiang river

從地貌發(fā)育期次定量統(tǒng)計數(shù)據(jù)看（表1），在岷江上游流域，處于早壯年期發(fā)育階段的亞流域盆地表現(xiàn)為數(shù)量多（70個）、總面積大（總面積為17 844.74 km2，占研究區(qū)面積的82.69%）、總周長長（總周長為3 411.70 km，占研究區(qū)亞流域盆地總周長的63.20%）等特征.從各發(fā)育期次亞流域盆地空間分布狀況看，岷江上游流域各亞流域盆地呈現(xiàn)出較明顯的區(qū)域分異特征.岷江斷裂和汶川—茂縣斷裂所圍合川西高原地區(qū)的大型亞流域盆地（小姓溝盆地、松坪溝盆地、黑水河盆地、雜谷腦河盆地和漁子溪盆地等）均處于早壯年期發(fā)育階段（這亦是該發(fā)育階段亞流域盆地在總面積、總周長均占優(yōu)勢的主要原因）.總體來看，斷裂圍合區(qū)亞流域盆地發(fā)育期次比斷裂外側(cè)（岷山斷塊和龍門山?jīng)_斷帶地區(qū)）亞流域盆地發(fā)育期次相對年輕.進一步分析可知，岷江斷裂和汶川—茂縣斷裂均為逆斷層，其斷層面傾向分別為正西向和北西向，即斷裂所圍合川西高原地區(qū)在上述兩條斷裂帶上均為上盤塊體，斷裂活動驅(qū)動的逆沖作用導致斷裂上盤區(qū)域沿斷層面垂直隆升，且兩條斷裂帶自第四紀以來活動均較為強烈；因而，位于川西高原之上的亞流域盆地地表抬升迅速、水流侵蝕劇烈、水系不斷擴展分支，亞流域盆地處于不均衡發(fā)育狀態(tài).反之，由于逆斷層下盤區(qū)域地表相對穩(wěn)定，水流侵蝕程度和支流水系擴展程度均不及上盤區(qū)域，亞流域盆地發(fā)育趨于平穩(wěn)化和老齡化；因而，處于老年期發(fā)育階段的亞流域盆地幾乎都位于斷裂帶外側(cè)的斷層下盤塊體之上（圖4）.

圖4 岷江上游流域亞流域盆地發(fā)育期次Fig.4 Geomorphologic development stage of the upper reaches of Minjiang drainge basin

4 結(jié)論

基于DEM數(shù)據(jù)，以ArcGIS水文分析模塊為技術(shù)平臺，系統(tǒng)提取岷江上游流域盆地地表水系網(wǎng)絡(luò)，并對岷江干流東、西兩側(cè)亞流域盆地的周長、面積、演化階段進行統(tǒng)計，得到如下結(jié)論：

1）隨著集水面積閾值由小到大變化，所提取的河網(wǎng)結(jié)構(gòu)由復(fù)雜趨于稀疏，通過均值變點法獲取河網(wǎng)密度值變化拐點，計算出岷江上游流域最佳集流閾值的柵格數(shù)目為4 000，以此為基礎(chǔ)生成研究區(qū)域水系網(wǎng)絡(luò).

2）將都江堰市定義為整個岷江上游流域盆地的匯水點，將各支流與岷江干流河道的交匯點作為各亞流域盆地的匯水點，在研究區(qū)域系統(tǒng)提取136個亞流域盆地，與干流東岸亞流域盆地相比，岷江干流西岸亞流域盆地表現(xiàn)為盆地面積大、周長和水系總長度長等特征.究其原因，主要是干流西岸的川西高原為各支流水系發(fā)育提供了廣闊演化空間.

3）由于岷江斷裂和汶川—茂縣斷裂所圍合川西高原地區(qū)的大型亞流域盆地（小姓溝盆地、松坪溝盆地、黑水河盆地、雜谷腦河盆地、漁子溪盆地等）主要位于斷裂上盤塊體區(qū)域，上盤區(qū)域的逆沖活動導致位于其上的亞流域盆地處于不均衡發(fā)育狀態(tài).因此，斷裂所圍合川西高原區(qū)域的亞流域盆地發(fā)育期次總體上比斷裂外側(cè)岷山斷塊和龍門山?jīng)_斷帶地區(qū)的亞流域盆地發(fā)育期次相對年輕.