王新棟,漆文邦,張江潮

(四川大學(xué),成都 610065)

0 前 言

隨著大壩觀測(cè)技術(shù)發(fā)展的深入、大壩測(cè)量設(shè)備精度的提高、大壩監(jiān)測(cè)理論的逐漸成熟，在大壩監(jiān)測(cè)原始觀測(cè)數(shù)據(jù)的基礎(chǔ)上，利用反演分析來(lái)推測(cè)壩體混凝土及基巖材料參數(shù)的方法業(yè)已成熟。反演分析法的原理基于現(xiàn)場(chǎng)測(cè)量的位移、應(yīng)力、應(yīng)變等數(shù)據(jù)，通過(guò)反演模型推出一些參數(shù)(比如彈模、泊松比等)的方法[1]。其最終結(jié)果是基于大壩實(shí)測(cè)數(shù)據(jù)，建立理論預(yù)測(cè)模型擬合大壩的各種特征，以此來(lái)反饋控制大壩的安全。

位移量是描述物體受力變形形態(tài)的一類(lèi)重要物理量，目前已被用作反演分析確定巖土工程材料特性參數(shù)的主要依據(jù)[2-3]。

1 工程概況

某水電站是一座以發(fā)電為單一目標(biāo)、具有日調(diào)節(jié)水庫(kù)的高水頭引水式電站。電站總裝機(jī)52 MW，設(shè)計(jì)水頭233.00 m，最大水頭272.00 m，最小水頭233.00 m，加權(quán)平均水頭258.00 m，最大引用流量24 m3/s，保證出力11.57 MW，年利用小時(shí)數(shù)5 110 h，年發(fā)電量2.45億kWh。電站水庫(kù)正常蓄水位1 090.00 m，設(shè)計(jì)洪水位1 084.60 m，校核洪水位1 088.85 m，電站正常蓄水位以下庫(kù)容26萬(wàn)m3，調(diào)節(jié)庫(kù)容21.2萬(wàn)m3，水庫(kù)具有日調(diào)節(jié)能力。本工程屬Ⅳ等工程，永久建筑物按4級(jí)次要建筑物按5級(jí)設(shè)計(jì)。相應(yīng)防洪標(biāo)準(zhǔn)，首部樞紐永久性建筑物按50年一遇洪水設(shè)計(jì)、200年一遇洪水校核。

某水電站首部建有一座有長(zhǎng)126.6 m、最大閘高31 m的攔河閘壩日調(diào)節(jié)水庫(kù)。首部樞紐從左至右依次布置有左岸錨索邊坡、進(jìn)水閘、泄洪閘(1、2、3號(hào)閘)、表孔閘(4號(hào)閘)、右岸擋水壩段及右岸護(hù)岸工程。進(jìn)水閘側(cè)向布置在左岸，其軸線與攔河閘壩軸線夾角為50°。閘室前為設(shè)有靜水柵的靜水段和沉礫段，此段與隧洞進(jìn)水口之間設(shè)高程為1 077.50 m的攔沙砍。隧洞進(jìn)水閘前設(shè)有攔污柵，攔污柵底板下布置3孔沖沙孔，孔底高程1 074.00 m，3孔沖沙孔匯集于后面的1.0 m×1.5 m(寬×高)沖砂廊道。3孔泄洪沖砂閘(1、2、3號(hào)閘)緊靠進(jìn)水閘右側(cè)布置，為平底板閘，孔口尺寸7.0 m×8.0 m，每孔設(shè)弧形工作閘門(mén)各1道，共用1道平板檢修門(mén)。閘底板高程1 070.50 m，閘頂高程1 091.50 m，閘室長(zhǎng)39 m，墩厚分別為2.60 m(邊墩，與進(jìn)水閘分縫)、3.6 m、3.2 m(2號(hào)與3號(hào)閘墩，中間縫寬2 cm)、3.2 m(3號(hào)閘與表孔閘中間縫寬2 cm)。1、2、3號(hào)閘基礎(chǔ)建基面1 060.50 m。表孔閘(4號(hào)閘)緊靠3孔泄洪沖砂閘布置，閘孔寬7.0 m，閘室長(zhǎng)39 m。堰頂采用WES曲線，堰頂高程1 083.00 m，溢流面采用400號(hào)硅粉混凝土防沖，4號(hào)閘基礎(chǔ)建基面1 060.50 m。

2 反演分析方法

水利及巖土工程的反演分析可以分為以下幾類(lèi)：位移反分析法、應(yīng)力反分析法以及二者結(jié)合綜合反分析法。位移反分析法是目前應(yīng)用最為普遍的方法，它又分為解析法及數(shù)值法。解析法更多的適用于求解較為簡(jiǎn)單邊界條件下的線粘彈性問(wèn)題，在解決這類(lèi)問(wèn)題時(shí)解析法優(yōu)點(diǎn)突出，主要表現(xiàn)在概念明確、計(jì)算速度快。相比之下，數(shù)值法解決范圍更加廣泛，且主要用在解決復(fù)雜邊界及復(fù)雜地質(zhì)條件下的非線性問(wèn)題。根據(jù)數(shù)值法的求解過(guò)程不同，可以分為正反分析法、逆反分析法等[4-5]。

正反分析法指的是首先指定待反演參數(shù)的初始值，通過(guò)大型計(jì)算軟件進(jìn)行正分析得到典型測(cè)點(diǎn)的變位或者應(yīng)力，之后將計(jì)算值與實(shí)測(cè)值進(jìn)行擬合比較，當(dāng)誤差較大時(shí)通過(guò)有效方式調(diào)整待反演參數(shù)，再重復(fù)以上操作直至計(jì)算值與實(shí)測(cè)值之間的誤差達(dá)到允許的范圍內(nèi)，此時(shí)取定的待反演參數(shù)值即作為反演計(jì)算的最終結(jié)果。正反分析方法程序編制簡(jiǎn)單，計(jì)算方法靈活、適用非線性問(wèn)題的參數(shù)反演[6-7]。

逆反分析法與正反分析法的計(jì)算程序相反，它是指根據(jù)待反演參數(shù)與實(shí)測(cè)值之間的某種關(guān)系，直接求解出由該關(guān)系式組成的方程組，即可得到反演參數(shù)。該方法原理簡(jiǎn)單直觀，且步驟簡(jiǎn)單，計(jì)算用時(shí)短，但是需要進(jìn)行程序編制，且更多情況下適用于線性問(wèn)題的反演[8]。

考慮到某水電站大壩已進(jìn)入運(yùn)行期，經(jīng)模型有效性初步驗(yàn)證，多元線性回歸及逐步回歸效果良好，因此本文采用基于多重響應(yīng)面的位移正反分析法來(lái)進(jìn)行材料參數(shù)的反演分析。

基于響應(yīng)面方法的原理，提出多重響應(yīng)面方法來(lái)擬合壩體參數(shù)與大壩在荷載作用下產(chǎn)生的順河向位移之間的非線性映射關(guān)系，相比于其他反分析方法存在的大量的數(shù)值計(jì)算過(guò)程，該方法可以有效減少數(shù)值計(jì)算次數(shù)，節(jié)省反分析時(shí)間，提高計(jì)算效率。

為了保證計(jì)算效率和準(zhǔn)確度，在本文中，運(yùn)用不帶交叉項(xiàng)的 2 次響應(yīng)面函數(shù)，其表達(dá)式為：

(1)

式中：x=(x1,x2,…,xn)為待反分析矢量；a、bi、ci(i=1,2,…,n)為待定系數(shù)。式(1)中右端共有2n+1個(gè)待定系數(shù)，需要有2n+1個(gè)方程才能求解。在本文的反分析中s(x)代表位移，x代表待反演參數(shù)。

采用響應(yīng)面法對(duì)大壩的各個(gè)參數(shù)進(jìn)行反分析，將計(jì)算軟件內(nèi)部的計(jì)算過(guò)程用響應(yīng)面方程來(lái)近似，可大大減少軟件的計(jì)算次數(shù)，提高反分析效率。

3 有限元模型

3.1 變形監(jiān)測(cè)

2013年9月，根據(jù)電站監(jiān)測(cè)設(shè)施施工設(shè)計(jì)圖及電站現(xiàn)場(chǎng)實(shí)際情況，在右岸布置二等單三角控制網(wǎng)，其中BL1、BL3作為日常水平位移監(jiān)測(cè)的工作基點(diǎn)，BL2作為檢查點(diǎn)。高程控制點(diǎn)，需由原監(jiān)測(cè)點(diǎn)聯(lián)測(cè)至BL2和BL3，BL2作為垂直位移監(jiān)測(cè)的工作基點(diǎn)，BL3作為檢查點(diǎn)。垂直位移監(jiān)測(cè)均采用二等水準(zhǔn)進(jìn)行測(cè)量。

垂直位移點(diǎn)和水平位移點(diǎn)均布設(shè)在同一位置，參照原設(shè)計(jì)的點(diǎn)位布設(shè)，右岸擋水壩測(cè)點(diǎn)編號(hào)：LD1，壩頂上游側(cè)測(cè)點(diǎn)編號(hào)：LD2、LD3，進(jìn)水口的閘墩測(cè)點(diǎn)編號(hào)：LD4，下游閘墩測(cè)點(diǎn)編號(hào)：LD5、LD6。首部樞紐變形監(jiān)測(cè)點(diǎn)布置見(jiàn)圖1。

圖1 首部樞紐變形監(jiān)測(cè)點(diǎn)布置圖

3.2 有限元模型的建立

電站首部樞紐從左至右依次布置有左岸錨索邊坡、進(jìn)水閘(側(cè)向布置在左岸)、泄洪閘(1、2、3號(hào)閘)、表孔閘(4號(hào)閘)、右岸擋水壩段及右岸護(hù)岸工程。本文采用基于多重響應(yīng)面的位移正反分析法來(lái)對(duì)電站大壩泄洪閘段某些參數(shù)進(jìn)行反演，分析軟件采用大型通用有限元軟件ANSYS來(lái)進(jìn)行分析。并將反演得到的參數(shù)進(jìn)行正分析運(yùn)算，通過(guò)應(yīng)力及變形結(jié)果來(lái)驗(yàn)證反演參數(shù)的合理性。

模型采用D-P準(zhǔn)則，壩體與基礎(chǔ)網(wǎng)格采用八節(jié)點(diǎn)六面體SOLID45單元。

模型采用計(jì)算坐標(biāo)系如下：

X軸,順河流方向，指向下游為正方向。

Y軸,鉛垂方向，鉛垂向上為正方向。

Z軸,垂直于水流方向，指向右岸為正方向。

模型范圍：

(1) 泄洪閘。泄洪閘壩段由3個(gè)閘室組成，其中靠近右岸的2個(gè)閘室連為一體，與靠近左岸的閘室用沉降縫分離。泄洪閘三維計(jì)算模型僅選取連為一體的2個(gè)閘室，閘高31 m。模型邊界確定為鉛垂方向取至閘室底邊界以下62 m；順河向上游取至閘室上游邊界以上50 m，向下游取至閘室下游邊界以下50 m；垂直水流方向以閘室左右邊界起向左右各延伸50 m。模型共劃分51 582個(gè)單元，56 210個(gè)節(jié)點(diǎn)。

(2) 表孔閘。表孔閘閘高31 m。模型邊界確定為：鉛垂方向取至閘室底邊界以下62 m；順河向上游取至閘室上游邊界以上50 m，向下游取至閘室下游邊界以下50 m；垂直水流方向以閘室左右邊界起向左右各延伸50 m。模型共劃分36 156個(gè)單元，39 710個(gè)節(jié)點(diǎn)。

模型邊界基礎(chǔ)四周和底面采用法向約束，結(jié)構(gòu)頂部及四周采用自由邊界。

4 參數(shù)反演分析計(jì)算

電站泄洪閘壩段及地基共7種材料，每種材料有5個(gè)參數(shù)，如果全部反演，則工程量十分浩大。從大壩建成到現(xiàn)在已經(jīng)運(yùn)行了17 a，隨著時(shí)間的推移，混凝土的彈性模量也在逐步的變化，因此本文針對(duì)影響大壩變形的6個(gè)主要參數(shù)進(jìn)行反演分析，依次是基巖t1彈模、基巖t2彈模、基巖t1泊松比、基巖t2泊松比、C10彈模、C20彈模。反演過(guò)程中通過(guò)2個(gè)工況的位移計(jì)算值差值與實(shí)測(cè)差值的擬合度來(lái)判定反演參數(shù)是否合理。

溫度因子是影響混凝土大壩變形的最主要因素，為了盡量減小誤差及消除溫度場(chǎng)對(duì)位移造成的影響，本文選取溫差最小、上游水位差值最大且受時(shí)效因素影響最小的2015年1月24日(溫度6 ℃，上游水位1 084.00 m，下游水位1 065.60 m，定義為工況1)、2015年2月25日(溫度6.5 ℃，上游水位1 088.51 m，下游水位1 065.60 m，定義為工況2)2種工況下的實(shí)測(cè)數(shù)據(jù)來(lái)進(jìn)行參數(shù)的反演。反演過(guò)程中通過(guò)工況2與工況1的位移計(jì)算值差值y2′-y1′與實(shí)測(cè)差值y2-y1的擬合度來(lái)判定反演參數(shù)是否合理。根據(jù)電站水庫(kù)運(yùn)行方式可知：首部樞紐泄洪沖沙閘及表孔溢流壩在枯期(12月-翌年3月)均關(guān)閉，可保持水位在1 090.00 m運(yùn)行。若遇水位超過(guò)1 090.00 m時(shí)可開(kāi)啟表孔閘門(mén)進(jìn)行調(diào)節(jié)，以保證水庫(kù)運(yùn)行水位在1 081.50～1 090.00 m內(nèi)運(yùn)行。而2種工況下的水位上游水位均小于1 090.00 m 。綜上可知2種工況下閘門(mén)全部關(guān)閉。以此來(lái)作為荷載施加的依據(jù)。

(1) 計(jì)算中選用不帶交叉項(xiàng)的二次多項(xiàng)式作為響應(yīng)面函數(shù)，共有6 個(gè)變量，所以響應(yīng)面函數(shù)有 13個(gè)未知量。根據(jù)土工試驗(yàn)及地質(zhì)報(bào)告得到的各個(gè)物理力學(xué)參數(shù)的建議值，通過(guò)SPSS軟件進(jìn)行正交試驗(yàn)得到13組參數(shù)組[9]。

(2) 分別取由第(1)步得到的13 組參數(shù)組里面的每一組參數(shù)，用ANSYS軟件做 13次正分析計(jì)算，得到13組泄洪閘壩段的LD2、LD3、LD5、LD6測(cè)點(diǎn)順河向及垂直位移數(shù)據(jù)。

(3) 對(duì)于LD2監(jiān)測(cè)點(diǎn)，將第一組參數(shù)值與位移值代入式(1)得到含有13個(gè)未知量的一次方程。以此類(lèi)推，由另外 12 組參數(shù)值及位移值可以得到相應(yīng)的12個(gè)方程，總共可以得到由13個(gè)線性方程組成的方程組。經(jīng)求解得到在LD2測(cè)點(diǎn)處建立的響應(yīng)面函數(shù)系數(shù)。同樣的，對(duì)LD3、LD5、LD6監(jiān)測(cè)點(diǎn)，采用同樣的方法，得到其對(duì)應(yīng)的響應(yīng)面函數(shù)的二項(xiàng)式系數(shù)。

將響應(yīng)面函數(shù)的系數(shù)值帶入式(1)，可得LD2測(cè)點(diǎn)的響應(yīng)面函數(shù)。同樣方法可得到LD3、LD5、LD6測(cè)點(diǎn)的響應(yīng)面函數(shù)。

(4) 建立8個(gè)響應(yīng)面函數(shù)方程后，需要找到 1 組參數(shù)值，使得由這8個(gè)響應(yīng)面得到的2種工況下的計(jì)算位移差值與實(shí)測(cè)差值之間的差值的平方和最小，即使得目標(biāo)函數(shù)最小。本文采取matlab中的優(yōu)化工具箱來(lái)進(jìn)行多參數(shù)的優(yōu)化求解，得到的6個(gè)待反演參數(shù)值如表1所示。

η= [S2H(x)-(0.39)]2+[S2V(x)-(0.2)]2+

[S3H(x)-(1.33)]2+[S3V(x)-(0.1)]2+

[S5H(x)-(0.27)]2+[S5V(x)-(0.4)]2+

[S6H(x)-(0.32)]2+[S6V(x)-(0.2)]2

(2)

表1 反演參數(shù)數(shù)據(jù)表

將6個(gè)反演參數(shù)值在2種工況下通過(guò)ANSYS軟件進(jìn)行正分析，得到LD2、LD3、LD5、LD6測(cè)點(diǎn)的順河向及垂直位移數(shù)據(jù)計(jì)算值的差值與實(shí)測(cè)值位移數(shù)據(jù)差值，如表2所示。

可以看出：4個(gè)測(cè)點(diǎn)中LD2、LD5、LD6的順河向位移計(jì)算差值與實(shí)測(cè)位移差值都比較接近，只有LD3相差較大。4個(gè)測(cè)點(diǎn)中LD5、LD6垂直位移計(jì)算差值與實(shí)測(cè)位移差值較接近，LD2、LD3的垂直位移計(jì)算差值的符號(hào)為負(fù)，而實(shí)測(cè)垂直位移差值符號(hào)為正。規(guī)定向下沉降為正，向上抬升為負(fù)。壩體受到指向下游的推力大于指向上游的推力；同時(shí)壩體自身結(jié)構(gòu)分布均勻，底板所受的揚(yáng)壓力從上游往下游逐漸減小。兩者共同作用導(dǎo)致壩體所受順時(shí)針?lè)较驈澗?。而工況2的上游水位比工況1高4.5 m左右，2種工況的下游水位相等，工況2大壩產(chǎn)生的向順時(shí)針?lè)较虻膹澗刈饔酶?，因此工況2位于閘墩上游側(cè)的LD2、LD3測(cè)點(diǎn)垂直位移較工況1小，位于閘墩下游側(cè)的LD5、LD6測(cè)點(diǎn)垂直位移較工況1大，LD2、LD3垂直位移2種工況下計(jì)算值的差值為負(fù)，LD5、LD6垂直位移2種工況況下計(jì)算值的差值為正是合理的。綜上可知，小部分實(shí)測(cè)數(shù)據(jù)在觀測(cè)過(guò)程中可能存在一定的誤差，大部分實(shí)測(cè)數(shù)據(jù)差值與位移計(jì)算差值接近，反演參數(shù)基本合理。泄洪閘位移云圖如圖2、3所示。

表2 計(jì)算值差值與實(shí)測(cè)值差值數(shù)據(jù)對(duì)比表

圖2 泄洪閘順河向位移云圖

5 典型工況下大壩變形及應(yīng)力分析

5.1 判定標(biāo)準(zhǔn)

(1) 沉降：根據(jù)SL265-2001《水閘設(shè)計(jì)規(guī)范》規(guī)定，天然土質(zhì)地基上水閘地基最大沉降量不宜超過(guò)15 cm，相鄰部位最大沉降不宜超過(guò)5 cm。

(2) 應(yīng)力：依照GB50010-2002《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》規(guī)定，混凝土的抗壓抗拉強(qiáng)度滿足混凝土強(qiáng)度設(shè)計(jì)值的要求，計(jì)算所用到混凝土強(qiáng)度設(shè)計(jì)值：C20，抗壓9.6 MPa、抗拉1.10 MPa；C40，抗壓19.1 MPa、抗拉1.71 MPa。

圖3 泄洪閘垂直位移云圖

5.2 選取的計(jì)算工況

(1) 正常蓄水位工況：上游1 090.00 m，下游無(wú)水。

(2) 設(shè)計(jì)洪水位工況：上游水位1 084.60 m，下游水位1 073.80 m。

(3) 校核洪水位工況：上游水位1 088.85 m，下游水位1 072.20 m。

利用反演出的參數(shù)，對(duì)泄洪沖沙閘及表孔閘在正常蓄水位工況、設(shè)計(jì)洪水位工況、校核洪水位工況3種工況下進(jìn)行了安全性論證，結(jié)果表明：

(1) 3種工況下泄洪閘的位移值均較小，且對(duì)比可知正常蓄水位工況下泄洪閘的位移值最大，最大順河向位移產(chǎn)生于閘墩頂部，位移值為1.46 mm，最大豎向位移發(fā)生于閘墩下游頂部，位移值為-2.75 mm。表孔閘在對(duì)應(yīng)工況下變形規(guī)律與泄洪閘類(lèi)似，在正常蓄水位工況下表孔閘的位移值最大，最大順河向位移產(chǎn)生于閘墩頂部，位移值為0.75 mm，最大豎向位移發(fā)生于閘墩下游頂部，位移值為-2.6 mm。兩壩段在3種典型工況下變形均在規(guī)范允許范圍之內(nèi)，是安全的。

(2) 3種工況下將大小主應(yīng)力進(jìn)行對(duì)比發(fā)現(xiàn)，泄洪閘在正常蓄水位工況下應(yīng)力最大，最大壓應(yīng)力為-1.63 MPa，出現(xiàn)在閘墩與閘底板相接處，最大拉應(yīng)力為1.99 MPa，出現(xiàn)在牛腿與邊墩的接觸面上。表孔閘在校核工況下應(yīng)力最大，最大壓應(yīng)力為-1.4 MPa，出現(xiàn)在閘墩與閘底板相接處，最大拉應(yīng)力為0.233 MPa，出現(xiàn)在溢流面處及閘室下游的末端位置。兩壩段在3種典型工況下應(yīng)力均在規(guī)范允許范圍之內(nèi)，是安全的。

綜上所述，反演參數(shù)基本合理且安全，反演方法正確可行。

6 結(jié) 語(yǔ)

本文通過(guò)分析研究多種反演分析方法，結(jié)合電站監(jiān)測(cè)資料，采用非線性問(wèn)題的位移正反分析方法。位移正反分析法基于正演分析，方法靈活，程序簡(jiǎn)單，適用于電站大壩參數(shù)的反演。通過(guò)對(duì)某電站大壩泄洪沖沙閘及表孔閘進(jìn)行了有限元模擬，采用基于多重響應(yīng)面的位移正反分析法，對(duì)影響大壩位移的基巖t1彈模、基巖t2彈模、基巖t1泊松比、基巖t2泊松比、C10彈模、C20彈模進(jìn)行了反演分析。結(jié)果表明：計(jì)算值與實(shí)測(cè)值中大部分值擬合的較接近，參數(shù)反演基本合理且參數(shù)反演方法正確可行。同時(shí)利用反演分析得到的參數(shù)通過(guò)對(duì)泄洪沖沙閘及表孔閘2個(gè)閘壩段在正常蓄水位、設(shè)計(jì)洪水位、校核洪水位3種工況下進(jìn)行了變形及應(yīng)力分析。位移分布符合實(shí)測(cè)規(guī)律，位移值及應(yīng)力值在規(guī)范允許范圍內(nèi)，參考類(lèi)似工程，反演參數(shù)值合理[10]。

參考文獻(xiàn):

[1] 周萍.基于神經(jīng)網(wǎng)絡(luò)和遺傳算法的巖體參數(shù)反分析研究[D].南京：河海大學(xué),2004.

[2] 吳瓊.隧道動(dòng)態(tài)設(shè)計(jì)中圍巖計(jì)算參數(shù)的應(yīng)力反分析研究[D].成都：西南交通大學(xué),2006.

[3] 周麗君.時(shí)空效應(yīng)下巖質(zhì)邊坡位移反分析研究[D].重慶：重慶大學(xué),2009.

[4] 徐韜.某巖質(zhì)高邊坡穩(wěn)定性計(jì)算與力學(xué)參數(shù)反分析[D].長(zhǎng)沙：中南大學(xué),2013.

[5] 裴建春.混凝土面板堆石壩原型觀測(cè)資料反饋分析研究[D].天津：天津大學(xué),2007.

[6] 聞世強(qiáng).茅坪溪瀝青混凝土心墻堆石壩反演分析[D].南京：河海大學(xué),2004.

[7] 班宏泰,覃文文,趙昕,柏洪雪,謝鵬,費(fèi)文平.巖體參數(shù)的反演分析方法[J].科技創(chuàng)新導(dǎo)報(bào),2008(08):7-8.

[8] 唐光華.圍巖位移反分析及智能化灰色預(yù)測(cè)與模糊分類(lèi)的應(yīng)用[D].武漢：武漢理工大學(xué),2005.

[9] 寧昕揚(yáng),劉曉青,牛志偉,齊慧君,潘世洋.基于響應(yīng)面法的土石壩參數(shù)隨機(jī)反演分析[J].水電能源科學(xué),2017(01):73-76.

[10] 雷曉鋒,馬雙科,鄭婭娜.某水電工程地下廠房開(kāi)挖的巖體參數(shù)反分析[J].中國(guó)農(nóng)村水利水電,2017(07):140-143.