氣候因子對地表水資源量變化影響的定量分析

陳立華，王　焰，關昊鵬

(1.廣西大學土木建筑工程學院，南寧 530004；2.廣西防災減災與工程安全重點實驗室，南寧 530004)

0　引　言

2017年《北部灣城市群發展規劃》中欽州市作為“雙軸”城市，區域內主要為源短流急中小入海河流，汛期降雨量占全年80%以上，全市各類水庫共393座，總興利庫容僅為4.85 億m3[1]。2012年《中國近海海洋調查評價》指出近90%的城市存在不同程度的缺水問題。在區域人口增長、經濟發展與水資源短缺間矛盾日益凸顯境況下，合理計算濱海區地表水資源量及其未來趨勢變化，剖析在降雨、蒸發、氣溫及其他因子影響下，各因子對徑流變化的貢獻率大小，是分析水資源潛力和正確評價與配置的前提。

深入分析地表水資源量在各因子影響下趨勢變化首要問題是徑流還原。因濱海區不同于內陸，具有河網密布且水文站點不足等特點，導致計算難度偏大。因此在滿足其精度條件下，有必要探究其他徑流還原方法在濱海地區的適用性。第二次全國水資源綜合規劃的水資源調查評價中規定選用的徑流還原方法為：降雨徑流變化趨勢法、雙累積曲線(Double Mass Curve，簡稱DMC)模型、蒸發差值法、分項調查法等，其中DMC模型[2]和蒸發差值法[3]相比于其他方法，具有原理簡單、操作簡便、資料要求低、忽略要素少等特點，因此采用以上兩種方法對徑流進行還原計算，并結合欽州市已有簡化分項調查法對其結果合理性、方法適用性方面進行評價。關于環境變化下各因子對徑流影響的貢獻率分析，目前在流域尺度上，可定量分析氣候變化與其他因子對水文影響的研究方法主要有水文模型法[4-6]和定量評估法。其中定量評估方法主要有氣候彈性系數[7-8]、敏感性分析法[9]、降雨-徑流雙累積曲線法等。水文模型法具有較好的物理基礎，但參數敏感性存在一定的不確定性，若對模擬結果不進行驗證，極可能使氣候變化對徑流的貢獻值偏高[10]；盡管定量評估法所需數據較少，但所需較長數據序列的同時，長序列中噪音會對評估結果造成干擾[8]。王隨繼[11]等于2012年提出累積量斜率變化率比較法(Slope Changing Ratio of Cumulative Quantity，簡稱SCRCQ)法，可有效剔除噪音，較簡便地分離出氣候變化和其他因子對徑流變化的影響比重，在洞庭湖三口河系[12]、長江荊南三口[13]、南水北調中線工程[14]、湟水川流域[15]、黑河流域中上游[16]等地區均得到較好的應用，2017年Wu[17]等人更是將其應用于巖溶區域的氣候變化與人類活動研究。但上述針對SCRCQ法的研究，所假定的基準期均為實測徑流序列，未對其進行還原計算以消除早期人類活動影響，未充分考慮氣溫對降雨與蒸發的影響，未深入探討降雨、蒸發、氣溫對徑流變化影響復雜的內在聯系，這將導致氣候因子貢獻比重偏小，結果誤差偏大。針對研究中存在的問題，有必要對SCRCQ法進行適當改進：基于SCRCQ法的基礎上重新定義基準期，考慮氣溫對降雨與蒸發的影響，并結合水量平衡原理探討各氣候因子對徑流變化影響的內在聯系，并間接獲取其他因子貢獻率值，使SCRCQ法分析結果更貼近實際情況。

1　基礎數據

欽州市地處桂南濱海，屬水資源四級區，境內河網密布，多為中小入海河流。將研究區域劃分為19個子流域如圖1所示，欽州市境內河流四級分區如表1所示。研究所選用資料包括32個雨量站、5個水文站點、3個蒸發站(陸屋、黃屋屯、坡朗坪站)以及1個氣象站(欽州站)分布如圖1所示，研究所采用的降雨、徑流、蒸發、氣象資料統一采用1956-2013年月平均序列。欽江為欽州市最大流域，如圖2所示，是欽南、欽北區及靈山縣最重要的取水來源河流。欽江流域干流全長191 km，集水面積2 326 km2，干流坡降0.36%，地勢東北部高、西南部低(欽江下游濱海平原)。

圖1　欽州市子流域及站網圖Fig.1 Station network and sub basin in Qinzhou

圖2　欽江流域地理位置Fig.2 geographical location Qinjiangbasin

2　分析方法

2.1　氣候變化對水文影響貢獻率分析方法

在SCRCQ法基礎上重新定義基準期，為實測徑流序列首個突變點前的還原徑流序列，研究期為突變點后實測徑流量序列。基準、研究期累積降水量～年份線性方程的斜率分別為Kpa(mm/a)、Kpb(mm/a)，斜率變化率為Sp；累積蒸散發量～年份線性方程的斜率分別為Kra(mm/a)、KEb(mm/a)，斜率變化率為SE；累積徑流量～年份線性方程的斜率分別為KRa(億m3/a)、KRb(億m3/a)，斜率變化率為SR；累積氣溫～年份線性方程的斜率分別為KTa(℃/a)、KTb(℃/a)，斜率變化率為ST。

SP=(KPb-KPa)/KPa

(1)

SE=(KEb-KEa)/KEa

(2)

SR=(KRb-KRa)/KRa

(3)

ST=(KTb-KTa)/KTa

(4)

根據式(1)～式(4)，可求降雨、蒸發、氣溫對徑流變化的貢獻率，分別為CP、CE、CT：

CP=100×SP/SR

(5)

CE=100×SE/SR

(6)

CT=100×ST/SR

(7)

依據所得CP、CE、CT，我們可間接獲取其他因子對徑流變化的貢獻率CO。根據水量平衡原理如式(9)，可知T為自變量，R為因變量，T能夠直接引起R變化，貢獻率為CT；不僅如此，T的改變可引起P、E變化，間接影響到R，貢獻率同樣為CT，但該部分CT對CP、CE進行了重復計算。因此，CO并非定值，而應處于是否考慮T影響下的值域即上下限內變化，如式(8)所示。

CO∈[100-CP-CE-CT, 100-CP-CE]

(8)

2.2　突變點檢測方法

DMC模型與SCRCQ法均需確定基準期，其關鍵在于分析水文序列的突變點位置。水文序列突變檢測方法眾多，其中Mann-Kendall法[1]與累積距平法[20]原理簡單，M-K法毋需要素樣本遵從特定分布，不受異常值的干擾，檢測范圍寬以及理論基礎較強等特點獲得廣泛應用；累積距平法的核心是判斷離散數據對其均值的離散幅度，若累積距平增值大，表明離散數據大于其平均值，反之則小于其平均值[11]。但以上兩種方法檢測突變點存在檢測結果不一致等問題，因此輔以滑動T檢驗法[21]檢測時間序列突變點的顯著性，若檢測超過|t|≥tα=0.05=2.003顯著性水平的時間節點即為突變點，否則為轉折點。

2.3　地表水資源總量計算方法

DMC模型原理如圖3所示，若徑流序列在1986年發生突變，則1986年前序列為基準期①(實心方框)；1986年后為還原期②(空心圓框)。若處于天然條件下則k2=k1，而圖中k1>k2，表明1986年后徑流序列受外界環境改變影響，致使k1減小，因此以①期降雨-徑流累積量線性關系式替代②期進行還原計算。

圖3　DMC模型原理圖Fig.3 Schematic diagram of DMC model

流域蒸發差值法遵循水量平衡原理[18]，地表水資源量R變化與降雨P、流域下墊面變化前的流域蒸發E前、土壤蓄水量W緊密相關，W又可表示為氣溫T與P的函數如式(9)。

R=P-E前-W=f(P,E前,T)

(9)

其中E前可表示為：

E前=[(A-A1)E陸+αA1E水]/A

(10)

式中：A為設計流域總面積，km2；A1為人類活動前流域內的水庫、湖、塘等水面面積，km2；α為水面蒸發折算系數；E水為蒸發皿實測水面蒸發，mm；E陸為陸面蒸發。

研究區域屬南方濱海區，因此陸面蒸發參考凱江公式[19]：

E陸=(αT+β)θ/L

(11)

θ=θ0[0.65-0.35(n-s)]

(12)

式中：θ為太陽總輻射，cal/cm2；T為日平均氣溫，℃；L為蒸發潛熱；α和β為系數和常數；θ0為碧空無云時的太陽總輻射，cal/cm2；n為月平均云量；s為月平均日照量。

3　應用實例分析

3.1　基準期分析

對欽州市內有水文控制站流域徑流量序列采用累積距平法及M-K法進行突變點檢測，主要以欽江流域為例進行詳細介紹。

由圖4(a)可知，UF與UB曲線相交于1986、1993、1994、1996、1997、2001、2003年，累積距平曲線在1964、1986、1992、2003年為極值點，在以上時間節點處徑流量極可能發生突變，但M-K法所檢測出突變點均在α=0.05(±1.96)顯著性水平，累積距平法顯著性無法判斷，因此需輔以滑動T對其進一步顯著性水平檢驗。由滑動T結果可知，僅有1986年與2003年t值分別為2.47和2.55，均超過95%顯著性水平，表明1986、2003年為陸屋站徑流量序列突變年份；對于其他節點并未超過95%顯著性水平，僅為轉折點。采用滑動T對圖4(b)、圖4(c)以及圖4(d)中降雨量、蒸發量、氣溫序列可能的突變節點進行顯著性分析，降雨量序列在整個時間域內并未出現顯著突變點；蒸發量、氣溫序列均僅有一個突變點超過95%顯著性水平，分別為1973年和1996年。

經分析，黃屋屯、坡朗坪、西塘、大江口站徑流量均于1986年發生顯著突變，而面降雨量序列均不存在明顯突變點，因此將1956-1986年假定為DMC模型徑流還原基準期，對1987-2013年序列進行徑流還原并計算欽州市地表水資源總量。

3.2　地表水資源總量及其趨勢分析

3.2.1欽州市地表水資源總量

DMC模型以1956-1986年為基準期，經Pearson相關分析及T統計量相關顯著性分析可知，基準期內5個水文站降雨～徑流相關系數r均超過99%顯著性水平如表1所示，可對1987-2013年徑流序列進行還原計算。蒸發差值法則對1956-2013年整個時段在降雨量的基礎上逐年剔除流域年蒸發量，從而得到水文站天然徑流量。在水文站天然徑流量基礎上采用代表站法分析不同保證率下各流域地表水資源量，結果如表1所示。DMC模型與蒸發差值法所得欽州市地表水資源總量分別為106.573、111.533 億m3。

表1　各典型流域水文控制站降雨～徑流Pearson相關分析Tab.1 Analysis of Pearson correlation between rainfall and runoff in typical watershed hydrological control station

圖4　時間序列M-K法及累積距平法突變點檢測Fig.4 Time series M-K method and cumulative anomaly detection

根據《廣西欽州市水資源綜合規劃報告(2008年)》所采用《廣西壯族自治區地表水資源(1984年)》中1956-1980年實測～天然徑流相關系數，對1980-2000年實測徑流進行還原，所得地表水資源量為104.41 億m3。蒸發差值法嚴格遵循水量平衡原理，對整個序列均進行還原計算，從而消除了早期人類活動對徑流所產生的影響，使地表水資源量更為接近于天然情況。綜上，表明兩種徑流還原方法在欽州市內不僅能夠很好的應用，而且計算簡單、所需資料較少、結果與天然情況較為接近，在南方濱海地區水資源總量計算中值得推廣。

從表2中可知，欽州市地表水資源量最為豐富的為獨流入海河流水資源分區，其地表水資源量占全市總量62.7%～63.3%。由于獨流入海河流源短流急，實際可利用水量不多，如欽江流域多年平均入海水量為17.40 億m3，其實際可利用水量僅為3.082～4.601 億m3。

3.2.2徑流序列的趨勢變化分析

采用3 a滑動平均與Morlet小波綜合對欽州市地表水資源量進行分析如圖5所示，58年序列中存在1960-1973年、1974-1987年、1988-2002年3個平均周期為14 a的完整豐-枯變化過程，1986年開始，降雨、徑流波動幅度較之前明顯增大。盡管整個水資源序列呈0.087 億m3/a的上升趨勢，但自20世紀70年代以來序列卻呈-0.400 億m3/a的下降趨勢。進一步對不同年代均值進行比較可知，各年代均值呈豐-枯轉換階梯式下降趨勢，降雨量序列趨勢性與其相同。綜上所述，可推測地表水資源量在2017-2023年左右將處于偏枯期，2024-2030年左右將處于偏豐期；根據天然與實測徑流量差值曲線如圖6可知，自20世紀60年代以來，各年代天然與實測徑流量差值序列趨向率高達0.406 億m3/a，總體呈明顯梯級上升趨勢，表明以人類活動為主的其他因子對實測徑流減少的影響持續增大，因此有必要選取欽州市典型流域，進一步定量分析降雨、蒸發、氣溫及其他因子對徑流減少影響的貢獻率研究。

圖5　欽州市水資源總量趨勢分析Fig.5 Trend analysis of surfacewater resources in Qinzhou

表2　欽州市各水資源分區地表水資源量Tab.2 Surface water resources quantity of each water resources division in Qinzhou

圖6　天然與實測徑流量差值曲線趨勢分析Fig.6 Analysis of difference between natural and measured runoff

3.3　氣候變化對主要河流徑流減少的貢獻的定量分析

采用改進SCRCQ法對欽州市主要河流定量分析降雨、蒸發、氣溫變化對徑流減少影響的貢獻率的同時，間接分析其他因子對徑流減少影響的貢獻率，其中以欽江流域為例進行詳細介紹，基準期天然徑流量序列則采用蒸發差值法還原結果。

欽江流域各因子突變年份如圖4所示，徑流序列(1986、2003)，降雨量序列無顯著突變點(為便于計算，假定其突變點與徑流序列一致)，蒸發量序列(1973、2003年)，氣溫序列(1996年)。若僅考慮徑流序列突變點情況，將1956-2013年序列統一劃分為3個時期(1956-1986年為AR期，即基準期；1987-2003年為BR期，2004-2013年為CR期，BR與CR期為研究期)，則將使結果偏差較大。因此根據突變點情況，蒸發量、氣溫序列的時期劃分需額外考慮：蒸發量序列AR期僅為1956-1973年，氣溫序列BR期僅為1997-2003年，如圖7所示。

BR與AR期相比，累積徑流深線性關系式斜率減少229.39 mm/a，減少率為23.04%；同期相比累積降雨量線性關系式斜率增加48.11 mm/a，增加率為2.73%；累積蒸發量線性關系式斜率減少74.73 mm/a，減少率為7.56%；累積氣溫線性關系式斜率增加0.97 ℃/a，增加率為4.41%。根據式(9)水量平衡可知地表水資源量變化是降雨量、蒸發量、氣溫的函數，均為氣候因子，而實測徑流量則存在其他因子影響，由其他因子影響較弱的AR期的累積徑流量與年份之關系(圖7)可知，假設BR期不存在其他因素影響，則上述氣候因子斜率變化率之和與徑流斜率變化率相等，而實際上徑流深變化是氣候及其他因子綜合作用的結果。根據式(1)～(7)計算結果可知，BR與AR期相比，降水、蒸發以及氣溫對徑流深減少的貢獻率為11.84%、32.80%和19.14%，由式(8)可知其他因子貢獻率并非定值，而應處于是否考慮氣溫影響下的值域，即[36.21%,55.35%]的上下限內變化。若依據原方法，在基準期仍為實測徑流序列、不考慮氣溫因素的影響，則其他因子對徑流減少的貢獻率為28.24%，氣候因子貢獻率明顯偏大。同理可知，CR與AR期相比，降水量、蒸發量的減少以及氣溫的升高對徑流減少的貢獻率分別為19.02%、5.32%和10.65%，其他因子對徑流減少的貢獻率為65.01%～75.66%的值域內，欽州市主要河流氣候因子對徑流減少的貢獻率如表4所示。

圖7　累積徑流深、降雨量、蒸發量、氣溫與年份之關系Fig.7 Relationship between cumulative runoff depth, rainfall, evaporation, temperature and year

表3　欽江流域氣候因子對徑流變化的貢獻率分析Tab.3 Analysis on contribution rate of climatic factors to runoff variation in Qinjiang Basin

表4　欽州市主要河流氣候因子對徑流變化的貢獻率　%Tab.4 The contribution of climatic factors to runoff change in Qinzhou

由表4可知，BR與AR期相比，位于欽州西南、西、西北部區域的欽江、茅嶺江，徑流減少的貢獻因子主要以其他因子與蒸發為主，降雨、氣溫的影響貢獻率相對較小，而欽州東、北、南以及中部區域的大風江、小江、武思江，徑流減少的貢獻因子則以蒸發為主導，降雨與其他因子貢獻率次之，不同流域受其他因子的貢獻率排序：欽江>茅嶺江>小江>武思江>大風江；CR與AR期相比，各流域徑流減少的主導貢獻為其他因子，具體為：其他因子>降雨>蒸發>氣溫，不同流域受其他因子的貢獻率為：欽江>茅嶺江>小江>大風江>武思江。

綜上所述，各流域徑流減少的主導貢獻因子由蒸發向其他因子轉移，以人類活動為主的其他因子從BR期的欽江流域逐步擴展至CR期的欽州市全境，而降雨、氣溫因子對徑流減少影響的貢獻率在整個研究階段(1956-2013年)內無明顯變化，且居于次要地位。

4　結　語

(1)為確定欽州市主要流域徑流序列基準期，采用M-K、累積距平法并輔以滑動T檢測各時間序列突變點。欽江、茅嶺江、大風江徑流量序列存在1986、2003年兩個突變點顯著，小江、武思江流域為1986年；各流域面雨量序列無明顯突變點，陸屋站蒸發量序列為1973年突變顯著，面蒸發量序列在此基礎上增加一個突變點于2003年。相對于上述突變點而言其他可能的突變點突變特征不顯著，僅在特定年份發生轉折或處于豐枯變化節點上，可視為轉折點。

(2)采用DMC模型及蒸發差值兩種簡化方法進行徑流還原，徑流還原結果偏差控制在9.2%以內，欽州市地表水資源量處于106.573～111.533 億m3范圍內，相對偏差僅為4.4%；自20世紀70年代以來地表水資源總量總體呈-0.400 億m3/a的趨勢下降，而天然與實測徑流量差值曲線呈明顯階梯上升趨勢；地表水資源量存在14 a左右的豐-枯轉換周期，且在2017和2024年左右將分別步入偏枯、偏豐階段。

(3)采用SCRCQ法定量分析了欽州市主要流域地表水資源量受各因子影響下徑流減少的貢獻率，其中欽江流域BR與AR期相比，降雨、蒸發和氣溫對徑流減少的貢獻率分別為11.84%、32.80%、19.14%，其他因子的貢獻率為36.21%～55.35%范圍內；CR與AR時期相比，降雨、蒸發和氣溫對徑流減少的貢獻率為19.02%、5.32%、10.65%，其他因子的貢獻率有所增加為65.01%～75.66%。欽州市主要流域各因子對徑流減少的貢獻研究表明，1986-2003年導致徑流減少的主要貢獻為氣候變化，其中又以蒸發因子最為顯著，而在2004-2013年間，其他因子對徑流減少的貢獻明顯高于氣候變化。

□

參考文獻：

[1]陳立華,王焰,易凱,等.欽州市降雨及入海河流徑流演變規律與趨勢分析[J].水文,2016,36(6):89-96.

[2]申紅彬,吳保生,鄭珊,等.黃河上游蘭州至頭道拐河段流量頻率變化分析[J].水力發電學報,2014,33(1):7-14.

[3]魏茹生.徑流還原計算技術方法及其應用研究[D].西安:西安理工大學,2008.

[4]王國慶,張建云,賀瑞敏.環境變化對黃河中游汾河徑流情勢的影響研究[J].水科學進展, 2006,17(6):853-858.

[5]夏智宏,劉敏,王苗.1990s以來氣候變化和人類活動對洪湖流域徑流影響的定量辨識[J].湖泊科學,2014,26(1):515-521.

[6]賀瑞敏,張建云,鮑振鑫.海河流域河川徑流對氣候變化的響應機理[J].水科學進展,2015,26(1):1-9.

[7]Fu Guobin, Charles Stephen P, Chiew Francis H S.A two-parameter climate elasticity of streamflow index to assessclimate change effects on annual streamflow[J].Water Resources Research, 2007,43(11):2 578-2 584.

[8]Sankarasubramanian A, Vogel R M, Limbrunner J F.Climate elasticity of streamflow in the United States[J].Water Resources Research, 2001,37(6):1 771-1 781.

[9]Jones R N, Chiew F H S, Boughton W C, et al.Estimating the sensitivity of mean annual runoff to climate change using selected hydrological models[J].Advances in Water Resources,2006,29(10):1 419.

[10]Legesse D，Vallet C C, Gasse F.Hydrological response of acatchment to climate and land use changes in Tropical Africa: Case study South Central Ethiopia[J].Journal of Hydrology,2003, 275(1-2):67-85.

[11]王隨繼,閆云霞,顏明,等.皇甫川流域降水和人類活動對徑流量變化的貢獻率分析-累積量斜率變化率比較方法的提出及應用[J].地理學報,2012,67(3):388-397.

[12]李楊,楊文堯.洞庭湖三口河系地區徑流演變情勢與農業用水面臨的挑戰[J].中國農村水利水電,2016,(11):86-92.

[13]帥紅,李景保,何霞,等.環境變化下長江荊南三口徑流變化特征檢測與歸因分析[J].水土保持學報,2016,30(1):83-88.

[14]李凌程,張利平,夏軍,等.氣候變化和人類活動對南水北調中線工程典型流域徑流影響的定量評估[J].氣候變化研究進展,2014,10(2):118-126.

[15]張調風,朱西德,王永劍,等.氣候變化和人類活動對湟水河流域徑流量影響的定量評估[J].資源科學,2014,36(11):2 256-2 262.

[16]何旭強,張勃,孫力煒,等.氣候變化和人類活動對黑河上中游徑流量變化的貢獻率[J].生態學雜志,2012,31(11):2 884-2 890.

[17]WuL H, Wang S J, Bai X Y, et al.Quantitative assessment of the impacts of climate change and humanactivities on runoff change in a typical karstwatershed, SWChina[J].Science of the Total Environment,2017, 601:1 449-1 465.

[18]李林,汪青春,張國勝,等.黃河上游氣候變化對地表水的影響[J].地理學報,2004,59(5):716- 722.

[19]金棟梁,楊世才.用普通氣象資料計算土壤蒸發量的方法(及凱江蒸發公式)[J].人民長江,1981,(4):47-52.

[20]Ran L S, Wang S J, Fan X L.Channel change at Toudaoguai station and itsresponses to the operation of upstream reservoirs in the upper Yellow River[J].Journal of Geographical Sciences, 2010,20(2):231-247.

[21]余予,孟曉艷,張欣.1980-2011年北京城區能見度變化趨勢及突變分析[J].環境科學研究,2013,26(2):129-136.