基于分位數(shù)回歸的娘子關(guān)泉降水及徑流變化分析

王永珂，范永輝

（天津師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院，天津300387）

娘子關(guān)泉域地處山西省陽泉市，是當(dāng)?shù)刂饕墓┧矗瑢ζ浣邓皬搅髯兓?guī)律的研究具有重要的理論與實(shí)踐意義.現(xiàn)已有許多文獻(xiàn)對不同水域降水、徑流之間的響應(yīng)關(guān)系進(jìn)行了研究，主要采用趨勢分析[1-2]或水文模型[3-4]等方法.分位數(shù)回歸模型最早由Koenker等提出[5]，之后隨著對該模型的研究和推廣，其廣泛應(yīng)用于經(jīng)濟(jì)、醫(yī)學(xué)、教育、生物學(xué)等領(lǐng)域.相比其他方法，分位數(shù)回歸能夠更充分反映自變量對不同部分因變量的分布產(chǎn)生的不同影響，即它能在很大程度上提供數(shù)據(jù)的全部信息.分位數(shù)回歸的估計(jì)量不易受到異常值的影響，具有很強(qiáng)的穩(wěn)健性，且分位數(shù)回歸不需要對模型中的隨機(jī)干擾項(xiàng)做任何分布的假定，使模型具有很強(qiáng)的適應(yīng)性.近些年，該方法被逐步應(yīng)用于水文氣象領(lǐng)域.文獻(xiàn)[6]將分位數(shù)回歸方法和變異點(diǎn)分析結(jié)合起來，對美國南部月降水的可能變化情況進(jìn)行了分析.文獻(xiàn)[7]運(yùn)用分位數(shù)回歸方法研究了波羅的海海平面的長時(shí)間變化情況，結(jié)果表明，斜率在極大值處表現(xiàn)出更大的趨勢性.文獻(xiàn)[8]運(yùn)用分位數(shù)回歸模型分析了灤河流域降水及徑流變化特性.

對于娘子關(guān)泉域降水、徑流變化特征的研究，文獻(xiàn)[9]采用廣義極值Pareto分布模型分析月最大降水量，并在剔除人類活動(dòng)造成的泉水流量下降趨勢（平穩(wěn)化）后，研究不同重現(xiàn)期的極大月降水量與平穩(wěn)化后的泉水年平均流量的關(guān)系，以獲得極端降水與泉水流量的相關(guān)關(guān)系.為了更詳細(xì)地分析降水、徑流的變化情況，本文采用分位數(shù)回歸方法，在選取的高、低分位數(shù)水平下對年降水、徑流序列的變化特征進(jìn)行研究.

1 研究區(qū)概況及數(shù)據(jù)處理

1.1 研究區(qū)概況

娘子關(guān)泉位于山西省與河北省交界處的平定縣娘子關(guān)附近，是我國北方最大的巖溶泉之一.其地理位置介于北緯 36°55′～37°50′、 東經(jīng) 112°20′～113°55′之間.涉及的行政區(qū)有陽泉市郊區(qū)、平定縣、盂縣，以及晉中市的和順縣、左權(quán)縣、昔陽縣和壽陽縣等7個(gè)區(qū)域.娘子關(guān)泉出露于桃河與溫河匯集地段，由11個(gè)主要泉組組成，分布自程家村至綿河的葦澤關(guān)約7 km長的河漫灘及階地上，出露標(biāo)高372～392 m，泉域面積為4 667 km2.

娘子關(guān)泉域?qū)倥瘻貛Т箨懶约撅L(fēng)氣候，多年平均降水量為534.3 mm，月降水量變化較大，最大月降水量達(dá)到341.83 mm（1963年8月）.多年平均徑流量為10.4 m3/s（1956—2003年），其中最小月流量為4.69 m3/s（1995年3月），最大月流量為17.30 m3/s（1963年10月）.

1.2 數(shù)據(jù)處理

考慮到不同地區(qū)降水對娘子關(guān)泉地下水補(bǔ)給的貢獻(xiàn)率[10]以及不同地區(qū)的面積，首先采用泰森多邊形加權(quán)平均的方法求得2個(gè)時(shí)間段內(nèi)的降水權(quán)重.1959—1971年6個(gè)雨量站的降水權(quán)重分別為：陽泉（0.38）、盂縣（0.35）、壽陽（0.09）、昔陽（0.10）、和順（0.05）、左權(quán)（0.03）.1971—2009年 7個(gè)雨量站的降水權(quán)重分別為：陽泉（0.30）、盂縣（0.15）、壽陽（0.09）、昔陽（0.10）、和順（0.05）、左權(quán)（0.03）、平定（0.28）.

2 研究方法

2.1 分位數(shù)回歸方法

給定一個(gè)連續(xù)型一元隨機(jī)變量Y，其分布函數(shù)為FY（y），則其τ分位數(shù)定義為

其中 τ∈（0，1）.

設(shè)給定X下Y的τ條件分位數(shù)[11]為

其中FY|X（y）為連續(xù)型隨機(jī)變量Y在給定X下的條件分布.顯然，Qτ（Y|X）是X的函數(shù)，記為

取g（X）為線性函數(shù)，即

其中 β（τ）是未知參數(shù).參數(shù) β（τ）滿足如下關(guān)系

批而未供和閑置土地大清查大處置工作雖然取得了階段性成果，但推進(jìn)批而未供土地消化利用和閑置土地處置工作不能松懈。要舉一反三，常抓不懈，以本次行動(dòng)為契機(jī)推進(jìn)長效機(jī)制的建設(shè)。下一步，重點(diǎn)抓好以下幾個(gè)方面工作。

其中：ρτ（z）=z（τ-I（z<0））稱為對號函數(shù)（check function），I（·）為示性函數(shù).容易看出，ρτ（z）是一個(gè)非對稱的損失函數(shù).當(dāng)τ=1/2時(shí)，式（5）即為最小一乘回歸.

給定觀測數(shù)據(jù)（x1，y1），（x2，y2），…，（xn，yn），分位數(shù)回歸通過最小化損失函數(shù)來求解參數(shù)估計(jì)值.由下式求解τ分位數(shù)的回歸估計(jì)：

式（6）也可以寫作：

參數(shù)β（τ）的估計(jì)可以通過下式得到：

2.2 分布變化分析方法

多條分位數(shù)回歸線表現(xiàn)出的回歸趨勢會(huì)引起不同年份概率分布的變化.因此，需要在得到不同分位點(diǎn)下回歸線的基礎(chǔ)上，對年降水和年徑流序列在不同年份的概率分布變化情況[12]進(jìn)行分析.具體的分布變化分析過程為：根據(jù)分位數(shù)回歸結(jié)果，可以得到不同τ值下的分位數(shù)值.選定一個(gè)年份，根據(jù)年降水和徑流序列在不同τ值下的分位數(shù)值求得分布函數(shù).然后，根據(jù)求得的分布函數(shù)，運(yùn)用自適應(yīng)核密度估計(jì)方法得到概率密度函數(shù).

3 結(jié)果與分析

3.1 年降水和年徑流的分位數(shù)回歸分析

為了能夠更明顯地顯示高低分位數(shù)水平下年降水指標(biāo)隨時(shí)間的變化趨勢，在0～1之間每間隔0.05選取一個(gè)分位點(diǎn)，其中，分位點(diǎn)為0.50時(shí)的回歸直線為中位數(shù)回歸.另外，為了便于對比不同的回歸線，擬合了一條線性回歸直線.結(jié)果見圖1～圖3，圖中自上而下的分位點(diǎn)分別為：0.95，0.90，0.85，0.80，…，0.05.

圖1 年降水量隨時(shí)間的變化趨勢Fig.1 Changes trend of annual precipitation over time

圖2 年徑流量隨時(shí)間的變化趨勢Fig.2 Changes trend of annual runoff over time

圖3 年降水量和年徑流量之間的關(guān)系Fig.3 Relationship between annual precipitation and annual runoff

觀察年降水量隨時(shí)間的變化趨勢（圖1）可見，中位數(shù)回歸直線和線性回歸直線幾乎重合.絕大多數(shù)分位點(diǎn)的回歸直線呈下降趨勢，當(dāng)分位點(diǎn)特別低時(shí)，年降水量隨時(shí)間變化不明顯，而當(dāng)分位點(diǎn)為0.05時(shí)，年降水量有顯著的下降趨勢.估算0.05和0.95這2個(gè)分位點(diǎn)處的年降水量可得，1959年的年降水量大約分布在406 mm和826 mm之間，2009年的近似最小和最大年降水量分別為264 mm和629 mm.

觀察年徑流量隨時(shí)間的變化趨勢（圖2）可見，中位數(shù)回歸直線在1995年之前一直位于最小二乘法所得到的回歸直線之下，這說明1995年之前的年流量是右偏的.另外，高分位點(diǎn)的分布隨時(shí)間的變化更為集中，而低分位點(diǎn)的分布隨時(shí)間的變化更為發(fā)散，分位數(shù)回歸直線左側(cè)和右側(cè)之間的間隙相差不大，這意味著數(shù)據(jù)點(diǎn)整體較為均勻.回歸直線均呈下降趨勢，說明年徑流量在逐年減少.估算0.05和0.95這2個(gè)分位點(diǎn)處的年徑流量可得，1959年的年徑流量大約分布在150 m3/s和193 m3/s之間，2009年的近似最小和最大年徑流量分別為37 m3/s和80 m3/s.

由年降水量和年徑流量之間的關(guān)系（圖3）可見，在絕大多數(shù)分位點(diǎn)處，隨著年降水量的增加，年徑流量呈增加的趨勢.而當(dāng)年徑流量特別低時(shí)，年徑流量隨著年降水量的增加反而呈下降趨勢.

3.2 年降水和年徑流的分布變化

為了解年降水和年徑流序列的分布變化規(guī)律，在百分位數(shù)線性回歸的基礎(chǔ)上，對序列在特定年份的經(jīng)驗(yàn)累積分布函數(shù)（CDF）以及經(jīng)驗(yàn)概率密度函數(shù)（PDF）進(jìn)行了分析.由于采用的方法是線性分位數(shù)回歸方法，所以其他年份的CDF和PDF的變化范圍介于第一年和最后一年之間.因此，選取時(shí)間序列中第一個(gè)年份（1959年）和最后一個(gè)年份（2009年）作為代表年份進(jìn)行研究.結(jié)果見圖4和圖5.

圖4 年降水量的經(jīng)驗(yàn)累積分布函數(shù)（CDF）和經(jīng)驗(yàn)概率密度函數(shù)（PDF）Fig.4 Empirical cumulative distribution function（CDF）and empirical probability density function（PDF）of annual precipitation

圖5 年徑流量的經(jīng)驗(yàn)累積分布函數(shù)（CDF）和經(jīng)驗(yàn)概率密度函數(shù)（PDF）Fig.5 Empirical cumulative distribution function（CDF）and empirical probability density function（PDF）of annual runoff

由年降水量的CDF和PDF曲線（圖4）可見，2009年的PDF與1959年的相比，不僅峰度發(fā)生了變化，從1個(gè)峰度變化為2個(gè)峰度，其分布的范圍也發(fā)生了改變.概率密度在1959年的最高值為0.003 5左右，而到2009年變?yōu)?.004 0左右，函數(shù)分布范圍從年降水量400～800 mm左右減少至300～600 mm左右.由圖1，年降水量分位數(shù)回歸線有匯聚并且下降的趨勢，因此PDF曲線向左移動(dòng).

由年徑流量的CDF和PDF曲線（圖5）可見，2009年的PDF與1959年的相比，不僅峰度發(fā)生了變化，從2個(gè)峰度變化為1個(gè)峰度，其分布的范圍也發(fā)生了改變.概率密度在1959年的最高值為0.03左右，而到2009年變?yōu)?.08左右，峰值顯著變大，函數(shù)分布范圍從年徑流量150～200 m3/s左右減少至20～80 m3/s左右.由圖2，年徑流量分位數(shù)回歸線有聚合的變化趨勢，且回歸直線斜率均為負(fù)值，下降趨勢明顯，因此PDF曲線有很大程度的左移.

4 結(jié)論

對娘子關(guān)泉域1959—2009年的年降水、年徑流資料運(yùn)用分位數(shù)回歸方法進(jìn)行分析，研究年降水和年徑流之間的變化特征及其響應(yīng)關(guān)系.在此基礎(chǔ)上，對序列在特定年份的分布函數(shù)及概率密度函數(shù)進(jìn)行分析，從而了解娘子關(guān)泉年降水、年徑流序列的分布變化，得出如下結(jié)論：

（1）年降水序列在絕大多數(shù)分位點(diǎn)的回歸直線有下降趨勢.高分位點(diǎn)的下降趨勢總體來說比低分位點(diǎn)的下降趨勢更為顯著.年徑流序列高分位點(diǎn)的回歸直線隨時(shí)間有收斂趨勢，低分位點(diǎn)的回歸直線隨時(shí)間有發(fā)散趨勢，但是整體來說，回歸直線都呈下降趨勢，說明年徑流量在逐年減少.

（2）經(jīng)估算，1959年的年降水量大約分布在406～826 mm之間，而2009年的年降水量大約分布在264～629 mm之間，1959年的年徑流量大約分布在150～193 m3/s之間，而2009年的年徑流量大約分布在37～80 m3/s之間.

（3）在絕大多數(shù)分位點(diǎn)處，年徑流量隨著年降水量的增加有增加的趨勢.而當(dāng)年徑流量特別低時(shí)，年徑流量隨著年降水量的增加反而呈下降的趨勢，這可能是因?yàn)槭艿搅巳祟惢顒?dòng)的影響.

（4）2009年與1959年相比，年降水序列的CDF分布范圍縮減，PDF曲線向左移動(dòng)，且從1個(gè)峰度變?yōu)?個(gè)峰度.函數(shù)分布范圍從年降水量400～800 mm左右減少至300～600 mm左右.年徑流序列的PDF曲線左移，從2個(gè)峰度變?yōu)?個(gè)峰度，且峰值顯著變大，函數(shù)分布范圍從年徑流量150～200 m3/s左右減少至20～80 m3/s左右.