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

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

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

(1)

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

1 模型結(jié)構(gòu)研究

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

1.1 變異變量定義與公式

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

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

(2)

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

(3)

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

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

(4)

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

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

(5)

向前滑動平均變量SF

(6)

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

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

(7)

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

1.2 變異模型推導(dǎo)

考慮降雨、水流、植被覆蓋度和水利工程指標(biāo)與泥沙的一般關(guān)系：

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

(8)

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

有微分表達(dá)

(9)

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

S=b0Qb1Cb2Wb3Pb4

(10)

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

(11)

類似地有

(12)

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

(13)

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

(14)

和

(15)

由此可得泥沙變異與影響因素變異的關(guān)系：

(16)

式中：bQ,bC,bW,bP為常系數(shù)。類似地可以推導(dǎo)得水流的變異公式：

(17)

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

(18)

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

2 模型檢驗與結(jié)果與分析

2.1 流域特征

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

表1 研究流域基本特征

2.2 精度指標(biāo)

本研究采用的模型計算精度評價指標(biāo)有有效性系數(shù)DC：

(19)

相關(guān)系數(shù)RR：

RR=

(20)

和絕對相對誤差比Re：

(21)

2.3 變異模型結(jié)構(gòu)合理性與模擬有效性分析

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

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

表2 黃河中上游流域水沙變異關(guān)系模型參數(shù)值

表3 溫家川站泥沙變異計算結(jié)果

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

(22)

式中：YN為年份。據(jù)這兩個公式可計算得植被覆蓋率C的年變化,進(jìn)而計算出相應(yīng)的變異值等。

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

(23)

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

圖2 溫家川站實際泥沙年SB與計算泥沙年SBc關(guān)系線

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

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

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

(24)

(25)

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

(26)

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

表5 黃河中上游流域平均變異貢獻(xiàn)比例

3 結(jié) 論

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