黃河中上游流域水沙變異關系模型研究

落全富, 包為民, 陳文波, 汪勤海, 孫青雪, 金 樊, 王正華

(1.杭州市水庫管理服務中心, 浙江 杭州 311305; 2.河海大學 水文水資源學院, 江蘇 南京 210098)

全球氣候變化和人類活動影響加劇，導致水文系列均值的變化，使得系列具有不一致性和變異性[1]。特別黃河流域泥沙銳減現象全球矚目，但其減少原因、成因機理至今罕有深入研究[2]，以致泥沙減少是由氣候變化還是水保措施作用所致或各因素的貢獻比例多大等問題一直難以解決[3-8]。黃土高原，由于水土保持措施、退耕還林、生態保護等技術與政策措施作用，河流水沙產生了顯著變化，水沙變化原因已成為研究熱點。目前研究的方法主要有概念水文模型為基礎的分類定量評估法、時間系列變化趨勢分析法和驅動因素分析法[9-10]。用水文概念性模型對泥沙變化原因進行分類定量評估是研究較早的方法，概念性模型具有的物理成因關系，定量確定影響因素變化與泥沙變化量的關系[11-14]。這類方法的優點是物理成因機理清楚，缺點是要求資料多，實際常難以滿足模型要求。如黃河流域泥沙概念模型，涉及超滲產流計算，要求有1 h甚至5 min間隔的時段雨量資料[14-15]，這對于變化原因分析要求有50 a以上的長系列計算，是非常困難的。趙羽西和謝平等的變異分級方法、Mann-Kendall突變檢驗法、Pettitt法、有序聚類分析法和Bernaola-Galvan分割法是趨勢分析法有代表性的成果[1,16-21]。這些方法雖然要求資料較少，但只能分析是否存在變異和確定系列的突變點及變異程度。驅動因素分析法，直接建立影響因素與泥沙變異之間的關系，是近期開始為水流泥沙變異規律研究所重視的。如Wang等[2]根據輸沙量(S)等于降雨量(P)、徑流系數(a)和含沙量(s)相乘的特殊關系S=Pas，推導了輸沙量變化與3者間變化的微分關系：

(1)

該研究開辟了變異規律研究的新途徑，但這關系還存在局限性。首先不能用于一般問題的分析中，特別是人類活動因素對泥沙變化的影響不能直接反映；其次含沙量s與輸沙量S類似都受降雨、徑流系數、人類活動因素的影響；還有徑流系數和含沙量已知相當于輸沙量已知，所以這兩因素同時作為輸沙量關系的自變量相當于自己與自己建立關系，是不合適的，通常含沙量不能作為輸沙量模型的自變量。為此，本研究從水沙變異影響因素分析入手，研究黃河中上游流域水沙變異與影響因素變異間的物理成因關系，定義了變異變量，推導了以微分為基礎的水沙變異關系模型，并用實際水沙系列進行了模擬論證和變化原因定量分析。

1 模型結構研究

水文系列變化通常受氣候和人類活動兩個層面的因素變化影響，氣候因素變化具有周期性，其均值和方差等統計量隨周期變化不大，而水文系列變異是指由均值和方差等特征統計量改變的變化。所以，時間系列變異規律和模型研究，首先要考慮如何消除氣候周期性變化的影響，使得變異規律突顯并簡化變異模型結構。

1.1 變異變量定義與公式

要研究泥沙變異規律和構造變異模型，首先要定義變異變量。對于泥沙時間系列的變異或突變點檢測等研究問題，變異或突變通常是指均值的變化。對于自然因素的均值變化，通常都是連續變化的，只不過變化速度突然加快，給人感覺產生了突變。因此，可以用一定時期的滑動平均構建變異變量，以消除周期內的氣候變化因素影響。

設有時間樣本系列x1,x2,x3,…,xL，可以構建向后滑動平均變量B

(2)

式中：T為時間系列變化周期；i表示時間序列；L為序列長度。類似地有向前滑動平均變量F

(3)

反映均值變化的變量V和變異量ΔV計算公式分別為

Vi=(Bi+Fi)/2， ΔVi=Bi-Fi

(4)

由公式(4)所表達的要素變異反映了滑動平均值的變化量，其值的絕對值越大表示變異越劇烈。

流域泥沙變化，影響因素有降雨、水流、植被、水利工程等。研究泥沙變異規律與模型，就是要研究這些類型因素變異與泥沙變異之間的因果關系。根據上述變異變量定義，對下面物理因素量有相應的表達。如泥沙S向后滑動平均變量SB為

(5)

向前滑動平均變量SF

(6)

反映均值變化的變量SV和變異量ΔSV如下

SVi=(SBi+SFi)/2， ΔSVi=SBi-SFi

(7)

其他植被因素C，降雨P，水流Q，水利工程因素W都可以類似地構建。C均值變化為CVi；C變異量為ΔCVi；P均值變化為PVi；P變異量為ΔPVi；Q均值變化為QVi；Q變異量為ΔQVi；W均值變化為WVi；W變異量為ΔWVi。

1.2 變異模型推導

考慮降雨、水流、植被覆蓋度和水利工程指標與泥沙的一般關系：

S=f(Q,C,W,P)

(8)

式中：P為降雨量(mm)；Q為流量(m3/s)；C為植被覆蓋度,表示植被占流域面積比例；W為水利工程指標,表示水利工程占流域面積比例。

有微分表達

(9)

該微分模型對于任意的泥沙關系都是通用的，關鍵是不知公式(8)的具體結構，所以無法獲得各項偏導數的具體表達函數。根據大量泥沙研究表明[22-26]，通用土壤流失方程是很有效的模型。該模型把土壤侵蝕量表達為：

S=b0Qb1Cb2Wb3Pb4

(10)

式中：b0,b1,b2,b3和b4為常參數。其中S對Q的偏導數據公式(10)可推導得：

(11)

類似地有

(12)

將S對Q，C，W和P求偏導的結果代入公式(9)，得到泥沙與其影響因素的微分模型如下

(13)

將上文定義的變異變量看作是上式中的dS，dQ，dC，dW，dP微分變量，則SV對變量QV，CV，WV，PV的偏導數就可以類似地表達為：

(14)

和

(15)

由此可得泥沙變異與影響因素變異的關系：

(16)

式中：bQ,bC,bW,bP為常系數。類似地可以推導得水流的變異公式：

(17)

式中：cE，cP，cC和cW為常系數； ΔEV為年蒸發變異量，類似地

(18)

公式(15)—(18)就是以微分理論為基礎的水沙變異模型。

2 模型檢驗與結果與分析

2.1 流域特征

本次檢驗選擇了黃河中上游的7個流域，流域特征及測站名稱詳見表1。這些流域面積變化范圍為8 515～682 166 km2，控制黃河流域面積的91%，資料系列較長，普遍在53 a以上，這些流域的模型檢驗，無論是流域尺度、空間地理位置和資料系列長度考慮，都具有很好的代表性。

表1 研究流域基本特征

2.2 精度指標

本研究采用的模型計算精度評價指標有有效性系數DC：

(19)

相關系數RR：

RR=

(20)

和絕對相對誤差比Re：

(21)

2.3 變異模型結構合理性與模擬有效性分析

2.3.1 模型計算結果 變異模型計算，首先要確定氣候周期T。氣候變化主要受太陽活動周期影響，太陽活動周期與太陽黑子爆發密切相關，據研究太陽黑子具有11.2 a的周期，所以這里選擇滑動平均周期T為11。7個斷面流域用最小二乘法確定的參數詳見表2，變異模型計算的結果和精度統計詳見表3—4及圖1—2。表2中SB和SBc是泥沙的觀測與模型計算向后滑動平均量，SF和SFc是泥沙的觀測與模型計算向前滑動平均量，e是計算量的相對誤差，下標s和q分別表示輸沙率和流量，如DCs和DCq分別表示輸沙率和流量模擬有效系數等。表2中植被覆蓋率C為面積比例。用于分析植被變化的遙感數據為GIMMS NDVI和SPOT VGT NDVI兩種數據集，其中GIMMS NDVI數據來自于寒區旱區科學數據中心(http:∥westdc.westgis.ac.cn/)提供的NOAA/AVHRR NDVI半月合成數據集，空間分辨率為8 km×8 km，時間跨度為1984年1月至2006年12月。信息的詳細提取過程參考文獻[27]。根據這些數據確定了植被覆蓋率和水面比例的年變化函數

圖1 溫家川站實際泥沙年變異ΔSV與計算泥沙年變異ΔSVc時間過程比較

表2 黃河中上游流域水沙變異關系模型參數值

表3 溫家川站泥沙變異計算結果

C=0.394e0.009 7(YN-1983)

(22)

式中：YN為年份。據這兩個公式可計算得植被覆蓋率C的年變化,進而計算出相應的變異值等。

這次模型計算中，由于水利工程的年變化信息和流域面平均蒸發資料難以獲得，并為了簡化模型檢驗，先不考慮蒸發和水利工程影響。簡化的模型為

(23)

2.3.2 模擬結果與分析 表3中觀測的年輸沙率向前滑動平均SF與向后滑動平均SB相減得實際輸沙年變異ΔSV，類似地計算的年輸沙率向前滑動平均SFc與向后滑動平均SBc相減得計算輸沙年變異ΔSVc。圖1為窟野河溫家川站實際泥沙年變異ΔSV與計算泥沙年變異ΔSVc的時間過程比較。從時間系列過程看，泥沙最大變異發生在1998年，其值為-732.4 kg/(s·a)，1979年為其次，泥沙變異值為-638.8 kg/(s·a)。從實際泥沙年變異ΔSV與計算泥沙年變異ΔSVc時間過程變化趨勢看，兩者十分吻合。圖2是實際泥沙年SB與計算泥沙年SBc關系，其兩者的關系線接近45°直線，有效性系數DCs為0.974，相關系數RRs為0.987，相對誤差只有6.9%。這些結果都說明微分變異模型的結構是合理的，可有效地表達泥沙的變異關系。

圖2 溫家川站實際泥沙年SB與計算泥沙年SBc關系線

表4是7個流域水流和泥沙變異模擬的精度統計。泥沙系列變異模擬的有效性系數最大的流域為0.981，最小的也有0.709，說明泥沙變異模擬的有效性在這7個流域是很高的；泥沙變異模擬與實際的系列相關系數大于0.900的有4個，占到了57%，最小的也有0.797，說明泥沙年變異ΔSV與計算泥沙年變異ΔSVc間的相關性在7個流域都是十分顯著的；相對誤差最小的為4.8%，最大的也只有23.4%，模擬誤差不同流域不僅都較小且變化不大，說明變異模型模擬關系很穩定。水流系列變異模擬的精度類似，這里限于篇幅不展開討論(表4)。

表4 黃河中上游流域泥沙變異模型計算精度

2.3.3 流域平均變異貢獻比例分析 為了進一步檢驗微分變異模型結構的合理性，這里把變異模型用于泥沙變化貢獻比例分析。由公式(23)—(24)組成的變異模型，有如下變異貢獻比例分析式：

(24)

(25)

式中：前兩項之和反映降雨因素改變導致的泥沙變化,后兩項之和是由植被覆蓋率變化引起的泥沙變異。則氣候因素和植被因素變異引起的貢獻比例分別為

(26)

式中：CHP，CHC分別表示氣候、植被因素變化引起的泥沙變化貢獻比例。具體計算結果詳見表5。從表5結果看,各個流域不同因素貢獻比例略有變化,氣候因素變化貢獻比例在0.106～0.155,植被變化貢獻比例在0.844～0.893,氣候、植被平均貢獻比例分別為0.129和0.870,這與潼關站的貢獻比例十分接近,因為其他6個流域散布在潼關流域內,這全部流域平均與潼關站接近,一方面說明這些流域和相應資料系列的代表性選擇合理；第二方面也進一步證明了變異模型結構的合理性和模擬結果的穩定性。

表5 黃河中上游流域平均變異貢獻比例

3 結 論

本文以微分為基礎，提出直接建立泥沙、水流變異模型，通過黃土高原和黃河中上游7個測站控制流域輸沙率和流量變異變量的模擬檢驗，獲得了7個流域輸沙率和水流變異模擬平均有效性系數為0.860和0.887，平均相對誤差分別為14.8%和6.7%。從水流變異、泥沙變異模擬精度和各因素變異貢獻比例3個層面檢驗了模型，初步證明了結構的合理性和模擬實際泥沙變異具有的高有效性。