李弘瑞 李新杰 張紅濤

摘 要：為了掌握黃河中游水沙關(guān)系的變化特點(diǎn)和變化趨勢，基于黃河頭道拐站1958—1990年33 a和1998—2017年20 a的水沙資料，通過P-Ⅲ型曲線分別建立年徑流量與年輸沙量的頻率曲線，對年徑流量和年輸沙量進(jìn)行豐平枯劃分;基于Copula函數(shù)對水沙關(guān)系進(jìn)行變異診斷，構(gòu)建年徑流量-年輸沙量聯(lián)合分布模型并通過擬合優(yōu)度檢驗(yàn)選取最合適的Copula函數(shù)，計(jì)算水沙豐平枯組合遭遇頻率。結(jié)果表明：頭道拐水文站在1958—1990年和1998—2017年的年徑流量和年輸沙量呈下降趨勢，1990年水沙關(guān)系發(fā)生變異;1958—1990年水沙關(guān)系存在一個(gè)變異點(diǎn)（1965年），1998—2017年水沙關(guān)系不存在明顯的變異點(diǎn)，整體水沙關(guān)系存在一定的平穩(wěn)性;1966—1990年和1998—2017年兩個(gè)時(shí)期的水沙豐平枯組合遭遇頻率中豐枯同頻高于豐枯異頻，1966—1990年的年徑流量和年輸沙量在較大值處相關(guān)性較強(qiáng)，1998—2017年年徑流量和年輸沙量整體相關(guān)性較強(qiáng)，水量和沙量的豐平枯變化趨勢基本一致。

關(guān)鍵詞：水沙關(guān)系;聯(lián)合分布;變異診斷;遭遇頻率;頭道拐水文站;黃河

Abstract： In order to grasp the characteristics and trends of the relationship between runoff and sediment at the upper and middle reaches of the Yellow River， based on the water and sediment volume data of the Yellow River Toudaoguai Station in the 33 years from 1958 to 1990 and the past 20 years from 1998 to 2017， the P-Ⅲ curve was adopted to establish the runoff and sediment discharge frequency and divide annual runoff and annual sediment transport into hight， medium and low flow. It carried out mutation diagnosis of runoff and sediment relationship， and joint distribution model of annual runoff and annual sediment transport was constructed based on copula function by the goodness of fit test. Finally， it calculated the encounter frequency of wet-normal-dry runoff and sediment combination. The results show that the annual runoff and sediment discharge of Toudaoguai Hydrological Station in 1958-1990 and 1998-2017 show a downward trend and the relationship between runoff and sediment is remarkable especially after 1990; There is a variation point （1965） in the relationship between water and sediment during 1958-1990. There is no obvious variability point in the water and sediment relationship between 1998 and 2017 and there is a certain stability in the overall runoff and sediment relationship; the synchronous frequency in the frequency of the combination of wet-normal-dry runoff and sediment is greater than that the asynchronous frequency during 1966-1990 and 1998-2017， the annual runoff and annual sediment transport in the previous period have a strong correlation at the annual sediment load is larger. The correlation between annual runoff and annual sediment transport in the latter period is stronger and the variation trend of annual runoff volume has similar variation tendency to the annual sediment.

Key words： runoff and sediment relationship; joint distribution; mutation diagnosis; encounter frequency; Toudaoguai Hydrological Station; Yellow River

1 引 言

隨著黃河上游水利工程的運(yùn)行和中游寧蒙河段灌區(qū)取用水等人類活動(dòng)的影響，進(jìn)入黃河中游的水沙過程發(fā)生了很大變化，對黃河中下游河道及水庫沖淤、流域生態(tài)環(huán)境等產(chǎn)生了重要影響[1]。因此，科學(xué)揭示水沙關(guān)系變化規(guī)律對黃河中下游水資源管理和高效利用具有重要意義。

針對黃河流域水沙關(guān)系，有關(guān)學(xué)者開展了相關(guān)研究[2-4]。金鑫等基于水沙序列的理論頻率曲線提出了黃河中游不同水文站的水沙頻率組合劃分標(biāo)準(zhǔn)[5];李新杰等從保證率的角度建立了潼關(guān)站1961—2018年58 a的年徑流量和輸沙量的頻率曲線，并確定了潼關(guān)站來水來沙豐平枯的劃分標(biāo)準(zhǔn)[6];馬雁等基于黃河上游蘭州站2014—2018年的水沙數(shù)據(jù)，總結(jié)了年徑流量和年輸沙量的關(guān)系和演變規(guī)律[7];考慮到環(huán)境變化和人類活動(dòng)對流域來水來沙量的影響，李艷玲等結(jié)合滑動(dòng)窗口和Copula函數(shù)，提出基于滑動(dòng)Copula函數(shù)的降水和徑流關(guān)系變異診斷法[8];郭愛軍等提出滑動(dòng)相關(guān)系數(shù)法用于涇河流域水沙關(guān)系的變異診斷，并基于Copula函數(shù)分析了徑河流域的水沙關(guān)系演變特征[9];姚曼飛等基于涇河干支流8個(gè)水文站實(shí)測水沙數(shù)據(jù)，采用Copula函數(shù)構(gòu)建了水沙聯(lián)合分布模型，并計(jì)算了各個(gè)水文站的水沙豐平枯遭遇頻率[10]。

鑒于頭道拐站水沙關(guān)系對黃河中下游水沙調(diào)控的重要作用，以及已有研究多集中于徑流量和輸沙量特點(diǎn)或水沙關(guān)系單方面性質(zhì)變化，對水沙組合遭遇情況分析稍有不足，筆者統(tǒng)計(jì)頭道拐站1958—1990年和1998—2017年共計(jì)53 a的年徑流量和年輸沙量原始數(shù)據(jù)，通過建立水沙頻率曲線，診斷水沙關(guān)系的變異特性，計(jì)算水沙豐平枯組合的遭遇頻率，基于二維Copula函數(shù)建立二維水沙聯(lián)合分布函數(shù)，分析黃河頭道拐站水沙關(guān)系整體的變化情況及豐平枯遭遇規(guī)律，以期為流域水沙調(diào)控提供技術(shù)支撐。

2 研究數(shù)據(jù)與方法

2.1 研究數(shù)據(jù)

頭道拐水文站位于內(nèi)蒙古托克托縣雙河鎮(zhèn)，是黃河上游的出口站，也是中游萬家寨水庫的入口站，其徑流主要來自蘭州站以上，泥沙主要來自蘭州至頭道拐區(qū)間支流，其水沙變化直接影響萬家寨水庫的水沙調(diào)度情況[11-12]。本研究選取頭道拐水文站1958—1990年和1998—2017年兩個(gè)時(shí)期的年徑流量和年輸沙量數(shù)據(jù)，資料來源于黃河流域水文年鑒。

2.2 研究方法

2.2.1 頻率曲線

通過對年徑流量和年輸沙量序列的計(jì)算得到基本參數(shù)、CV、CS，進(jìn)而確定參數(shù)α、β和b的值，即確定P-Ⅲ型分布函數(shù)。

本文運(yùn)用適線法對頻率曲線進(jìn)行尋優(yōu)，通過P-Ⅲ型頻率分布建立年徑流量和年輸沙量的理論頻率曲線，應(yīng)用數(shù)學(xué)期望公式計(jì)算年徑流量和年輸沙量的經(jīng)驗(yàn)頻率，計(jì)算理論頻率和經(jīng)驗(yàn)頻率擬合時(shí)的確定性系數(shù)R2，其計(jì)算公式如下[15]：

2.2.2 Copula函數(shù)

Copula函數(shù)是把隨機(jī)變量X1、X2、…、XN的聯(lián)合分布函數(shù)與各自的邊緣分布函數(shù)U1、U2、…、UN相連接的連接函數(shù)，即函數(shù)C（u1，u2，…，uN）[16]。

2.2.3 Copula函數(shù)類型及參數(shù)估計(jì)

Copula函數(shù)在將邊緣分布進(jìn)行聯(lián)合分布時(shí)包含了變量所有的相依信息，在建立聯(lián)合分布函數(shù)的過程中保持信息的真實(shí)度[17]，更好地反映出水沙序列邊緣分布之間的關(guān)系變化。Copula函數(shù)在總體上大致分為3種：正態(tài)型、t型和Archimedean型。在水文統(tǒng)計(jì)中，主要使用Archimedean型Copula函數(shù)對水沙序列進(jìn)行聯(lián)合分布，具體的表達(dá)式及參數(shù)范圍見表1[17]。

2.2.4 Copula函數(shù)擬合優(yōu)度檢驗(yàn)

為了檢驗(yàn)Copula函數(shù)對于水沙序列聯(lián)合分布的擬合優(yōu)度，引入經(jīng)驗(yàn)Copula函數(shù)，比較經(jīng)驗(yàn)累計(jì)頻率與理論累計(jì)頻率的擬合度。將樣本值代入Copula頻率公式計(jì)算出經(jīng)驗(yàn)累計(jì)頻率，同樣將樣本值代入3類Archimedean型Copula函數(shù)計(jì)算理論累計(jì)頻率。計(jì)算確定性系數(shù)R2并進(jìn)行比較，再根據(jù)離差平方和OLS最小準(zhǔn)則與AIC信息準(zhǔn)則做進(jìn)一步的擬合優(yōu)度檢驗(yàn)[18]。

通過擬合優(yōu)度檢驗(yàn)選取最合適的Copula函數(shù)建立水沙序列的聯(lián)合分布，從而得到接近實(shí)際值的理論分布。

2.2.5 變異診斷

本文基于Copula函數(shù)對水沙序列中水沙關(guān)系在時(shí)間分布上的變異情況進(jìn)行診斷[21]。構(gòu)造基于Copula函數(shù)的極大似然對數(shù)統(tǒng)計(jì)量公式如下：

原序列經(jīng)過一次診斷后，若發(fā)現(xiàn)變異點(diǎn)則依據(jù)變異點(diǎn)劃分為兩個(gè)子序列，對子序列進(jìn)行二次診斷，重復(fù)以上操作直至子序列的長度小于25。若原序列長度小于25，一次診斷即可。

對原水沙序列進(jìn)行變異診斷是為了發(fā)現(xiàn)水沙序列在時(shí)間分布上是否存在平穩(wěn)性，若存在，對原序列的水沙關(guān)系進(jìn)行整體分析;若不存在，則根據(jù)變異點(diǎn)將原序列劃分為若干個(gè)子序列，對每個(gè)子序列的水沙關(guān)系獨(dú)立分析。

2.2.6 水沙豐平枯組合遭遇頻率劃分

在《水文基本術(shù)語和符號標(biāo)準(zhǔn)》[23]中，將河流的豐平枯年頻率劃分為5個(gè)級別，通過不同的保證率確定豐平枯年的頻率劃分標(biāo)準(zhǔn)，見表2[24]。

由于徑流量和輸沙量之間相依度高、相關(guān)性強(qiáng)，一般認(rèn)為輸沙量頻率分布與徑流量頻率分布類似，因此輸沙量的豐平枯年的頻率劃分也可以采用表2的標(biāo)準(zhǔn)[25]。通過建立滿足P-Ⅲ型分布的年徑流量和年輸沙量頻率分布，依據(jù)頻率劃分標(biāo)準(zhǔn)確定豐平枯年年徑流量和年輸沙量的劃分標(biāo)準(zhǔn)，把年徑流量和年輸沙量均劃分為5個(gè)級別，水沙豐平枯組合具體劃分標(biāo)準(zhǔn)見表3。

3 數(shù)據(jù)分析

3.1 確定水沙頻率曲線

本文選取頭道拐站1958—1990年和1998—2017年兩個(gè)時(shí)段的年徑流量和年輸沙量序列，通過P-Ⅲ型分布函數(shù)建立年徑流量和年輸沙量序列的理論頻率分布，同經(jīng)驗(yàn)頻率分布進(jìn)行擬合，在滿足確定性系數(shù)R2≥0.95的前提下，得到擬合度較高的年徑流量和年輸沙量的理論頻率曲線，見圖1和圖2。兩個(gè)時(shí)段年徑流量和年輸沙量頻率變化情況見圖3。

3.2 水沙聯(lián)合分布及擬合優(yōu)度檢驗(yàn)

在確定合適的邊緣分布后，作為變量輸入Clayton、Gumbel和Frank Copula函數(shù)以及經(jīng)驗(yàn)Copula函數(shù)，分別建立年徑流量和年輸沙量的聯(lián)合分布與經(jīng)驗(yàn)聯(lián)合分布并進(jìn)行擬合，根據(jù)確定性系數(shù)R2、OLS和AIC信息準(zhǔn)則3個(gè)擬合優(yōu)度檢驗(yàn)指標(biāo)選取最合適的Copula函數(shù)，3類Archimedean型Copula函數(shù)的參數(shù)及擬合優(yōu)度檢驗(yàn)值的計(jì)算結(jié)果見表4。

根據(jù)表4的檢驗(yàn)值可知：1958—1990年的聯(lián)合分布Clayton Copula的R2值最大、OLS和AIC的值均最小，即Clayton Copula的擬合優(yōu)度最佳，表明Clayton Copula函數(shù)是用于頭道拐站1958—1990年水沙聯(lián)合分布最合適的Copula函數(shù);1998—2017年的聯(lián)合分布中Frank Copula的R2值最大、OLS和AIC的值均最小，說明Frank Copula的擬合優(yōu)度最佳，明顯優(yōu)于其他兩類，F(xiàn)rank Copula函數(shù)是用于頭道拐站1998—2017年水沙聯(lián)合分布最合適的Copula函數(shù)。

3.3 水沙關(guān)系的變異診斷

基于Gumbel Copula函數(shù)對1958—1990年和1998—2017年水沙關(guān)系進(jìn)行變異診斷（一次診斷），見圖4。

從圖4中可知，兩時(shí)段水沙關(guān)系的一次診斷中統(tǒng)計(jì)量Zn最大值（41.174 4）出現(xiàn)在1990年，在顯著水平α為0.05時(shí)Zn遠(yuǎn)大于閾值（3.2），故存在變異點(diǎn)（1990年）。說明1958—1990年和1998—2017年兩時(shí)段水沙關(guān)系存在差異。

分別基于Clayton Copula函數(shù)和Frank Copula函數(shù)對1958—1990年和1998—2017年年徑流量和年輸沙量關(guān)系進(jìn)行變異診斷，見圖5。

從圖5中可知，在1958—1990年的一次診斷中統(tǒng)計(jì)量Zn最大值（5.933 7）出現(xiàn)在1965年，大于顯著水平α=0.05時(shí)閾值（3.2），故1965年為變異點(diǎn)。在二次診斷中統(tǒng)計(jì)量Zn最大值（1.695 2）小于閾值，故1966—1990年水沙關(guān)系不存在變異點(diǎn)，具有一定的平穩(wěn)性;在1998—2017年的一次診斷中統(tǒng)計(jì)量Zn最大值（0.203 9）小于閾值，因此認(rèn)為在1998—2017年水沙關(guān)系不存在變異點(diǎn)，具有一定的平穩(wěn)性。

頭道拐1998—2017年水沙邊緣分布的基于Frank Copula的聯(lián)合分布如圖8所示。從圖8（a）中可見，1998—2017年水沙規(guī)律性較強(qiáng)，服從固定分布，年徑流量和輸沙量概率密度函數(shù)呈現(xiàn)尾部相關(guān)的特點(diǎn)，上尾部和下尾部相關(guān)性最強(qiáng)且具有對稱性，表明年徑流量的極大或極小值處對年輸沙量有較大影響;從圖8（b）中可見，F(xiàn)rank Copula的分布函數(shù)具有對稱性，則水沙聯(lián)合分布的累計(jì)頻率具有對稱性。

3.4 計(jì)算遭遇頻率

根據(jù)建立的頭道拐站年徑流量和年輸沙量的P-Ⅲ型頻率分布，確定頭道拐站水沙豐平枯組合中年徑流量（X）和年輸沙量（Y）的劃分標(biāo)準(zhǔn)。選取1966—1990年和1998—2017年兩段時(shí)期計(jì)算，這兩段時(shí)期已被驗(yàn)證水沙關(guān)系存在一定的平穩(wěn)性，見表5～表8。

通過Clayton Copula建立水沙聯(lián)合分布計(jì)算頭道拐站1966—1990年水沙豐平枯組合遭遇頻率，見表9。

由表9可知，頭道拐站1966—1990年水沙豐平枯組合遭遇頻率中，水沙同頻組合（同豐、同偏豐、同平、同偏枯和同枯）頻率總和為0.645 6，水沙異頻組合頻率為0.354 4，同頻組合頻率明顯高于異頻組合;水沙異頻頻率的高頻部分主要集中于偏豐水平沙組合、平水偏豐沙組合、平水偏枯沙組合和偏枯水平沙組合，其頻率均大于0.04;總體而言，同偏豐組合（0.186 7）最大，同平組合（0.137 4）次之，水沙豐枯相反組合的頻率最小，說明頭道拐站1966—1990年年徑流量和年輸沙量頻率的關(guān)聯(lián)性強(qiáng)，且在年徑流量和年輸沙量偏豐時(shí)關(guān)聯(lián)性最強(qiáng)。

通過Frank Copula建立水沙聯(lián)合分布計(jì)算頭道拐站1998—2017年水沙豐平枯組合遭遇頻率，見表10。由表10可知，頭道拐站1998—2017年水沙豐平枯組合遭遇頻率中，同頻組合頻率總和為0.646 1，異頻組合頻率總和為0.353 9，則同頻組合頻率明顯高于異頻頻率;異頻頻率的高頻部分主要集中于豐水偏豐沙組合、偏豐水豐沙組合、偏豐水平沙組合、平水偏豐沙組合、平水偏枯沙組合、偏枯水平沙組合、偏枯水枯沙組合和枯水偏枯沙組合，其頻率均大于0.04;總體而言，同偏豐（同偏枯）組合（0.162 2）最大，同平組合（0.156 8）次之，水沙同豐（水沙同枯組）組合（0.082 5）第三，而枯水豐沙組合（豐水枯沙組合）最小，說明頭道拐站1998—2017年年徑流量和年輸沙量頻率的關(guān)聯(lián)性強(qiáng)，水沙狀態(tài)完全相反的組合的遭遇頻率極小。

4 結(jié) 論

以頭道拐站1958—1990年和1998—2017年兩個(gè)時(shí)期的年徑流量和年輸沙量序列為研究對象，利用P-Ⅲ頻率曲線探究水沙豐平枯組合遭遇頻率，基于Copula函數(shù)構(gòu)建二維水沙聯(lián)合分布，分析水沙聯(lián)合分布的累計(jì)頻率，得到以下結(jié)論。

（1）頭道拐站年徑流量和年輸沙量呈下降趨勢，在1990年水沙關(guān)系發(fā)生變異。1958—1990年水沙序列中存在一個(gè)變異點(diǎn)（1965年），1966—1990年和1998—2017年序列水沙關(guān)系未發(fā)生明顯變異，具有一定平穩(wěn)性。

（2）由基于Copula建立的水沙聯(lián)合分布的概率密度函數(shù)和分布函數(shù)可知，1958—1965年和1966—1990年水沙子序列概率密度函數(shù)上尾部相關(guān)性較強(qiáng)，即水沙序列中年徑流量的較大值處變化趨勢同年輸沙量變化趨勢關(guān)聯(lián)性較強(qiáng);1998—2017年水沙子序列的概率密度函數(shù)上尾部和下尾部相關(guān)性較強(qiáng)，即水沙序列中年徑流量的高低變化趨勢同年輸沙量高低變化趨勢關(guān)聯(lián)性較強(qiáng)。

（3）由水沙豐平枯組合遭遇頻率分析可得，在1966—1990年和1998—2017兩個(gè)時(shí)段的水沙組合中水沙同頻明顯高于水沙異頻，其中水沙異頻組合發(fā)生的概率較小，頭道拐站1966—1990年年徑流量和年輸沙量在平豐時(shí)相關(guān)性強(qiáng)，1998—2017年年徑流量和年輸沙量整體相關(guān)性較強(qiáng)，水沙豐平枯變化趨勢基本一致。

參考文獻(xiàn)：

[1] 胡春宏.黃河水沙變化與治理方略研究[J].水力發(fā)電學(xué)報(bào)，2016，35（10）：1-11.

[2] 趙陽，胡春宏，張曉明，等.近70年黃河流域水沙情勢及其成因分析[J].農(nóng)業(yè)工程學(xué)報(bào)，2018，34（21）：112-119.

[3] 高宗軍，馮國平.黃河水沙變化趨勢及成因分析[J].地下水，2020，42（1）：147-151.

[4] 高航，姚文藝，張曉華.黃河上中游近期水沙變化分析[J].華北水利水電學(xué)院學(xué)報(bào)，2009，30（5）：8-12.

[5] 金鑫，郝振純，張金良.黃河中游水沙頻率關(guān)系研究[J].泥沙研究，2006，31（3）：6-13.

[6] 李新杰，郜國明，朱亮，等.黃河潼關(guān)站水沙關(guān)系頻率及協(xié)調(diào)度分析研究[J].人民黃河，2020，42（5）：47-51.

[7] 馬雁，賈生海，張彥洪.黃河上游蘭州水文站2004—2018年水沙特性分析[J].水利規(guī)劃與設(shè)計(jì)，2020（4）：52-54，89.

[8] 李艷玲，暢建霞，黃強(qiáng)，等.基于滑動(dòng)Copula函數(shù)的降水和徑流關(guān)系變異診斷[J].水力發(fā)電學(xué)報(bào)，2014，33（6）：20-24，60.

[9] 郭愛軍，黃強(qiáng)，暢建霞，等.基于Copula函數(shù)的涇河流域水沙關(guān)系演變特征分析[J].自然資源學(xué)報(bào)，2015，30（4）：673-683.

[10] 姚曼飛，黨素珍，孟美麗，等.基于Copula函數(shù)的涇河流域水沙豐枯遭遇頻率分析[J].水土保持研究，2019，26（1）：192-196，202.

[11] 冉大川，姚文藝，張攀，等.黃河頭道拐站水沙來源空間分布及其影響因素[J].泥沙研究，2015，40（1）：42-48.

[12] 任智慧，王婷，曲少軍.萬家寨水庫庫區(qū)沖淤特點(diǎn)分析[C]//中國大壩工程學(xué)會(huì).水庫大壩高質(zhì)量建設(shè)與綠色發(fā)展中國大壩工程學(xué)會(huì)2018學(xué)術(shù)年會(huì)論文集.北京：中國大壩工程學(xué)會(huì)，2018：172-177.

[13] 雷冠軍，王文川，殷峻暹，等.P-Ⅲ型曲線參數(shù)估計(jì)方法研究綜述[J].人民黃河，2017，39（10）：1-7.

[14] 黃繼文.P-Ⅲ型分布頻率分析在Excel中的實(shí)現(xiàn)及應(yīng)用[J].水資源研究，2006，27（4）：7-9.

[15] 郭生練，葉守澤.論水文計(jì)算中的經(jīng)驗(yàn)頻率公式[J].武漢水利電力學(xué)院學(xué)報(bào)，1992，25（2）：38-45.

[16] PATTON A J. A Review of Copula Models for Economic Time Series[J].Journal of Multivariate Analysis，2012，110：4-18.

[17] 杜懿，麻榮永.不同Copula函數(shù)在洪水峰量聯(lián)合分布中的應(yīng)用比較[J].水力發(fā)電，2018，44（12）：24-26，58.

[18] 李天元，郭生練，羅啟華，等.雙參數(shù)Copula函數(shù)在洪水聯(lián)合分布中的應(yīng)用研究[J].水文，2011，31（5）：24-28，46.

[19] 賈玉紅，宋松柏.陜北地區(qū)年降水量頻率分布參數(shù)估算研究[J].水資源與水工程學(xué)報(bào)，2012，23（5）：48-50.

[20] 宋喜芳，李建平，胡希遠(yuǎn).模型選擇信息量準(zhǔn)則AIC及其在方差分析中的應(yīng)用[J].西北農(nóng)林科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2009，37（2）：88-92.

[21] HUANG S Z， LI P， HUANG Q， et al. Copula-Based Identification of the Non-Stationarity of the Relation Between Runoff and Sediment Load[J].International Journal of Sediment Research， 2017，32（2）：221-230.

[22] COSTA D A D. Copula Inference for Finance and Insurance[D].Zurich，Switzer-Land：ETH，2004：123-141.

[23] 中華人民共和國水利部.水文基本術(shù)語和符號標(biāo)準(zhǔn)：GB/T 50095—2014[S].北京：中國計(jì)劃出版社，2014：2-51.

[24] 徐宇程，朱首賢，張文靜，等.長江大通站徑流量的豐平枯水年劃分探討[J].長江科學(xué)院院報(bào)，2018，35（6）：19-23.

[25] 林沫，劉穎，叢遠(yuǎn)飛，等.松遼流域主要河流水沙規(guī)律分析[J].東北水利水電，2009，27（12）：40-42，72.

【責(zé)任編輯 張 帥】