黃河源區水文要素變化及驅動因素分析

張 鳴，朱 奎，魯 帆，戴雁宇，宋昕熠，于嵩彬

（1.中國礦業大學 資源與地球科學學院，江蘇 徐州 221008；2.中國水利水電科學研究院 水資源研究所，北京 100038；3.長沙理工大學 水利與環境工程學院，湖南 長沙 410114）

0 引言

黃河是我國僅次于長江的第二大河，多年平均水資源總量為647 億m3，不到長江的7%，但水資源開發利用率高達80%，遠超40%的生態警戒線，水資源短缺是黃河流域最大的矛盾［1－3］。黃河源區唐乃亥水文站的多年平均徑流量占黃河水資源總量的34.5%［4］，黃河源區的水資源狀況不僅關系當地的經濟社會發展和生態環境保護，而且對整個黃河流域水資源開發利用有著重要影響。近年來，已有不少學者對黃河源區徑流變化進行了研究，蘇中海等［5］、劉希勝等［6］認為黃河源區上游徑流量呈現不顯著增大趨勢，而劉晶等［7］、Yan 等［8］認為黃河源區徑流量整體呈減小趨勢。各學者研究選取的時間段和數據不同，結論也有所不同。對于徑流驅動因素，商瀅等［9］采用雙累積曲線得出非降水因素是徑流變化的主要驅動因素的結論，張士鋒等［10］通過建立驅動模型得出降水、潛在蒸散發對徑流變化分別起正向、負向驅動作用的結論。上述學者研究徑流時多將研究期看作一個整體，較少考慮研究期內徑流變化過程，驅動因素分析大都只考慮部分氣象要素和人類活動的影響，較少綜合考慮氣象要素和下墊面的影響并定量分析。黃河源區位于素有“亞洲水塔”和“第三極”之稱的青藏高原，對氣候變化較敏感，在全球變暖的環境下以往的研究結論已不完全適用，需補充研究論證。

目前采用Budyko 法、水文模型、雙累積曲線法均能定量揭示氣候變化和人類活動對徑流的影響。采用Budyko 法便于從氣候變化中分離出不同影響因子。近年來，國內外學者采用Budyko 法建立多種經驗公式并在長江和黃河等流域廣泛應用［11－15］，但傳統Budyko法下墊面參數在不同時段會出現跳躍式的變化，這與實際情況存在一定偏差。本文基于1959—2019 年黃河源區水文氣象數據，分析黃河源區水文要素演變規律，運用累積距平曲線和Mann－Kendall 突變檢驗法確定徑流突變年份，采用變參數的Budyko 公式定量分離降水、潛在蒸散發和下墊面變化對徑流量變化的貢獻率，以期為黃河流域合理開發利用水資源、建設水利工程以及防汛抗旱提供參考。

1 研究區概況和數據來源

1.1 研究區概況

黃河流域位于東經96°—119°、北緯32°—42°，發源于青藏高原巴顏喀拉山北麓的約古宗列盆地，流域地勢西高東低，西部河源地區平均高程在4 000 m 以上，源區涉及青海、四川、甘肅3 個省的6 個州18 個縣。黃河源區唐乃亥水文站控制集水面積為12.2 萬km2，黃河源區是世界上高原水體與濕地集中分布的地區之一。研究區內水系、湖泊、10 個氣象站和唐乃亥水文站分布見圖1，為增加站網密度，選取的部分氣象站位于集水面積之外。

圖1 黃河源區范圍

1.2 數據來源

使用的唐乃亥水文站數據源自《黃河水文年鑒》《黃河水資源公報》《黃河泥沙公報》，10 個氣象站數據源自中國氣象數據網，1982—2019 年歸一化植被指數（精度為5 km）源自國家地球系統科學數據中心，1980—2020 年土地利用類型遙感監測數據（精度為1 km）源自中國科學院資源環境科學與數據中心。

2 研究方法

2.1 趨勢檢驗和突變檢驗

趨勢檢驗采用線性回歸分析法和Mann－Kendall檢驗法（簡稱M－K 檢驗）［16］。線性回歸分析法可以反映變量隨時間的變化趨勢。M－K 檢驗是一種非參數檢驗法，其優點是不需要樣本遵從一定的分布規律，也不受少數異常值的干擾。

突變檢驗采用M－K 檢驗法和累積距平法，累積距平法是根據時間序列累積值與均值的差值繪制累積距平曲線，能比較直觀地判斷樣本的變化趨勢和突變年份。

2.2 潛在蒸散發量與植被覆蓋度計算

采用FAO Penman－Monteith 方法計算流域潛在蒸散發量：

式中：E0為潛在蒸散發量，Tmean為日均氣溫，Rn為凈輻射，G為土壤熱通量，γ為干濕表常數，u2為2 m 高風速，es、ea分別為飽和、實際水氣壓，Δ為飽和水氣壓曲線斜率。

李苗苗等［17］基于像元二分模型建立植被覆蓋度計算公式：

式中：fc為植被覆蓋度，N為歸一化植被指數，Nsoil為裸土或者無植被覆蓋區域的歸一化植被指數，Nveg為完全被植被覆蓋區域的歸一化植被指數。

由于存在噪聲等其他因素的影響，因此取5%的N值為Nsoil，取95%的N值為Nveg。為更好地反映研究區內植被變化情況，歸一化植被指數可選取年內最大值。

2.3 徑流歸因分析

基于Budyko 假設的長時間序列中，水分輸入可看作降水量P，能量輸入可看作潛在蒸散發量E0，考慮水量平衡和能量平衡的水熱耦合方程為

式中：E為流域蒸散發量。

對于閉合流域來說，長期水量平衡方程為

式中：Q為徑流量，ΔS為流域水分蓄變量。

考慮到區域氣候、地理環境等的差異，國內外學者基于Budyko 理論陸續建立了各種考慮下墊面參數的經驗公式［18－21］，本文選用Zhang 等［18］建立的經驗公式：

式中：ε為干旱指數，ε＝E0/P；n為下墊面參數。

采用式（5）分析徑流變化時，下墊面參數和水文要素會在不同時段呈現跳躍性變化，這與實際情況存在偏差。Jiang 等［15］提出了變參數的Budyko 公式，選取氣象和人類活動等要素為協變量構建下墊面參數的相關函數，假定選取m個自變量，滑動平均序列為y1，y2，…，yt，第i（i＝1，2，…，t）年的m個自變量為，，…，，可選用一次函數、二次函數、指數函數等擬合下墊面參數和m個自變量之間的關系，第i年的徑流量Qi為

式中：Pi為第i年降水量，εi為第i年干旱指數，βi為第i年回歸系數，ei為第i年殘差。

假定ei符合正態分布N（0，σ2），則Qi符合正態分布，根據參數β和σ的似然函數求出極大似然估計量和，對應有ni＝Yi，ni為第i年的下墊面參數。

為定量描述氣象因子和下墊面參數對徑流量變化的影響，引入彈性系數φ來表示徑流量對某因子x變化的敏感程度：

假設變化期和基準期內降水量、潛在蒸散發量、下墊面參數變化分別為ΔP、ΔE0、Δn，根據基準期信息推求變化期內徑流變化量：

3 水文要素變化

基于唐乃亥水文站實測徑流數據，分析1959—2019 年黃河源區年徑流量變化情況，見圖2（a）。M－K趨勢檢驗結果顯示，整體上年徑流量呈現減小趨勢，線性減小速率約為2.36 億m3/10 a，多年實測年均徑流量為203.9 億m3，1959 年年徑流量（為328.08 億m3）和2002 年年徑流量（為105.74 億m3）分別為61 a樣本系列中徑流量的最大值和最小值。利用M－K 檢驗方法對徑流序列進行突變檢驗，見圖2（b），在0.05顯著性水平下，UBK曲線與UFK曲線在1987—1990 年置信區間內多次交匯。為明確徑流突變年份，引入年徑流量累積距平曲線，見圖2（c），年徑流量在1959—1989 年整體呈現上升趨勢，在1989—2019 年整體呈現下降趨勢，在1989 年發生一次明顯突變，綜合判斷1989 年為徑流突變年份。為進一步反映近年來全球變暖環境下黃河源區的徑流變化，繪制1990—2019 年年徑流量累積距平曲線，見圖2（d），2004 年為第2 個突變年份。總體上，從1959—1989 年、1990—2004 年和2005—2019 年3 個時段來看，年徑流量呈現先下降后上升趨勢。

圖2 1959—2019 年黃河源區徑流量變化

對1959—2019 年黃河源區年內徑流量進行分析，1981 年9 月月徑流量最大（為92.02 億m3），2003 年1月月徑流量最小（為2.33 億m3），每年7 月月均徑流量最大（為35.04 億m3），每年1 月月均徑流量最小（為4.58 億m3），7—9 月徑流量占全年總徑流量的46.97%，而1—3 月徑流量僅占全年總徑流量的7.32%，由此可見，黃河源區徑流量季節分配不均。徑流集中度指月徑流量以矢量形式累積，其各分矢量之和占年徑流量的百分比，反映年內徑流量的集中程度。集中期反映全年徑流量集中的月份［22］。唐乃亥水文站年內徑流集中度和集中期計算結果見圖3，集中度為35%～55%，整體上呈現下降趨勢，下降率為0.716%/10 a。集中期為7—9 月，整體上呈現略微提前的趨勢，對于唐乃亥水文站和下游水庫來說，洪峰有提前的趨勢。

圖3 唐乃亥水文站年內徑流集中度和集中期計算結果

基于10 個氣象站降水數據分析黃河源區1959—2019 年降水量變化趨勢，見圖4。1959—2019 年多年平均降水量為504.83 mm，1962 年降水量（為398.24 mm）和2018 年降水量（為674.21 mm）分別為61 a 樣本系列中的最小值和最大值。M－K 趨勢檢驗結果顯示，統計量Z為3.72＞2.56，降水量呈現顯著增加趨勢（顯著性水平α＝0.01），平均線性增加速率為9.1 mm/10 a。

圖4 1959—2019 年降水量變化

徑流系數代表流域產流能力，綜合反映流域內自然地理要素對降水—徑流關系的影響。1959—2019年徑流系數變化情況見圖5，1959—1989 年、1990—2004 年和2005—2019 年對應的徑流系數均值分別為0.36、0.28 和0.32，徑流系數呈現先減小后增大的趨勢，與徑流量變化趨勢相同。

圖5 1959—2019 年徑流系數變化情況

4 驅動因素

4.1 降水、氣溫和潛在蒸散發分析

基于黃河源區10 個氣象站數據對降水量進行Kriging 插值分析，得出多年平均降水量的空間分布，見圖6（a），黃河源區1959—2019 年多年平均降水量在空間上呈現由西北向東南遞增的趨勢。M－K 突變檢驗結果［見圖6（b）］和降水量累積距平曲線［見圖6（c）］顯示，2002 年和2016 年為降水序列的突變點。1959—2002 年、2003—2016 年和2017—2019 年多年平均降水量分別為490.67、525.97、613.84 mm，尤其是近年來降水量增加趨勢愈發明顯。

圖6 降水量變化

基于10 個氣象站實測數據分析黃河源區1959—2019 年年均氣溫變化趨勢，見圖7（a）。多年平均氣溫為－0.66 ℃，1965 年年均氣溫最低（為－2.02 ℃），2016 年年均氣溫最高（為0.73 ℃），M－K 趨勢檢驗統計量Z為11.40＞2.56，氣溫呈現顯著上升趨勢（α＝0.01），平均線性上升速率為0.3 ℃/10 a。由氣溫累積距平曲線［見圖7（b）］可以看出，氣溫累積距平值在1997 年發生轉折，結合趨勢檢驗結果，氣溫在1959—1997 年呈現不顯著上升趨勢，在1997 年之后開始顯著上升。

圖7 1959—2019 年氣溫變化趨勢與氣溫累積距平曲線

基于10 個氣象站實測氣溫、風速、日照和相對濕度數據，利用FAO Penman－Monteith 方法計算各氣象站潛在蒸散發量，采用Thiessen 多邊形方法計算出黃河源區年平均潛在蒸散發量，見圖8（a）。1959—2019年多年平均潛在蒸散發量為724.8 mm，1965 年潛在蒸散發量最小（為672.4 mm），2010 年潛在蒸散發量最大（為784.5 mm）。M－K 趨勢檢驗統計量Z為7.89＞2.56，潛在蒸散發量呈現顯著增加趨勢（α＝0.01），平均線性增加率為8.2 mm/10 a。潛在蒸散發量主要受日照時數、風速、氣溫和相對濕度的影響，因此利用多元線性回歸法對每個氣象站的潛在蒸散發量進行歸因分析。結果表明，氣溫是90%氣象站潛在蒸散發量的主要影響因素，風速是10%氣象站潛在蒸散發量的主要影響因素。M－K 突變檢驗結果［見圖8（b）］顯示，潛在蒸散發量在1997 年發生一次突變，突變時間點與氣溫的相同，說明氣溫與潛在蒸散發量之間具有較強的相關性。

圖8 1959—2019 年潛在蒸散發量變化趨勢與M－K 突變檢驗

為分析1959—2019 年黃河源區降水量、氣溫、潛在蒸散發量與徑流量的相關性，對其進行單因素相關性檢驗，結果顯示，徑流量與降水量之間的皮爾遜系數為0.778，具有極強的相關性，遠遠大于氣溫的和潛在蒸散發量的。

4.2 下墊面變化分析

下墊面主要指與大氣直接接觸的地球表面，本文主要分析黃河源區植被、冰雪和凍土的變化情況。植被覆蓋度是衡量地表植被變化的重要指標，黃河源區植被覆蓋度變化趨勢與累積距平曲線見圖9。M－K 趨勢檢驗統計量Z為2.01＜2.56，植被覆蓋度呈現不顯著上升趨勢（α＝0.01），這與管曉祥等［23］、劉啟興等［24］、韓思淇等［25］的研究結果一致。在空間上植被覆蓋度呈現由西北向東南遞增的趨勢，與降水量的空間分布一致，這與高思琦等［26］認為黃河源區植被的空間分布主要受降水影響的觀點一致。

為分析黃河源區內不同植被類型的變化情況，對1980—2020 年土地利用類型數據進行分析，見圖10。源區內土地利用類型主要為草地，城鎮用地和湖泊面積變化不大，而林地、草地分別在1990—1995 年、1996—2000 年呈現略微退化趨勢。2000—2013 年未利用土地面積呈現顯著減小趨勢，而林地面積和草地面積呈現顯著增大趨勢。植被變化是人類活動和氣候因子共同作用的結果，在生態環境較為敏感的青藏高原地區植被變化更為顯著。氣候變化對植被的影響主要通過降水量和氣溫的變化改變土壤濕度和溫度，調節植物生長過程中酶的活性，影響植物光合作用和呼吸作用［27－29］。降水量增加會促使高寒草原和荒漠植被生長，增加植被面積，氣溫升高則促使植被覆蓋度提高［30］，植被的改善會增強區域蒸散發，進而增加降水量。

圖10 1980—2020 年土地利用類型變化

歸一化植被指數是表征區域植被變化的重要指數，1982—2019 年黃河源區降水量與歸一化植被指數的相關性分析及顯著性檢驗結果見圖11。78%的區域降水量對歸一化植被指數起正向促進作用，兩者成顯著正相關的區域主要集中在降水相對貧乏的西北部地區和東北部地區；21%的區域降水量與歸一化植被指數負相關，兩者成顯著負相關的區域很少，植被類型以林地和耕地為主，降水較充沛，這些地區的植被變化可能受人類活動影響較大，這與史丹丹等［31］的研究結果一致。管曉祥等［23］研究發現黃河源區氣溫與歸一化植被指數顯著正相關，且春季的影響最大。黃河源區多年平均降水量為504.83 mm，降水較充沛，因此水量不是限制植被生長的主要條件。黃河源區高程普遍在4 000 m以上，熱量不足限制植被生長，因此氣溫是影響黃河源區植被變化的主要因子［23，27－28，32］。此外，黃河源區的人類活動也對植被變化產生一定影響，如過度砍伐和放牧、退耕還林、減畜工程和成立三江源國家自然保護區等。

圖11 1982—2019 年降水量與歸一化植被指數相關性分析

黃河源區冰雪面積變化情況見圖12，累積距平曲線顯示，1982—1997 年冰雪面積呈現增加趨勢，1997—2015 年冰雪呈現消融趨勢，冰雪面積變化的突變時間為1997 年，與氣溫的一致，表明氣溫是控制冰雪變化的主要因素［33］。

圖12 1982—2015 年冰雪面積占比累積距平曲線

關于黃河源區凍土的變化情況，曹慧宇［34］的研究表明2003—2019 年多年凍土面積由3.59 萬km2減小到3.42 萬km2，4.82%的多年凍土轉化為季節性凍土，季節性凍土最大凍結深度減小16 cm，活動層厚度增加9.28 cm。馮雨晴［35］的研究表明1961—2015 年黃河源區凍土融化夾層厚度和活動層厚度呈現增加趨勢，1997 年之后呈現顯著增大趨勢。多年凍土的退化以及活動層厚度的增加導致土壤滲透性增強，地表徑流更多入滲成為地下水，土壤含水量整體增加，蒸散發量增加，流域地下水儲量增加，從而減少夏季徑流量，增加冬季徑流量和基流量，調節徑流量的季節性分布［36－40］。此外，活動層底部的凍結層上水可能側向補給江河湖泊，黃河源區地下冰融水對地表徑流的補給比例為13.2%～16.7%，但其融水對徑流的改變較小［41－44］。

4.3 Budyko 假設定量分離

根據徑流突變結果，以1989 年和2004 年為突變年份將研究期分為3 個時段，3 個時段的多年平均徑流深、降水量、潛在蒸散發量、氣溫見表1。

表1 黃河源區各時段徑流深和協變量

為消除氣候變化及地下水儲量對水量平衡的影響，構建下墊面參數時滑動平均步長取9 a，利用式（5）計算滑動平均后的下墊面參數n。基于Jiang等［15］提出的變參數Budyko 公式，選取降水量、潛在蒸散發量和氣溫作為協變量，選用線性函數進行擬合，下墊面參數n與降水量、潛在蒸散發量和氣溫的關系式為

式（9）標準化系數顯示降水量是下墊面參數的主要影響因素，其次是氣溫。根據式（5）和式（4）利用氣象數據和徑流數據推求氣候變化造成的徑流變化量dQc，再根據dQh＝dQ－dQc計算分析人類活動對徑流的影響（dQh為人類活動對徑流的改變量，dQ為徑流總改變量）。結果顯示，在第1 個突變期（1959—2004年）氣候變化對徑流量減小的貢獻率為86.63%、人類活動對徑流量減小的貢獻率為13.37%，在第2 個突變期（1990—2019 年）氣候變化對徑流量增加的貢獻率為101.18%、人類活動對徑流量增加的貢獻率為－1.18%。在整個研究期內氣候變化是徑流量變化的主導因素，人類活動對徑流量的增大一直起負向作用，但對徑流量變化的影響程度明顯降低。

為定量分離不同氣象因子和下墊面變化的貢獻率，分別計算降水量、潛在蒸散發量和下墊面參數在不同時期對應的彈性系數，結果見表2。基準期和變化期（分別以1959—1989 年為基準期、1990—2004 年為變化期，1990—2004 年為基準期、2005—2019 年為變化期）降水量的彈性系數始終大于潛在蒸散發量的和下墊面參數的，說明徑流對降水的響應最為敏感。

表2 黃河源區下墊面參數和各協變量對應的彈性系數

根據式（8）計算突變期各因子對徑流量變化的貢獻率，在第1 個突變期內下墊面參數n對徑流量變化的貢獻率達到60.69%，是徑流量減小的主要驅動因素，降水其次（為31.75%），潛在蒸散發的貢獻率最小（為7.56%）；第2 個突變期內降水對徑流量變化的貢獻率達到103.91%，是徑流量增加的主要驅動因素，潛在蒸散發其次（為－10.37%），下墊面變化的貢獻率最小（為6.46%）。此外，計算得出2 個突變期內各因子引起的模擬徑流變化量之和與實測徑流變化量的誤差分別為3.98%和0.1%，誤差較小，表明結果具有準確性。

5 結論

1）黃河源區1959—2019 年徑流整體呈現先下降后上升的變化趨勢，多年實測年均徑流量為203.9 億m3，年內徑流集中度為35%～55%，集中期位于7—9 月，年內徑流量季節分配不均，1989 年和2004 年為徑流量突變年份。1959—2019 年年降水量和潛在蒸散發量都呈現顯著增加趨勢，氣溫在1997 年之后顯著上升。1982—2019 年植被覆蓋度呈現不顯著上升趨勢，在空間上植被覆蓋度呈現由西北向東南遞增的趨勢，與降水量的空間分布一致。

2）基于變參數的Budyko 方法定量分離不同因子對徑流量變化的影響，氣候變化是黃河源區1959—2019 年徑流量變化的主要影響因素，徑流對降水的響應最為敏感。1959—2004 年下墊面參數n對徑流量變化的貢獻率達到60.69%，是徑流減小的主要驅動因素，降水其次，潛在蒸散發的貢獻率最小。1990—2019年降水對徑流量變化的貢獻率達到103.91%，是徑流增大的主要驅動因素，潛在蒸散發次之，下墊面的貢獻率最小。

3）徑流量與降水量之間的皮爾遜相關系數為0.778，遠遠大于氣溫的和潛在蒸散發量的，因此1959—2019 年黃河源區徑流量改變最主要的影響因素是降水，其次為氣溫。