黃河中上游流域水沙變異關(guān)系模型研究

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

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

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

(1)

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

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

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

1.1 變異變量定義與公式

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

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

(2)

式中：T為時(shí)間系列變化周期；i表示時(shí)間序列；L為序列長(zhǎng)度。類似地有向前滑動(dòng)平均變量F

(3)

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

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

(4)

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

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

(5)

向前滑動(dòng)平均變量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)

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

S=b0Qb1Cb2Wb3Pb4

(10)

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

(11)

類似地有

(12)

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

(13)

將上文定義的變異變量看作是上式中的dS，dQ，dC，dW，dP微分變量，則SV對(duì)變量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 模型檢驗(yàn)與結(jié)果與分析

2.1 流域特征

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

表1 研究流域基本特征

2.2 精度指標(biāo)

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

(19)

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

RR=

(20)

和絕對(duì)相對(duì)誤差比Re：

(21)

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

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

圖1 溫家川站實(shí)際泥沙年變異ΔSV與計(jì)算泥沙年變異ΔSVc時(shí)間過(guò)程比較

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

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

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

(22)

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

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

(23)

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

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

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

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

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

(24)

(25)

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

(26)

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

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

3 結(jié) 論

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