三峽水庫運行后宜昌-城陵磯河段沖刷重心下移與時空演變

王華琳，鄭 珊，談廣鳴，李凌云

（1. 武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072；2. 長江科學院 河流研究所，湖北 武漢 430012）

1 研究背景

大壩修建與大型水庫運用會對下游河道的沖淤演變產(chǎn)生重要影響，如我國的三峽水庫［1-4］、丹江口水庫［5］、小浪底水庫［6-7］、三門峽水庫［8-10］，國外的Parker 水庫［11］、Wloclawek 水庫［12］等。壩下河道演變的響應表現(xiàn)在多個方面，在垂向和縱向上，河床呈現(xiàn)不同程度的沖刷下切［13-14］，河道縱比降減?。?3］；在橫向上，河道可向?qū)挏\或窄深方向發(fā)展［11，15］，如波蘭的Wloclawek大壩下游河道發(fā)生沖刷，河床粗化后沖刷受到抑制，河岸侵蝕加劇，河道展寬［11］，而美國的特里尼蒂河（Trinity River）壩下游河道因水庫削峰、岸坡植被生長導致河寬減小20%～60%［15］。在平面形態(tài)上，河道內(nèi)的洲灘格局及河型也可能改變，例如，國內(nèi)外研究發(fā)現(xiàn)壩下河道沖刷使得洲灘面積和數(shù)量減少［16-18］，壩下河道可能由辮狀向彎曲型轉(zhuǎn)化［19］。

三峽水庫運行后，長江中下游河道普遍發(fā)生沖刷，其中宜昌至城陵磯河道（以下簡稱宜城河道）沖刷最為明顯。以往研究顯示［20］，宜昌至城陵磯、城陵磯至漢口、漢口至湖口2002—2016年累計沖刷量分別為11 億、4.6 億和5.2 億m3，單位河長沖刷強度分別為18.7 萬、12.6 萬和12.2 萬m3/（km·a），可見宜城河段沖刷強度最大。該河段的演變特征有：沖刷主要集中在中枯水河槽［13，21-22］，深泓普遍下切［23］，床沙粗化且沿程粗化程度趨于減弱［24］，粗顆粒泥沙在監(jiān)利站附近恢復至建壩前水平，而細顆粒懸沙由于河床補給不足，恢復距離較長［22］，不同河型河段的演變發(fā)生較大變化，如彎曲河段多發(fā)生崩岸與“撇彎切灘”，斷面由單槽向“W”型雙槽轉(zhuǎn)化［25］，分汊河段發(fā)生支汊沖刷發(fā)展等演變現(xiàn)象［26］。

以往研究表明沖刷具有向壩下游傳播、發(fā)展及強度減弱的變化特點［27-28］。例如，董炳江等［20］對2002—2016年宜昌至湖口河段沖淤研究發(fā)現(xiàn)，2012年后城陵磯以下河道河床沖刷明顯，認為壩下沖刷已發(fā)展至城陵磯以下；許全喜等［27］研究發(fā)現(xiàn)三峽壩下游粗沙補給帶逐漸由宜昌至枝城河段下移至沙市至監(jiān)利河段。

本研究關(guān)注沖刷強度最大的河段（稱為沖刷重心），其位置隨時間的變化對壩下河道的穩(wěn)定與防洪安全十分重要。然而，沖刷重心的位置與壩下沖刷延伸所達到的位置有所不同，沖刷重心的位置可能向下游遷移，也可能在上、下游河道之間擺動，由于上荊江河道沖刷強度較大，已有研究認為沖刷重心或主要沖刷帶目前仍位于上荊江［37］。羅方冰等［28］對比宜昌至陳家灣（壩下130 km長河段）不同子河段在2002—2014年沖淤量隨時間的變化，認為三峽壩下游主沖刷帶的發(fā)展速率約8 km/a，然而，其研究河段較短，對主沖刷帶位置的界定較為模糊，亟需引入定量方法，研究三峽壩下游河道沖刷重心的遷移規(guī)律。

由于壩下不同河段的水沙及邊界條件不同，河道時空沖淤變化十分復雜，缺乏壩下游沖刷重心變化規(guī)律的量化研究。本文以宜城河道為研究對象，基于2002—2018年的水沙與斷面實測資料，在分析河道形態(tài)變化的基礎(chǔ)上，構(gòu)建河道沖淤變化速率等級的時空矩陣，引入機器學習中kmeans++聚類方法［29］，研究壩下游河道沖淤演變及沖刷重心的遷移規(guī)律，為壩下游河道演變的趨勢預測與科學管理提供參考。

2 研究區(qū)域

如圖1 所示，宜城河段位于長江中游，由長約60 km 的宜昌至枝城河段（簡稱宜枝河段）和長約340 km的荊江河段組成，其中荊江又常分為上荊江（枝城-藕池口）和下荊江（藕池口-城陵磯）。以上荊江楊家垴為界，宜昌至楊家垴為砂卵石河床河段，楊家垴至城陵磯為沙質(zhì)河床河段［13］。宜枝河段右岸有清江入?yún)R，荊江河段右岸有松滋口、太平口、藕池口和調(diào)弦口（已建閘控制）分流向洞庭湖，左岸有沮漳河入?yún)R，洞庭湖出流在城陵磯附近匯入長江。

圖1 研究區(qū)域

三峽水庫運用后宜昌、沙市和監(jiān)利三站（圖1）的年均徑流量減少約4%～6%，年均來沙量分別比建庫前（1991—2002年）減少91%、85%和79%（圖2），來沙量減幅沿程減弱，反映了河床沖刷沿程補給泥沙的作用。此外，宜昌站床沙有所粗化，而監(jiān)利站床沙在建庫后波動變化，趨勢尚不明顯（圖3）。

圖2 水沙條件

圖3 床沙級配變化（資料來自文獻［27，31，36］）

綜合考慮三峽水庫的不同運用方式和三峽上游梯級水庫投入運用的時間，將研究時段分為3段［31］：Ⅰ）2003—2007年，三峽水庫蓄水初期；Ⅱ）2008—2012年，三峽水庫175 m 蓄水期；Ⅲ）2013—2018年，三峽上游梯級水庫運用后。3個時段宜昌站對應的年均徑流量分別為3936億、4020億和4191 億m3，與時段Ⅰ相比，時段Ⅱ和Ⅲ的年均徑流量分別增加2.1%和6.5%；年均來沙量分別為0.67億、0.30億和0.11億t，與時段Ⅰ相比，時段Ⅱ和Ⅲ的年均來沙量分別減少55.2%和83.6%。

空間上，將研究河段劃分為32個子河段（圖1）。根據(jù)河道平面特征和162個觀測資料較全的斷面（宜34 至荊181 斷面）的分布劃分子河段，考慮河段為單一或分汊河型，或為順直、微彎或彎曲河型，可分為6種不同河型（單一順直、單一微彎、單一彎曲、順直分汊、微彎分汊和彎曲分汊），分汊型子河段包含完整的江心灘，彎曲型子河段包含彎頂段及部分上、下游銜接段，在考慮河型的同時，盡量避免各子河段的河長差別太大。河道曲率大于1.2為彎曲型，小于1.2為順直微彎型［30］，為進一步區(qū)分順直與微彎河型，以1.05為界，認為曲率介于1.05 ～1.2為微彎河型，小于1.05為順直河型，需要說明的是，本文采用1.05為界僅為區(qū)分不同子河段的相對彎曲程度，并不一定適用于其它河道。對多汊河段，采用主汊計算河道曲率。劃分子河段時不可避免存在主觀因素的影響，各子河段沖淤量與劃分方法有關(guān)，但河道總體的時空沖淤規(guī)律不會受到影響。

3 研究方法

3.1 河床形態(tài)特征量計算基于實測資料，計算各斷面的形態(tài)特征量，包括深泓高程、深泓橫向擺動速率、平灘水深和平灘面積。平灘面積指平灘高程以下的橫斷面面積，本文平灘高程根據(jù)斷面灘唇位置選?。?2］，對于無明顯灘唇的斷面，根據(jù)上、下游有灘唇可確定平灘高程的斷面資料，擬合平灘高程與河長的關(guān)系式，插值得到無明顯灘唇的斷面處的平灘高程。根據(jù)平灘高程確定平灘河寬W和平灘面積A，計算得到平灘水深H。

采用下式計算各斷面的年際深泓擺動速率Δx（m/a）：

式中xi和xi+1分別為第i和第i+1年斷面深泓點的起點距，m。

根據(jù)各斷面的河床形態(tài)特征量，采用式（2）計算河段平均形態(tài)特征量：

式中：G泛指某一斷面的河床演變特征量，為河段尺度的特征量值，可將各斷面的深泓高程的變化、深泓橫向擺動速率Δx、平灘面積和平灘水深的變化作為G分別代入式（2）計算其河段平均值；Gj和Gj+1分別為第j和第j+1個斷面的形態(tài)特征量；Lj，j+1為第j與第j+1個斷面之間的河長，km；L1，M為河段長度，km；M為河段的斷面數(shù)。

根據(jù)各斷面平灘面積，采用式（3）計算河段沖淤量V（m3）：

式中ΔAj和ΔAj+1分別為第j和第j+1個斷面平灘面積的變化，m2。

3.2 沖淤等級劃分根據(jù)2003—2018年32個子河段的平均平灘水深和平灘面積的變化（共512個數(shù)據(jù)），對平灘面積和平灘水深的變化速率進行分級。除去平灘面積和平灘水深變化最小的前5%的數(shù)據(jù)（該數(shù)據(jù)歸為0 等級）之后，對剩余數(shù)據(jù)按出現(xiàn)頻率進行等分（各等級出現(xiàn)頻率各占23.75%），得到±1、±2、±3和±4四類等級，負數(shù)等級代表沖刷，即平灘面積或水深增大，正數(shù)代表淤積，即平灘面積或水深減小，各等級對應的平灘面積和平灘水深變化的絕對值范圍如圖4所示，需要注意的是，由于平灘面積和水深的變化幅度較大，±3和±4等級對應的平灘面積和水深變化的跨度較大，本文按頻率劃分等級可保證-4等級具有一定的樣本點數(shù)量，以便于對沖刷強度最大的-4等級對應的樣本點進行聚類分析。

圖4 子河段平灘面積和平灘水深變化的等級劃分

3.3 沖淤等級的時空聚類根據(jù)上述方法得到各子河段每年平灘水深和平灘面積變化的等級，選取沖刷最強烈的-4等級，采用空間聚類中的kmeans++方法［29］，研究沖刷重心（即平灘水深沖刷最強和平灘面積增大最多的子河段區(qū)域）隨時空的變化規(guī)律。

kmeans++方法采用距離作為相似性的評價指標，即兩個對象距離越近，相似度越大，由此得到緊湊且獨立的簇。對圖5所示樣本點，以k=2為例說明采用簡化kmeans++算法的聚類過程［29］：（a）在N個樣本點中隨機選擇第一個初始中心點A；（b）計算其余N-1點到已有中心點的最短距離D（x），計算各點被選為下一個中心的概率P：

圖5 簡化kmeans++算法示意

取概率最大的點為第二個中心點，離上一個中心值最遠的點B具有最大機會作為下一次選取的中心值；（c）計算各點到兩個聚類中心A和B點的距離，按照距離就近原則分類，圓點屬于一類，其中心點為A點，三角形點屬于另一類，其中心點為B點；（d）按現(xiàn)有分類計算各類重心位置，得到新的聚類重心，重復（c）得到新類簇；迭代至重心位置不再變化。kmeans++方法在步驟（b）時采用輪盤賭算法確定下一個中心點位置［29］，本文對這一方法進行了簡化，取概率P最大的點為中心點，這與kmeans++方法的思想核心即離上一個中心值較遠的點更有機會作為下一次選取的中心值基本一致。

4 結(jié)果分析

4.1 河道沖淤演變

4.1.1 縱向變化 2003—2018年宜城河段河床普遍沖刷下降，河床沖深在空間呈不均衡分布（圖6）。平灘水深與深泓高程的變化趨勢基本一致（圖7），在3個研究時段內(nèi)，全河段年均深泓調(diào)整速率分別為-0.41、-0.17和-0.20 m/a，平灘水深變化速率分別為0.20、0.12和0.14 m/a，沖刷的子河段河長分別占全河段總長的93%、69%和88%，由此可見，蓄水初期（時段Ⅰ）河床沖刷強度最大，175 m蓄水運用后（時段Ⅱ）沖刷減弱，三峽上游梯級水庫運用后（時段Ⅲ）沖刷增強，但強度弱于蓄水初期。時段Ⅱ沖刷減弱與河床粗化有關(guān)，上游梯級水庫運用后沖刷增強與來沙量較時段Ⅱ大大減少有關(guān)。

圖6 深泓縱剖面及其變化

圖7 河床沖淤變化

4.1.2 深泓橫向擺動 由于分汊河段包含多個汊道，較難采用統(tǒng)一方法對比單一和分汊型河道深泓的橫向擺動速率，因此，本節(jié)僅對單一型河段的深泓擺動進行分析。需要注意的是，子河段8（上荊江董2—荊13）受河道心灘影響，以往研究將其歸為單一［17］或分汊型［18］河道，本節(jié)暫不對這一子河段進行分析。

宜枝河段和下荊江以單一河型為主，宜枝河段受河道邊界條件控制橫向擺動較弱，下荊江河道深泓擺動強度較大（圖8）。除個別子河段外（如子河段18和28），2013—2018年深泓擺動強度明顯弱于前期兩個時段，說明單一型河道的橫向擺動有所減弱。圖9 顯示了不同典型形態(tài)斷面的形態(tài)變化，位于第4個子河段的宜63斷面主槽呈V型，左側(cè)邊灘沖刷導致深泓橫向擺動，位于第28個子河段的荊149斷面趨于W型，主流在心灘兩側(cè)頻繁擺動，橫向擺動速率較大。

圖8 深泓橫向擺動速率

圖9 深泓擺動的典型斷面

4.1.3 河道沖淤 采用式（3）計算得2003—2018年宜城河道累計沖刷量約11.9 億m3。對比本文方法與其他學者計算不同時段河道沖淤量，考慮所采用的方法有所不同，但相對誤差較?。ū?），認為本文方法計算的沖刷量結(jié)果基本可信。

表1 本文方法與以往研究計算沖淤量結(jié)果對比

宜城河道在3個時段內(nèi)年均沖刷量分別約0.86億、0.59億和0.78億m3，再次證明河道在梯級水庫運行后沖刷速率增大，但沖刷強度仍弱于三峽水庫運行初期。2003—2018年宜枝、上荊江和下荊江河段單位河長沖刷速率分別為17.4 萬、23.6 萬和16.8 萬m3/（km·a），上荊江的沖刷強度最大，下荊江最?。▓D10（a））。

3個時段內(nèi)沖刷強度最大的峰值分別出現(xiàn)在第5、第7和第13子河段（見圖10（b）中箭頭），即壩下游約85、106和178 km，相應沖刷強度分別為70萬、65萬和100萬m3/（km·a）。沖刷強度最大的峰值隨時間不斷下移，但第3個時段沖刷峰值最大，可能與梯級水庫運用后河道沖刷增強及其它因素（如河道邊界條件、采砂等［27］）有關(guān)。沖刷強度最大峰值的位置與平灘水深變化最大峰值的位置（如圖7（b））一致，說明河道垂向沖深與斷面面積增大基本同步，這與以往研究認為宜城河道沖刷總體為垂向沖深、橫向展寬居次要地位的結(jié)論一致［11，30］。

4.2 沖淤等級時空分布與沖刷重心變化規(guī)律

4.2.1 沖淤等級的時空矩陣 2003—2018年32個子河段平灘面積變化速率等級的時空矩陣如圖11所示，圖中藍色代表沖刷（平灘面積增大），紅褐色代表淤積（平灘面積減?。?，顏色越深代表沖淤越劇烈。平灘水深變化速率等級與平灘面積變化速率等級的時空分布相似。

計算2003—2018年32 個子河段的平灘面積和平灘水深變化等級的均值和方差，結(jié)果如圖12 所示，由圖可知，除砂卵石河床的宜枝段以外，平灘面積和水沙變化等級的均值和方差具有較顯著的向下游增大的趨勢，變化等級的均值沿程增大說明越往下游河段沖刷越弱或越易發(fā)生淤積，方差沿程增大說明下游河段沖淤等級差異性越大，或越易發(fā)生沖淤交替。圖11沖淤等級的時空分布矩陣也顯示下荊江子河段相對上荊江更易出現(xiàn)沖淤波動。需要說明的是，由于宜枝河段屬于山區(qū)到平原河流過渡段，河床抗沖性較強，三峽水庫運用后沖刷強度不大，甚至在一些年份發(fā)生淤積，因此沖淤等級的均值也較大。

圖11 平灘面積變化速率等級的時空分布與沖刷重心

圖12 2003—2018年各子河段沖淤等級和方差

4.2.2 沖刷重心的時空分布與遷移 采用kmeans++算法對圖11中沖刷等級為-4等級（沖刷速率最大）的樣本點進行聚類分析。首先需確定聚類中心個數(shù)k，k越大，用聚類中心代表總樣本時信息量的丟失程度越小，但當k值太大（如等于樣本點個數(shù)）時，則達不到分類效果，因此，需針對不同樣本點個數(shù)選擇合適大小的k值。本文對平灘面積變化的-4等級樣本點，計算k由1逐漸增大至樣本點個數(shù)對應的畸變函數(shù)（即式（4）右側(cè)分母部分），隨著k增大，畸變函數(shù)值逐漸減小，當k達到13時，畸變函數(shù)值趨近平穩(wěn)，因此本文取k=13。此外，kmeans++算法隨機給出初始聚類中心，本文通過6次隨機聚類得到78個聚類中心，其位置較鄰近，甚至有重疊，說明初始中心值的隨機選取對聚類結(jié)果的影響不大。需要說明的是，本文采用的空間聚類算法得到的沖刷重心位置，雖然在個別年份這個位置可能發(fā)生淤積，但其代表數(shù)年和數(shù)個子河段在時空平均意義上的沖刷最強區(qū)域的中心，應從較長的時空尺度上理解沖刷重心位置的遷移變化。根據(jù)平灘面積變化最大的-4 等級的沖刷重心的聚類結(jié)果，可以劃分如下4個時空沖淤區(qū)域（圖11）。

Ⅰ區(qū)：相對穩(wěn)定區(qū)，紅色虛線以上、矩陣右上角的三角形區(qū)域，在這一區(qū)域內(nèi)河道沖淤幅度較小，河道相對穩(wěn)定，并且沖淤幅度較小的河段或壩下相對穩(wěn)定區(qū)的河長隨時間不斷增加，至2018年壩下游約140 km長河段均處于壩下游相對穩(wěn)定區(qū)。

Ⅱ區(qū)：沖刷重心下移區(qū)，紅色與黑色虛線內(nèi)的傾斜帶狀區(qū)域，在這一區(qū)域內(nèi)分布著多個沖刷重心聚類點，這些沖刷重心聚類點具有隨時間向下游遷移的趨勢。圖11紅色矩形框顯示4.1節(jié)中3個時段沖刷峰值所在子河段（5、7和13）的位置，與沖刷重心的聚類結(jié)果基本吻合。

Ⅲ區(qū)：弱沖弱淤區(qū)，黑色虛線與黃色虛線內(nèi)的傾斜帶狀區(qū)域，這一區(qū)域內(nèi)河道沖淤幅度相對較小，約2010年前河床以弱淤為主，之后以弱沖為主，但沖刷強度明顯弱于沖刷重心下移區(qū)。

Ⅳ區(qū)：沖淤交替區(qū)，沖刷和淤積交替出現(xiàn)的概率較大，沖刷重心上、下游波動明顯。該區(qū)域主要集中于下荊江河段（第20—28子河段），如圖13所示，典型子河段（荊108—荊120、荊123—荊133和荊147—荊166）沖淤交替頻繁，但沖淤幅度隨時間不斷降低（圖13）。需要說明的是，與其它3個區(qū)相同，Ⅳ區(qū)內(nèi)河道沖刷仍占主導，河道整體表現(xiàn)為沖刷，4個區(qū)域并沒有嚴格的界限，圖11所示分界線僅為示意。平灘水深變化最大的-4等級的聚類結(jié)果與平灘面積的聚類結(jié)果基本一致。

圖13 沖淤交替區(qū)代表子河段年沖淤量變化

選取平灘面積和平灘水深聚類所得Ⅱ區(qū)中的沖刷重心，得到?jīng)_刷重心隨時空的分布，如圖14所示。蓄水初期沖刷重心下移速率較快，2008年后逐漸減慢，上游梯級水庫運用后，沖刷重心下移速率明顯加快，至2014年后又有所減慢。采用直線擬合沖刷重心在2003—2018年的遷移變化（對平灘面積R=0.92，對平灘水深R=0.96），擬合線的斜率可作為三峽水庫運用以來沖刷重心的平均下移速率，約為7.5 km/a。

圖14 沖刷重心隨時間下移

上述沖刷重心的遷移特點與以往研究結(jié)論基本一致。例如，本文得到2003—2018年沖刷重心的平均下移速率約為7.5 km/a，與羅方冰等［28］得到三峽壩下游主沖刷帶的發(fā)展速率約8 km/a 基本一致，但本文沖刷重心下移速率是基于量化的數(shù)據(jù)處理與分析得到的，并且顯示了不同時段沖刷重心下移的快慢變化。本文結(jié)果顯示沖刷重心在2010年左右下移至宜枝河段末端（圖11 與圖14），以往研究顯示三峽水庫運用后至2010年左右宜枝河段床沙粗化明顯、之后變化不大，河道沖刷減弱［14，34］，這一變化規(guī)律與沖刷重心下移至宜枝河段以下有關(guān)。此外，圖11 和圖14 顯示2013年沖刷中心下移至宜昌下游約100 km 的位置，相對穩(wěn)定區(qū)基本涵蓋松滋口及其上游河段范圍，這與岳紅艷等［35］的研究結(jié)論即松滋口以上河段在2013年以后沖刷速率明顯減緩、河道平面形態(tài)及河勢基本穩(wěn)定一致。

圖15分別給出了4個分區(qū)對應的典型子河段累計沖刷量的變化過程。子河段1（宜34—宜45）位于相對穩(wěn)定區(qū)（Ⅰ區(qū)），累計沖刷量隨時間變化不大；子河段4（宜60—宜69）于2010年左右由沖刷重心下移區(qū)（Ⅱ區(qū)）進入相對穩(wěn)定區(qū)（Ⅰ區(qū)），沖刷速率由快變慢；子河段13（荊30—荊36）于2012年左右由弱沖弱淤區(qū)（Ⅲ區(qū)）進入沖刷重心下移區(qū)（Ⅱ區(qū)），沖刷速率由慢變快；子河段28（荊147—荊166）位于沖淤交替區(qū)（Ⅳ區(qū)），河道以沖為主，但累計沖刷量波動較明顯。

圖15 4個區(qū)典型子河段的累計沖刷量變化

5 結(jié)論

三峽水庫運用后，宜昌至城陵磯河道經(jīng)歷了蓄水初期沖刷較強（2003—2007）——正常蓄水后沖刷減弱（2008—2012）——上游梯級水庫運用后沖刷增強（2013—2018）三個階段；在3 個時段內(nèi)，平灘水深變化速率分別為0.20、0.12 和0.14 m/a，年均沖刷量分別約0.86 億、0.59 億和0.78 億m3。單一河型河段在2013年后橫向擺動有所減弱。建立了2003—2018年32 個子河段平灘面積和平灘水深沖淤變化等級的時空分布矩陣，發(fā)現(xiàn)沖淤等級的均值和方差沿程增大，表明往下游河道沖刷逐漸減弱、年際沖淤交替越頻繁。以順直微彎分汊型為主的上荊江河道普遍沖刷，而以彎曲和彎曲分汊為主的下荊江河道沖淤交替頻繁。針對平灘面積和平灘水深增速最大等級對應的子河段進行聚類分析，發(fā)現(xiàn)沖刷重心具有向下游遷移的變化特點，在三峽蓄水初期沖刷重心向下遷移較快，之后減緩，在上游梯級水庫運用后，沖刷重心下移速率明顯加快，之后再次減緩，2003—2018年沖刷重心平均下移速率約7.5 km/a；至2018年沖刷重心已下移至壩下游約140 km處。未來沖刷重心可能繼續(xù)沿壩下游河道遷移和發(fā)展，然而，其變化受到與來水沙條件和河道邊界調(diào)整的影響，具有較強的不確定性，考慮三峽水庫運行時間較短，后續(xù)需在更長時空尺度上深入研究沖刷重心的時空遷移規(guī)律與影響機理。