黃河源區水沙變化特征分析

董曉寧，田世民，張 麗

（1.華北水利水電大學，河南 鄭州450046；2.黃河水利科學研究院 河南省黃河水生態環境工程技術研究中心，河南 鄭州450003；3.黃河水利科學研究院 河南省黃河流域生態環境保護與修復重點實驗室，河南 鄭州450003）

受全球氣候變暖的影響，水循環過程發生了劇烈變化。黃河源區是黃河流域的重要產流區和水源涵養區，多年（1956—2017年）平均年徑流量約198.2億m3，占黃河年徑流量的34.1%，多年平均年輸沙量為0.1億t，水資源狀況對整個黃河流域有極大的影響，且其輸沙量對上游梯級水庫的壽命有極大影響。黃河源區地處青藏高原腹地，地形復雜，海拔落差較大，且流域內還分布有冰川、凍土、湖泊、濕地等，生態系統多樣且脆弱，人類活動痕跡較少，幾乎處于自然狀態，水文過程對氣候變化敏感程度較高，氣候變暖導致的水循環劇烈變化在這里得到放大，多種影響因素疊加導致水沙過程較為復雜。

由于黃河源區對西部地區及黃河流域的重要性，因此許多專家學者深入分析了其徑流演變機制。流域內水沙過程在時空變化上主要受降水影響，且降水轉化率偏低［1］，此觀點得到較多學者的認同［2-4］。除此之外，氣溫變化引起的蒸發變化［5-6］、冰川消退［7］和凍土的凍融過程［8］也對驅動水沙變化有一定的貢獻。地表改造和引水等活動也是水沙特征改變的一大驅動因素［9］。湖泊、濕地的調蓄作用對流域內部分地區的水沙過程也存在影響［10］。降水與水沙變化成正比，氣溫與水沙變化成反比［11］。有學者對黃河源區近幾十年徑流減少各影響因素作用占比進行分析，發現20世紀90年代到2005年左右徑流減少主要受人類活動影響，而2008年以后徑流減少則主要受氣候因素影響［12］。也有學者預測隨著冰川積雪的減少，融雪徑流將會減少，部分河流徑流量會出現由增轉減的“拐點”［13］，流域內多年凍土也會持續退化為季節凍土［14］。總之，黃河源區氣候、下墊面、人類活動等與水文過程之間存在著相當復雜的關系。

要探究黃河源區水沙變化規律，僅從氣候變化和人類活動等外在因素分析黃河源區徑流、泥沙變化特征還不夠全面，空間上流域內各區域的水沙變化特征及其對整個流域的水沙變化貢獻率也值得研究，藍永超等［15］認為在對黃河源區徑流變化進行研究時，不應忽略流域內不同區域徑流變化的差異。2019年9月習近平總書記提出要推進黃河流域生態保護和高質量發展，并對三江源、祁連山、秦嶺等重點區域生態保護建設提出要求，以三江源、祁連山、甘南黃河上游水源涵養區等為重點，推進實施一批重大生態保護修復和建設工程，提升水源涵養能力。筆者運用統計學分析方法分析黃河源區把口站唐乃亥站徑流泥沙演變趨勢及特征，并將黃河源區分為黃河沿以上、黃河沿—吉邁、吉邁—瑪曲、瑪曲—唐乃亥4個區間，從空間上分析徑流泥沙在各區間不同時間段水沙變化規律及區間產流產沙對唐乃亥站水沙變化的貢獻率，以期為黃河源區生態系統恢復和水源涵養能力提高的區域性規劃提供參考。

1 資料來源和研究方法

1.1 資料來源

選取黃河源區黃河沿、吉邁、瑪曲、唐乃亥4個水文站的徑流和泥沙數據，這4個水文站具有較長的數據系列，且資料少有丟失，可以代表黃河源區流域內的水沙變化。由于歷史原因，黃河沿站1968—1975年水文資料缺失，下游吉邁站與黃河沿站的水文資料具有較好的相關關系，因此使用吉邁站資料對黃河沿站缺失資料進行插補。徑流、泥沙資料來自水利部黃河水利委員會刊印的黃河流域水文年鑒。

1.2 研究方法

1.2.1 趨勢檢驗

趨勢檢驗使用Mann-Kendall非參數趨勢檢驗法，根據統計量Z值的大小來確定序列變化趨勢以及趨勢是否顯著。Z值大于0，序列呈上升趨勢；Z值小于0，序列呈下降趨勢。當Z的絕對值大于1.28、1.64、2.32時，分別表示趨勢通過了顯著性水平為90%、95%、99%的檢驗。Sen’s斜率估計是Sen于1968年提出并發展的一種非參數檢驗法，估計n個樣本中N對數據的趨勢斜率，其值可以反映數據趨勢的陡峭程度。

1.2.2 突變檢驗

綜合累積距平法和Mann-Kendall突變檢驗確定突變點，再用t檢驗確定突變點是否合理。距平是離散數據與其均值的偏差，累積距平是距平值的疊加。累積距平值的起伏變化可以用來判斷序列的長期變化趨勢以及持續時間，起伏拐點處則可能為序列突變點。Mann-Kendall突變檢驗是根據時間序列數據構造秩序列，根據秩序列計算統計量UFk和UBk，給定顯著性水平，統計量UFk和UBk在臨界值內的交點即為序列突變點，若UFk或UBk的值大于0，則表明序列呈上升趨勢，小于0則表明呈下降趨勢。當它們超過臨界直線時，表明上升或下降趨勢顯著。

1.2.3 周期檢驗

周期性檢驗使用Morlet小波分析方法，繪制小波系數實部等值線圖和小波方差圖，峰值對應的點即為周期，最大峰值所對應的周期即為主周期，以此確定水沙序列的變化周期。

1.2.4 趨勢持續性

用重標極差法（R/S分析法）計算Hurst指數、預測徑流泥沙在未來一段時間內的變化趨勢。當H=0.5時，表示時間序列不存在記憶性；當H＞0.5時，時間序列具有狀態持續性，未來趨勢與過去相似；當H＜0.5時，時間序列具有反持久性，未來趨勢與過去相反。Hurst指數分級見表1。

表1 Hurst指數分級

2 研究結果分析

2.1 水沙變化趨勢

唐乃亥站水沙變化曲線及累積距平曲線見圖1。1956—2017年黃河源區唐乃亥站多年平均流量為628.58 m3/s，多年平均輸沙率為369.94 kg/s。流量呈波動下降趨勢，下降速率為1.21 m3/（s·a），年均流量最大值和最小值出現在1989年和2002年，其值分別為1 035.57 m3/s和334.51 m3/s；輸沙率變化同樣呈波動下降趨勢，下降速率為1.62 kg/（s·a），年均輸沙率最大值與年均流量最大值出現年份相同，為1989年，其值為1 298.49 kg/s，而最小值則出現在2008年，其值為87.29 kg/s。累積距平曲線表明徑流和泥沙變化過程均經歷了1989年之前的上升階段和之后的下降階段。

圖1 唐乃亥站水沙變化曲線及累積距平曲線

計算得黃河源區水沙變化相關系數為0.85，說明徑流在很大程度上影響泥沙過程。用變差系數CV來表示徑流、泥沙的變異程度，泥沙過程的變差系數CV和極值比均大于徑流過程（見表2）。使用Mann-Kendall趨勢檢驗法對黃河源區徑流、泥沙過程進行分析，得到統計量Z分別為-1.02和-0.98，均未通過置信度為90%的顯著性檢驗，說明黃河源區徑流泥沙雖一直呈下降趨勢，但并未達到顯著性水平。

表2 徑流泥沙變異程度

2.2 水沙突變和周期性檢驗

2.2.1 突變檢驗

使用Mann-Kendall突變檢驗法（以下簡稱MK突變檢驗）對黃河源區1956—2017年徑流、泥沙過程進行分析，如圖2所示。結合前文所述累積距平曲線所展現出的特征，確定幾個較為可能的突變點，再對可能突變點進行t檢驗，根據統計量t值的大小決定突變的年份。此處給定顯著性水平α=0.05，MK突變檢驗上下臨界值為±1.96。

圖2 黃河源區水沙MK突變檢驗

由圖2可以看出，徑流的拐點在1989年、2008年，泥沙的拐點在1989年。MK突變檢驗中徑流過程UFk與UBk曲線的交點在1993年、2008年、2012年，泥沙過程UFk與UBk曲線的交點在2006年、2013年。將兩種方法所得的可能突變年份設為基準年，對其前后兩個序列做t檢驗，得到統計量t并判斷其顯著性，檢驗結果見表3。

表3 突變點顯著性t檢驗結果

由于本文研究的時間尺度基于年代際變化，因此上述可能突變點距離時間系列終止年份不足10 a的不再考慮，相鄰兩個突變年份相差不超過10 a的，取統計量t值高的突變點為突變年份。總結表3檢驗結果，1989年、1993年均檢測到徑流突變，但兩個突變年份僅相差5 a，取突變更顯著的1989年為徑流突變年份；1989年、2013年均檢測到泥沙突變，但2013年距資料系列終止年份僅有4 a，取1989年為泥沙突變年份。徑流、泥沙過程的突變年份均在1989年，區別在于徑流過程相對泥沙過程突變結果更為顯著［16］。

2.2.2 周期檢驗

對1956—2017年唐乃亥站徑流、泥沙過程進行小波分析，分別繪制小波系數實部等值線和小波方差，如圖3所示。徑流、泥沙過程存在周期性，從圖3可以看出，徑流和泥沙在周期性變化上存在相似性，徑流過程存在7、12、25 a的周期和45 a的準周期，泥沙過程存在5、12、25 a的周期和34、44 a的準周期。徑流過程的7 a周期主要發生在1966—1980年，12 a周期發生在1956—2000年，25 a周期發生在1980—2010年。泥沙過程的5 a周期主要發生在1980—1995年，12 a周期發生在1956—2000年，25 a周期主要發生在1980年以后，從圖3可以看出泥沙過程在1965—1980年也存在7 a的周期。結合小波系數實部等值線圖的周期密集性和小波方差圖的峰值大小來看，在5～7 a和12 a周期，泥沙過程震蕩強度大于徑流過程的，尤其在12 a周期的震蕩強度，而在25 a周期尺度上泥沙過程的震蕩強度則小于徑流過程的。

圖3 水沙小波分析

2.3 水沙變化的區間趨勢

黃河源區水文站點稀少，有50 a以上資料的站點僅有黃河沿、吉邁、瑪曲、唐乃亥4個，因此本文選用這4個站點1959—2017年的水沙資料來研究黃河源區的水沙變化特征，并將黃河源區分為黃河沿以上、黃河沿—吉邁、吉邁—瑪曲、瑪曲—唐乃亥4個區間，按區間上下游順序記為區間一、二、三、四，將相鄰上、下游兩站實測徑流泥沙數據之差作為區間產流產沙量。由表4可以看出，黃河源區的最大產流區為區間三，多年平均產流量占總產流量的50.80%；最大的產沙區為區間四，多年平均產沙量占總產沙量的62.53%。使用MK趨勢檢驗法對各區間水沙趨勢進行分析，由表5得出，4個區間除區間二徑流過程表現出不顯著上升趨勢外，其余3個區間都表現出下降趨勢，且主產流區區間三表現出置信度為99%的顯著下降趨勢；4個區間的泥沙量都表現出下降趨勢，區間一和區間四均通過置信度為90%的顯著水平檢驗，而區間二則通過置信度為99%的顯著水平檢驗。

表4 黃河源區水沙區間占比

表5 水沙MK趨勢檢驗

根據上述突變點選取標準，對表6可能突變點顯著性檢驗結果總結后得出，區間一徑流過程突變年份為1990年、2006年，泥沙過程突變年份為1994年；區間二徑流過程突變年份為1986年、2004年，泥沙過程突變年份為1985年；區間三徑流過程突變年份為1989年，泥沙過程突變年份為1994年；區間四徑流過程突變年份為1989年、2007年，泥沙過程突變年份為1989年。

2.4 不同時期的水沙區間變化

分析各區間水沙突變年份可以發現，黃河源區各區間水沙突變年份和變化趨勢具有很大的相似性，根據突變點年份再綜合水沙貢獻率區間占比，選擇1989年和2007年作為此次研究的“拐點”，由此將黃河源區水沙變化趨勢分為1959—1989年、1990—2007年和2008—2017年3個時段來研究各區間的水沙變化，按時間順序分別將3個時期定義為T1、T2、T3。各區間年均產流產沙量和不同時期年均水沙變化率分別見表7和表8。

表6 各區間水沙突變點顯著性t檢驗結果

從產流時期來看，所有區間都經歷了T2時期的徑流減少和T3時期的徑流回升，上游兩區間T3時期的年均產流量超過了T1時期，而下游兩區間T3時期徑流量雖有所回升，但并未達到T1時期的產流量。從T1—T3時期的年均產流量和變化率來看，區間三既是產流量減少最多又是變化率最大的區間。

表8 各區間不同時期年均水沙變化率 %

黃河源區產沙區域主要集中在區間三和區間四，其中約60%的泥沙來自區間四、約30%的泥沙來自區間三，而區間一和區間二的產沙量約占總產沙量的10%。與徑流變化相似，各區間產沙量也經歷了T2時期的減少，但在T3時期只有上游3個區間產沙量有所回升，區間四的產沙量及其占比均在持續下降。從年均產沙量來看，產沙量減少最多的是區間四，但從泥沙變化率來看，變化最大的是區間一。

2.5 水沙變化的區間貢獻率

將不同時期各區間對促進黃河源區水沙變化的作用占比定義為水沙變化區間貢獻率。黃河源區不同時期的年均水沙產量變化情況和區間貢獻率見表9。

根據表9，T1到T2時期黃河源區年均產流量減少48.14億m3，區間三和區間四對徑流減少的貢獻率較大，分別是45.60%和35.56%，區間一和區間二對徑流減少的貢獻率比較接近，均不到10%；T2到T3時期黃河源區年均產流量增加29.09億m3，區間四對徑流增加的貢獻率為55.92%，區間二和區間一對徑流增加的貢獻率分別為25.03%和17.98%，區間三對徑流增加的貢獻率則僅有1.07%；T1到T3時期黃河源區年均產流量減少19.05億m3，區間三對徑流減少的貢獻率達113.60%，區間一和區間二對徑流減少反而起到牽制作用。

表7 各區間年均產流產沙量

表9 不同時期水沙變化情況和區間貢獻率

T1到T2時期黃河源區年均輸沙量減少492.71萬t，區間四對輸沙量減少的貢獻率最大，超過其他3個區間對輸沙量貢獻率的總和；T2到T3時期黃河源區年均輸沙量減少2.28萬t，區間四對輸沙量減少的貢獻率最大，達2 872.81%，而其余3個區間則均對輸沙量減少起牽制作用；T1到T3時期黃河源區年均輸沙量減少494.99萬t，區間四對輸沙量減少的貢獻率為81.81%，其余3個區間對輸沙量減少的貢獻率總和不到20%。

2.6 區間輸沙模數

通常用輸沙模數來表示流域內水土流失程度，從表10可以看出，從上游至下游水土流失的嚴重程度逐漸加重，區間四水土流失最為嚴重。但可以看到區間四在整個時期內輸沙模數在減小，此結果與表7年均產沙量的變化一致。

表10 各區間不同時期輸沙模數 t/（km2·a）

2.7 水沙趨勢持續性

用R/S分析法估算Hurst指數以確定未來一段時間水沙變化趨勢，估算的H值見表11。

表11 黃河源區各區間H值

由表11可知，各區間水沙過程的H值均小于0.5，表示黃河源區水沙變化存在反持續性，且泥沙的反持續性稍強于徑流的，對比上文所得MK趨勢檢驗結果，未來一段時間黃河源區各區間水沙將會呈現出上升的趨勢，這一結果與藍永超等［17］所得未來徑流是先揚后抑的結果相一致。再將表11與Hurst指數分級表對照，得到黃河源區各區間徑流的反持續性的強度均為較弱，而泥沙則是區間一和區間四的反持續性較弱，其余區間的反持續性較強。

3 討 論

黃河源區水沙變化的影響因素主要可以歸結為氣候變化和人類活動兩種，氣候變化主要通過降水和氣溫影響水沙變化，人類活動帶來的影響主要是引水用水和對土地覆被的改造等。

徑流對降水變化最為敏感，從降水空間分布來看，黃河源區降水最豐沛的地區在東南部［18］，而區間三主要處于黃河源區的中部和東南部，降水空間上可以解釋區間三是黃河源區的主產流區間，且區間三同時有來自黑河和白河兩大支流的徑流補充，黑河多年平均產流量占區間三產流的16.1%，白河多年平均產流量占區間三產流量的34.3%。

黃河源區T1到T2時期徑流減少的原因已被多位學者探討且得出結論，徑流減少是降水減少、氣溫升高和人類活動加劇共同作用的結果［19］，其中降水是徑流變化的主導因素［20］，尤其在T2時期降水強度和空間變化對徑流減少起到很大作用［21-22］。

觀察到區間二產流量的MK趨勢檢驗結果與其余區間不同，呈微弱上升趨勢，且T1到T2時期年均產流量的變化率是最小的。對此也有學者給出相應的研究結果，黃河源區的冰川凍土主要集中在吉邁以上區域，阿尼瑪卿山冰川1966—2000年面積萎縮21.7 km2，冰川儲量虧損2.66 km3，冰川水資源損失23.9億m3［7，23］。20世紀60年代到90年代約有3.1萬km2的多年凍土轉化為季節凍土［24］，而多年凍土退化對徑流的影響在多年凍土覆蓋率大于40%的流域較為顯著［8］。因此，區間二產流量在T1到T2時期變化率較小是因為受到冰川凍土融水的補給。

區間四在T1到T2時期和T2到T3時期產沙量均在減少，雖然T1到T2時期氣溫升高導致的植被覆蓋度降低和人類活動對地表的破壞會導致產沙量增加［25-26］，但對輸沙量變化影響最大、最顯著的因素是徑流，因此導致T1到T2時期輸沙量減少的主要原因是徑流量減少。而在T2到T3時期區間四產流量回升，而輸沙量依然減少，是因為2005年國家三江源生態修復計劃的實施，2005—2012年三江源區草地覆蓋率明顯提高，尤其在興海縣北部和瑪多縣草地覆蓋率提升10%以上［27］，且有研究表明隨著地表植被覆蓋度的提高，流域產沙量會減少［28］。

4 結 論

（1）1956—2017年黃河源區水沙變化均呈下降趨勢，但并未達到顯著性水平，徑流下降速率為1.21 m3/（s·a），泥沙下降速率為1.62 kg/（s·a）。

（2）使用MK、累積距平、滑動t檢驗結合的方法，檢驗出水沙過程的突變點均為1989年。使用小波分析法對水沙過程進行周期性檢驗，水沙周期變化有很強的相似性，均存在5～7 a、12 a和25 a的周期變化，區別在于泥沙過程12 a的周期震蕩強度遠大于徑流過程的。

（3）對水沙過程進行區間分析，60多a來黃河源區各區間產流量均經歷了T1到T2時期的減少和T2到T3時期的回升；區間一、區間二和區間三輸沙量同樣經歷了T1到T2時期的減少和T2到T3時期的回升，而區間四在3個時期內產沙量在持續減少。

（4）T1到T2時期對徑流量減少貢獻率最大的是區間三，對輸沙量減少貢獻率最大的是區間四；T2到T3時期對徑流量增加貢獻率最大的是區間四，對輸沙量減少貢獻率最大的是區間四。

（5）結合R/S法和MK趨勢檢驗，徑流量、輸沙量在未來一段時間內將會呈上升趨勢。

【責任編輯 張 帥】

［20］ ZHENG H，YUAN J，CHEN L.Short-Term Load Forecasting Using EMD-LSTM Neural Networks with a Xgboost Algorithm for Feature Importance Evaluation［J］.Energies，2017，10（8）:1168-1172.