落全富, 包為民, 陳文波, 汪勤海, 孫青雪, 金 樊, 王正華
(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)行了模擬論證和變化原因定量分析。
水文系列變化通常受氣候和人類活動兩個層面的因素變化影響,氣候因素變化具有周期性,其均值和方差等統(tǒng)計量隨周期變化不大,而水文系列變異是指由均值和方差等特征統(tǒng)計量改變的變化。所以,時間系列變異規(guī)律和模型研究,首先要考慮如何消除氣候周期性變化的影響,使得變異規(guī)律突顯并簡化變異模型結(jié)構(gòu)。
要研究泥沙變異規(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。
考慮降雨、水流、植被覆蓋度和水利工程指標(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ǔ)的水沙變異模型。
本次檢驗選擇了黃河中上游的7個流域,流域特征及測站名稱詳見表1。這些流域面積變化范圍為8 515~682 166 km2,控制黃河流域面積的91%,資料系列較長,普遍在53 a以上,這些流域的模型檢驗,無論是流域尺度、空間地理位置和資料系列長度考慮,都具有很好的代表性。
表1 研究流域基本特征
本研究采用的模型計算精度評價指標(biāo)有有效性系數(shù)DC:
(19)
相關(guān)系數(shù)RR:
RR=
(20)
和絕對相對誤差比Re:
(21)
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)比例
本文以微分為基礎(chǔ),提出直接建立泥沙、水流變異模型,通過黃土高原和黃河中上游7個測站控制流域輸沙率和流量變異變量的模擬檢驗,獲得了7個流域輸沙率和水流變異模擬平均有效性系數(shù)為0.860和0.887,平均相對誤差分別為14.8%和6.7%。從水流變異、泥沙變異模擬精度和各因素變異貢獻(xiàn)比例3個層面檢驗了模型,初步證明了結(jié)構(gòu)的合理性和模擬實際泥沙變異具有的高有效性。