董 潔，譚秀翠，張慶華，刁艷芳，徐 妍

(山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東 泰安 271018)

0 引 言

水利建設(shè)規(guī)劃設(shè)計(jì)時(shí)，要確定一個(gè)合理的設(shè)計(jì)標(biāo)準(zhǔn)或工程規(guī)模。一般采用水文頻率分析方法計(jì)算設(shè)計(jì)值。目前主要有參數(shù)和非參數(shù)數(shù)理統(tǒng)計(jì)方法[1]。

用參數(shù)法推求設(shè)計(jì)值xp通常必須解決好兩個(gè)基本問(wèn)題：首先，必須確定水文變量的概率分布模型，即線型選擇；其次，估計(jì)所選線型中的未知參數(shù)即參數(shù)估計(jì)[1]。由于水文現(xiàn)象影響因素的復(fù)雜性，有些水文學(xué)者從統(tǒng)計(jì)理論上分析論證水文變量服從的分布，提出了很多線型供選擇，例如P-Ⅲ型分布、對(duì)數(shù)正態(tài)分布和K-M分布等線型，但是水文變量的真實(shí)分布是很難確切表達(dá)的。在實(shí)際水文工作中，常根據(jù)實(shí)測(cè)經(jīng)驗(yàn)點(diǎn)據(jù)和頻率曲線擬合的好壞選擇線型，由于實(shí)際應(yīng)用中評(píng)判擬合優(yōu)劣的標(biāo)準(zhǔn)各異，所得結(jié)論往往相差較大[2]。根據(jù)我國(guó)水文特點(diǎn)，一般水文變量采用P-Ⅲ型分布。每一種概率分布中都包含若干參數(shù)，選定線型后還必須確定其中的參數(shù)才能進(jìn)行頻率計(jì)算。但是這些參數(shù)是無(wú)法根據(jù)水文現(xiàn)象的物理機(jī)制確定的，必須利用實(shí)測(cè)資料加以估計(jì)。因此，國(guó)內(nèi)外水文學(xué)者研究了很多符合水文特點(diǎn)的參數(shù)估計(jì)方法，主要有適線法，權(quán)函數(shù)法和概率權(quán)重矩法等。由于水文系列長(zhǎng)度短，且所需推求的是稀遇的設(shè)計(jì)值等原因，參數(shù)統(tǒng)計(jì)估計(jì)方法并不理想。

非參數(shù)統(tǒng)計(jì)中的概率密度估計(jì)是指在給定實(shí)測(cè)樣本后，對(duì)其總體概率密度的估計(jì)。非參數(shù)概率密度估計(jì)始于直方圖法，后來(lái)發(fā)展為最近鄰法、核估計(jì)法等，其中理論發(fā)展最完善的是概率密度核估計(jì)法[3]。

本文把非參數(shù)回歸理論應(yīng)用于水文頻率計(jì)算中，建立了水文頻率分析的非參數(shù)回歸變換模型，通過(guò)變換把樣本變換成“最佳樣本”，“最佳樣本”對(duì)應(yīng)的“最佳函數(shù)”具有較高的估計(jì)精度。

1 非參數(shù)回歸變換理論

設(shè)(X1,Y1),…,(Xn,Yn)是R1×1維隨機(jī)變量，滿足模型：

Yi=m(Xi)+εi,i=1,2,…,n

(1)

式中：εi是均值為零的獨(dú)立隨機(jī)變量。

回歸的目的就是根據(jù)樣本(Xi,Yi)估計(jì)Y關(guān)于X的條件概率密度的均值：

m(x)=E(Y|X=x)

(2)

Watston給出了條件均值的核估計(jì)：

(3)

由于核回歸是兩個(gè)量的比，直接求核回歸的偏差比較困難，可以分別考慮上式分子分母的偏差。

(4)

其中：[m(·)f(·)]″(x)=m″(x)f(x)+2m′(x)f′(x)+m(x)f″(x)。

將原樣本Xi變換為均勻分布的新樣本Ui，目的是對(duì)Yi進(jìn)行加權(quán)平均時(shí)，權(quán)的大小只與樣本的鄰近次序有關(guān)，與距離無(wú)關(guān)，這體現(xiàn)了各樣本互相獨(dú)立的原則。只要知道Xi的分布函數(shù)F(x)，就能進(jìn)行變換Ui=F(Xi)。

這里假定Xi已經(jīng)變換成均勻分布的樣本Ui，現(xiàn)在只討論因變量樣本Yi的變換。設(shè)Vi是均值為1的隨機(jī)樣本，可以進(jìn)行如下“加法”迭代。

加法迭代的基本原理是將回歸曲線分解成兩條曲線的和，第一條曲線接近回歸曲線，第二條曲線接近常數(shù)0。迭代方法如下：

(5)

式中：l是迭代次數(shù)。

2 水文頻率分析非參數(shù)回歸變換模型的建立

2.1 分布函數(shù)的非參數(shù)變換估計(jì)

分布函數(shù)估計(jì)與密度估計(jì)一樣有許多重要的作用。當(dāng)樣本數(shù)量足夠大時(shí)，經(jīng)驗(yàn)分布函數(shù)可以近似逼近總體的分布函數(shù)，但是如果樣本數(shù)量較少，表現(xiàn)為階梯函數(shù)，使得經(jīng)驗(yàn)分布函數(shù)估計(jì)受到限制。如果想由密度函數(shù)估計(jì)得到分布函數(shù)估計(jì)，密度函數(shù)估計(jì)應(yīng)是無(wú)偏的，方差大點(diǎn)不要緊，因?yàn)榉e分時(shí)會(huì)消除方差的影響。若要從分布函數(shù)得到密度函數(shù)估計(jì)，分布函數(shù)應(yīng)連續(xù)光滑，有點(diǎn)偏差不要緊，數(shù)值微分是用鄰近幾個(gè)點(diǎn)估計(jì)一個(gè)點(diǎn)的微分值，對(duì)非光滑點(diǎn)(方差大的點(diǎn))很敏感。因此，一個(gè)合理的想法是利用經(jīng)驗(yàn)分布的無(wú)偏性加以平滑就得到偏差和方差都較小的估計(jì)。這種方法是非參數(shù)回歸的一個(gè)特例，不同之處在于因變量不是事先得到的樣本，而是分布函數(shù)的無(wú)偏估點(diǎn)計(jì)。因此分布函數(shù)估計(jì)的第一步，是構(gòu)造無(wú)偏估計(jì)點(diǎn)作為回歸的因變量。下面采用經(jīng)驗(yàn)分布函數(shù)作為分布函數(shù)的無(wú)偏估計(jì)。

2.2 樣本經(jīng)驗(yàn)分布函數(shù)

設(shè)x1,x2,…,xn,是隨機(jī)樣本，有分布函數(shù)F(x)，考慮水文的特點(diǎn)，取經(jīng)驗(yàn)分布函數(shù)是：

(6)

式中：n為樣本的個(gè)數(shù)；m為樣本x大于等于某個(gè)給定值的個(gè)數(shù)。

對(duì)上式的經(jīng)驗(yàn)分布函數(shù)定義：除x為樣本點(diǎn)外的Fn(x)仍和式(6)一樣，而x=xi處的分布函數(shù)估計(jì)值為：

(7)

x1≥x2≥…≥xn

其中，α是[0，1]區(qū)間的常數(shù)，通常取α=1。

由式(6)和式(7)定義的函數(shù)稱(chēng)為樣本經(jīng)驗(yàn)分布函數(shù)。

提出樣本經(jīng)驗(yàn)分布函數(shù)的理由如下：先考慮只有一個(gè)樣本即n=1，樣本出現(xiàn)對(duì)應(yīng)一次貝努利試驗(yàn)，出現(xiàn)的概率是α/2，也就是樣本位于α/2分位點(diǎn)。若樣本取自中心較高且具有對(duì)稱(chēng)密度函數(shù)的母體，α/2=0.5意味著單個(gè)樣本位于分布的中點(diǎn)，對(duì)應(yīng)最大的概率密度。此時(shí)的經(jīng)驗(yàn)分布函數(shù)符合最大似然準(zhǔn)則。對(duì)于n個(gè)樣本，由于各樣本是獨(dú)立同分布，樣本分布函數(shù)值將分布函數(shù)等分。樣本應(yīng)位于等分后的每小塊中的α/2分位點(diǎn)上，即各樣本等間隔地位于(i-1-α)/(n+1)分位點(diǎn)上。這個(gè)想法與離散分布中的古典概率論一致。

顯然，樣本經(jīng)驗(yàn)分布函數(shù)的性質(zhì)與經(jīng)驗(yàn)分布函數(shù)一致，由式(6)知，樣本經(jīng)驗(yàn)分布函數(shù)是小于等于樣本值的樣本數(shù)與樣本總數(shù)加一之比，因此nF(xi)是二項(xiàng)分布隨機(jī)變量，則有樣本經(jīng)驗(yàn)分布函數(shù)的均值和方差為：

E[Fn(xi)]=F(xi)

varFn(xi)=F(xi)[1-F(xi)]/n

(8)

所以樣本經(jīng)驗(yàn)分布函數(shù)是無(wú)偏估計(jì)，這為回歸提供了好的先決條件。在樣本經(jīng)驗(yàn)分布的基礎(chǔ)上回歸，可使分布函數(shù)的方差減小。另外，以樣本經(jīng)驗(yàn)分布函數(shù)為固定點(diǎn)進(jìn)行內(nèi)插，也可得到比較光滑連續(xù)的曲線，但是這樣得到的曲線保留了樣本經(jīng)驗(yàn)分布函數(shù)的方差，這個(gè)方差的最明顯的證明是對(duì)內(nèi)插得到的曲線進(jìn)行微分，得到的密度函數(shù)估計(jì)存在許多毛刺和突起，微分對(duì)方差的敏感是檢驗(yàn)方差存在的直觀方法。

2.3 樣本經(jīng)驗(yàn)分布函數(shù)的變換回歸

由前面定義知，樣本經(jīng)驗(yàn)分布函數(shù)是：

(9)

i=1,2,…,n

記yi=Fn(xi) ，則(xi,yi)構(gòu)成R1×1維隨機(jī)變量，可以進(jìn)行非參數(shù)回歸計(jì)算。但這里的yi不是通?；貧w中事先給出的樣本，是分布函數(shù)在xi點(diǎn)的無(wú)偏估計(jì)。將原樣本(xi,yi)變成新樣本(ui,vi)，其中vi是均值為1的隨機(jī)變量，ui是均勻分布隨機(jī)變量。

由于樣本Xi不是等間隔的固定點(diǎn)，變換回歸包含自變量X和因變量Y的兩步變換，將原樣本(xi,yi)變成新樣本(ui,vi)。在分布函數(shù)估計(jì)中，這兩步變換可以一起完成，因?yàn)閅i和Xi有一定的關(guān)系。

(10)

此時(shí)的Yi相當(dāng)于U1,i近似于45°斜線，則：

V1,i=Yi-U1,i+1

(11)

變換回歸的偏差很小，基本上是無(wú)偏回歸，因樣本經(jīng)驗(yàn)分布函數(shù)也是無(wú)偏估計(jì)，所以最后得到的分布函數(shù)估計(jì)是無(wú)偏估計(jì)。

3 蒙特卡羅統(tǒng)計(jì)試驗(yàn)分析

蒙特卡羅統(tǒng)計(jì)試驗(yàn)分析[6,7]可以很好的說(shuō)明模型的穩(wěn)健性問(wèn)題。假設(shè)水文變量分別來(lái)自P-Ⅲ型分布和對(duì)數(shù)正態(tài)分布總體，用參數(shù)統(tǒng)計(jì)的適線法和非參數(shù)核回歸變換法比較分析模型對(duì)總體和樣本數(shù)量的穩(wěn)健性。

3.1 試驗(yàn)方案

為了考察模型對(duì)樣本數(shù)量和不同設(shè)計(jì)頻率的穩(wěn)健性，分別選取了樣本容量n1=40；n2=50；n3=60；設(shè)計(jì)頻率p1=0.01,p2=0.001。

3.2 計(jì)算方案

用蒙特卡羅方法生成P-Ⅲ型分布和對(duì)數(shù)正態(tài)分布LN樣本，為了考察模型對(duì)樣本總體和計(jì)算方法的穩(wěn)健性，設(shè)計(jì)了108種方案。分別采用非參數(shù)的核回歸變換法DRM；參數(shù)統(tǒng)計(jì)的適線法CFM估計(jì)設(shè)計(jì)值。計(jì)算結(jié)果見(jiàn)表1和2。(只列出樣本容量為n=60的結(jié)果)

3.3 評(píng)價(jià)指標(biāo)

(1)均方誤差σ和相對(duì)均方誤差δ：

(2)絕對(duì)誤差θ和相對(duì)誤差?。

3.4 統(tǒng)計(jì)試驗(yàn)結(jié)果

統(tǒng)計(jì)試驗(yàn)結(jié)果見(jiàn)表1、表2。

3.5 估計(jì)方法的比較分析

表1和表2分別用4個(gè)評(píng)價(jià)指標(biāo)計(jì)算了P-Ⅲ型分布和對(duì)數(shù)正態(tài)分布(LN)，擬合理論總體是P-Ⅲ型分布和對(duì)數(shù)正態(tài)分布(LN)的穩(wěn)健性，分析結(jié)果如下：

(1)總體是P-Ⅲ型分布。當(dāng)用P-Ⅲ型分布曲線擬合P-Ⅲ型分布總體時(shí)，?xp1≤1%；?xp2≤1%δxp1≤15%δxp2≤20%。

當(dāng)用對(duì)數(shù)正態(tài)分布(LN)擬合P-Ⅲ型分布，總體適線的穩(wěn)定性較差，?xp1≥120；?xp2≥20%δxp1≥25%δxp2≥25%。

當(dāng)用非參數(shù)的密度變換法DEM估計(jì)設(shè)計(jì)值，百年一遇和千年一遇的設(shè)計(jì)值相對(duì)誤差一般介于以上兩種情況之間，?xp1≤15%；?xp2≤15%δxp1≤15%δxp2≤20%。

(2)總體是對(duì)數(shù)正態(tài)分布(LN)。當(dāng)用對(duì)數(shù)正態(tài)分布(LN)擬合對(duì)數(shù)正態(tài)分布總體，?xp1≤8%；?xp2≤8%δxp1≤15%δxp2≤20%。

當(dāng)用P-Ⅲ型分布擬合對(duì)數(shù)正態(tài)分布LN總體適線的穩(wěn)定性較差，一般求出的設(shè)計(jì)值偏小，很不安全。δxp1≥15%δxp2≥20%。

表1 P-Ⅲ型分布理論總體計(jì)算成果表(p1=0.01,EX=1 000,p2=0.001)Tab.1 Calculation result of P-Ⅲ theory distribution (p1=0.01, EX=1 000, p2=0.001)

注：各組真值xp如下：Cv=0.5，Cs=1.5時(shí)，xp1=2 665,xp2=3 617；Cv=0.5,Cs=2.0時(shí)，xp1=2 803,xp2=3 954；Cv=1，Cs=2.5時(shí)，xp1=4 875，xp2=7 548。

表2 對(duì)數(shù)正態(tài)分布理論總體計(jì)算成果表(p1=0.01,EX=1 000,p2=0.001)Tab.2 Calculation result of Logarithmic normal theory distribution (p1=0.01, EX=1 000, p2=0.001)

注：各組真值xp如下：Cv=0.5，Cs=1.5時(shí)，xp1=2 665，xp2=3 755；Cv=0.5，Cs=2.0時(shí)，xp1=2 760，xp2=4 120；Cv=1，Cs=2.5時(shí)，xp1=4 672，xp2=7 877。

非參數(shù)的回歸變換法DEM估計(jì)設(shè)計(jì)值，百年一遇和千年一遇的設(shè)計(jì)值相對(duì)誤差一般介于以上兩種情況之間，?xp1≤15%；?xp2≤15%δxp1≤15%δxp2≤20%。

當(dāng)樣本容量由40、50增加到60，相對(duì)誤差和相對(duì)均方誤差減小，說(shuō)明樣本容量增加，穩(wěn)定性增加。

4 應(yīng)用實(shí)例研究

為了更好地檢驗(yàn)該模型的適用性,以下對(duì)具有代表性的小浪底站77 a；宜昌站132 a；崗南站57 a；平山站30 a和黃壁莊站56 a的年最大洪峰流量系列，用P-Ⅲ適線法和變換回歸模型進(jìn)行了頻率分析計(jì)算,各站年最大洪峰流量系列資料表詳見(jiàn)文獻(xiàn)[8]，計(jì)算成果見(jiàn)表3。通過(guò)以上計(jì)算可以看出，用非參數(shù)回歸變換模型計(jì)算出的設(shè)計(jì)值比適線法計(jì)算的值普遍偏大，頻率P大于百年以上的設(shè)計(jì)值較參數(shù)適線法計(jì)算的設(shè)計(jì)值大3%～7%；頻率P小于百年以下的設(shè)計(jì)值較參數(shù)適線法計(jì)算的設(shè)計(jì)值大2%～3%。

表3 小浪底站、平山等站年最大洪峰流量頻率分析結(jié)果Tab.3 Flood frequency analysis of the XIAOLANGDI station and PINGSHAN station

5 結(jié) 語(yǔ)

(1)我國(guó)水文頻率計(jì)算一般采用參數(shù)統(tǒng)計(jì)的“適線法”。需要事先假定水文變量服從某種線型分布，如果假設(shè)線型和實(shí)際分布不符就會(huì)影響計(jì)算結(jié)果。非參數(shù)統(tǒng)計(jì)方法不需要假定線型，但是理論上要求適用于大樣本。本文建立的非參數(shù)回歸變換方法通過(guò)對(duì)小樣本進(jìn)行變換，不僅拓展了非參數(shù)統(tǒng)計(jì)的應(yīng)用，還提高了小樣本的估計(jì)精度。

(2)用蒙特卡羅統(tǒng)計(jì)試驗(yàn)方法考察模型的穩(wěn)健性問(wèn)題。通過(guò)對(duì)水文計(jì)算中常用的P-Ⅲ型和對(duì)數(shù)正態(tài)分布兩個(gè)總體，用參數(shù)統(tǒng)計(jì)的“適線法”和非參數(shù)統(tǒng)計(jì)的“變換回歸”兩種估計(jì)方法比較分析了三組參數(shù)、三組樣本。結(jié)果表明，模型是穩(wěn)健的，方法是合理的。

(3)以我國(guó)幾個(gè)有代表性的測(cè)站的年最大洪峰流量為例，比較研究變換回歸模型和參數(shù)適線法在水文頻率分析計(jì)算中的結(jié)果特征，通過(guò)計(jì)算得到用非參數(shù)方法計(jì)算的水文頻率結(jié)果較適線法偏大，多種計(jì)算結(jié)果可以互相參考，為水利工程規(guī)劃提供了一種有效方法。