邊坡位移預(yù)警指標(biāo)的實(shí)時(shí)估計(jì)與診斷

任　杰，蘇懷智，楊　孟，周志杰

(1. 河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇南京　210098； 2. 河海大學(xué)水利水電學(xué)院，江蘇南京　210098)

?

邊坡位移預(yù)警指標(biāo)的實(shí)時(shí)估計(jì)與診斷

任杰1,2，蘇懷智1,2，楊孟1,2，周志杰2

(1. 河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇南京210098； 2. 河海大學(xué)水利水電學(xué)院，江蘇南京210098)

摘要：邊坡預(yù)警指標(biāo)的實(shí)時(shí)估計(jì)與診斷是實(shí)現(xiàn)邊坡安全監(jiān)控、預(yù)防邊坡演變?yōu)榛碌闹匾侄?。安全預(yù)警注重考慮極端事件，其關(guān)鍵在于能否準(zhǔn)確進(jìn)行效應(yīng)量序列分位數(shù)分析及刻畫(huà)其尾部特征?；跇O值理論中超閾值(Peaks over Threshold, POT)模型，分時(shí)段考慮邊坡位移監(jiān)測(cè)序列，借助Hill圖法擬定合理的閾值，再利用廣義帕累托分布(Generalized Pareto Distribution, GPD)對(duì)超閾值序列進(jìn)行擬合分析，漸進(jìn)地刻畫(huà)位移序列分布的尾部特征，分析其抵御已經(jīng)歷荷載的能力。假定邊坡失事概率，得到位移預(yù)警指標(biāo)實(shí)時(shí)估計(jì)序列，挖掘或評(píng)估出該邊坡面臨可能發(fā)生極端荷載時(shí)的抵御能力。最后基于突變理論中尖點(diǎn)突變模型對(duì)位移預(yù)警指標(biāo)實(shí)時(shí)估計(jì)序列進(jìn)行突變?cè)\斷，結(jié)合具體邊坡工程實(shí)例，驗(yàn)證了POT-突變理論混合模型的合理性。

關(guān)鍵詞：邊坡； 位移預(yù)警指標(biāo)； 極值理論； POT模型； 尖點(diǎn)突變理論

根據(jù)邊坡工程性態(tài)及原始監(jiān)測(cè)資料，借助先進(jìn)的數(shù)學(xué)、力學(xué)工具建立各種效應(yīng)量的預(yù)警指標(biāo)估計(jì)模型，是實(shí)現(xiàn)邊坡工程安全預(yù)警的重要手段。邊坡工程安全預(yù)警的主要目的在于通過(guò)實(shí)時(shí)監(jiān)測(cè)并分析邊坡變形及發(fā)展態(tài)勢(shì)，以便能及時(shí)采取工程或非工程措施，達(dá)到防止或減少邊坡破壞所帶來(lái)的災(zāi)害損失。對(duì)邊坡效應(yīng)量預(yù)警指標(biāo)的估計(jì)是實(shí)現(xiàn)邊坡安全預(yù)警的前提[1]。預(yù)警機(jī)理在于對(duì)抵御已經(jīng)歷荷載能力的分析，評(píng)估和預(yù)測(cè)抵御可能發(fā)生荷載的能力。預(yù)警指標(biāo)估計(jì)重點(diǎn)考慮極端事件，借助極值理論可實(shí)現(xiàn)對(duì)其的估計(jì)和研究[2]。

本文基于極值理論中POT模型進(jìn)行建模。其流程通常是依據(jù)一定的規(guī)則擬定合理的閾值，以超閾值序列作為子樣本序列，利用GPD分布擬合子樣本序列，得到其分布函數(shù)，根據(jù)統(tǒng)計(jì)學(xué)條件概率公式繼而得到母樣本的分布函數(shù)，結(jié)合邊坡失事概率完成對(duì)預(yù)警指標(biāo)的估計(jì)。A. A. Balkema以及J. Pickands等指出，來(lái)自同一總體的簡(jiǎn)單隨機(jī)樣本，若閾值擬定充分大，超閾值的隨機(jī)變量的極限分布即為GPD分布[3-4]。大量學(xué)者研究表明，基于POT模型的極值理論方法在擴(kuò)大樣本容量，提高樣本質(zhì)量的同時(shí)比其他方法能夠更準(zhǔn)確地描述序列分布的尾部特征，尤其在采樣數(shù)據(jù)較少的情況下，是一種較為準(zhǔn)確的分位數(shù)分析和預(yù)測(cè)手段[2,5-7]。

通過(guò)對(duì)不同時(shí)段位移監(jiān)測(cè)序列進(jìn)行POT建模分析，得到預(yù)警指標(biāo)實(shí)時(shí)估計(jì)序列，再利用尖點(diǎn)突變理論對(duì)其變化趨勢(shì)進(jìn)行性態(tài)突變?cè)\斷，結(jié)合某邊坡工程實(shí)例，對(duì)該法的合理性進(jìn)行了驗(yàn)證。

1邊坡位移預(yù)警指標(biāo)估計(jì)極值POT模型

1.1GPD分布與POT模型

利用GPD分布擬合超閾值序列來(lái)漸進(jìn)地刻畫(huà)分布的尾部特征是建立POT模型的基礎(chǔ)。GPD分布定義為：

(1)

式中：ε為形狀參數(shù)；當(dāng)ε≥0時(shí)，x≥0；當(dāng)ε<0時(shí)，u

1.2超閾值分布函數(shù)(CEDF)和總體分布函數(shù)

(2)

式中：0≤y

根據(jù)條件概率公式，可得總體分布函數(shù)：

(3)

(4)

(5)

1.3位移監(jiān)測(cè)序列厚尾檢驗(yàn)

1.4閾值u擬定的Hill圖法

利用GPD分布擬合超閾值序列，需解決閾值u擬定這個(gè)關(guān)鍵問(wèn)題，同時(shí)閾值u也是準(zhǔn)確估計(jì)參數(shù)ε和σ的前提。若閾值u的擬定過(guò)高，超閾值序列樣本容量少，使得參數(shù)估計(jì)的方差過(guò)大；若閾值u的擬定過(guò)低，則有可能不滿(mǎn)足PBDH定理成立的條件，使得超閾值序列樣本分布與GPD分布產(chǎn)生較大差異。閾值擬定方法分作圖法與數(shù)值計(jì)算法兩類(lèi)[6,9]。本文借助Hill圖法進(jìn)行閾值u的擬定。

(6)

1.5GPD分布參數(shù)估計(jì)

GPD分布參數(shù)估計(jì)的方法有極大似然估計(jì)法、概率權(quán)矩法以及L矩估計(jì)法等。本文采用最常用的極大似然估計(jì)法對(duì)GPD分布參數(shù)進(jìn)行估計(jì)。

(7)

(8)

(9)

至此，根據(jù)上文公式及計(jì)算所得參數(shù)便可確定式(5)的分布函數(shù)。據(jù)分布函數(shù)反函數(shù)，結(jié)合邊坡失事概率，可以求得邊坡位移預(yù)警指標(biāo)估計(jì)值。

2邊坡位移預(yù)警指標(biāo)性態(tài)診斷突變理論

2.1尖點(diǎn)突變模型

突變理論是以奇點(diǎn)理論、穩(wěn)定性理論等數(shù)學(xué)理論為基礎(chǔ)，用以研究微積分方法不能解釋的不連續(xù)變化現(xiàn)象的理論[11-12]。突變理論通過(guò)給出系統(tǒng)在突變過(guò)程的勢(shì)函數(shù)，研究其突變流形和分叉集進(jìn)而確定系統(tǒng)所屬狀態(tài)。根據(jù)邊坡失穩(wěn)特征，采用尖點(diǎn)突變模型判別其穩(wěn)定性[13-14]。

(10)

式中：x為邊坡?tīng)顟B(tài)變量；u，v為邊坡控制變量。

(11)

突變流形的實(shí)根個(gè)數(shù)由下式Δ符號(hào)判斷：Δ=8u3+27v2

(12)

Δ>0時(shí)有一個(gè)實(shí)根，邊坡處于穩(wěn)定狀態(tài)；Δ<0時(shí)有3個(gè)互異實(shí)根，邊坡處于非穩(wěn)定狀態(tài)；Δ=0時(shí)，邊坡處于臨界狀態(tài)。

2.2邊坡預(yù)警指標(biāo)尖點(diǎn)突變模型

本文采用經(jīng)驗(yàn)方法[11,13]建立尖點(diǎn)突變模型。對(duì)于邊坡預(yù)警指標(biāo)估計(jì)序列，采用最小二乘法進(jìn)行數(shù)據(jù)擬合，得到突變流形結(jié)果如下式所示：

(13)

(14)

(15)

計(jì)算所得u和v，代入式(12)中可計(jì)算Δ，進(jìn)而判別預(yù)警指標(biāo)序列的穩(wěn)定狀態(tài)。分析過(guò)程中，往往需確定突變點(diǎn)位置，本文借鑒文獻(xiàn)[15]的方式確定出在何時(shí)發(fā)生了突變。年份為ti時(shí)，取前i年的預(yù)警指標(biāo)估計(jì)序列進(jìn)行式(13)～(15)計(jì)算，并獲取Δi，當(dāng)Δi=0時(shí)，所對(duì)應(yīng)年份即為突變發(fā)生臨界年份，若前i年數(shù)不足以進(jìn)行式(13)擬合，則通過(guò)線(xiàn)性插值使點(diǎn)數(shù)增加。

3邊坡工程實(shí)例

為闡述建模方法，驗(yàn)證模型可靠性，本文以黃河流域某邊坡工程為例，依據(jù)該邊坡某測(cè)點(diǎn)2009—2013年垂直方向位移監(jiān)測(cè)數(shù)據(jù)，以年為區(qū)間，獲得2009年73個(gè)監(jiān)測(cè)數(shù)據(jù)，2010年170個(gè)監(jiān)測(cè)數(shù)據(jù)，2011年165個(gè)監(jiān)測(cè)數(shù)據(jù)，2012年261個(gè)監(jiān)測(cè)數(shù)據(jù)，2013年70個(gè)監(jiān)測(cè)數(shù)據(jù)，結(jié)合邊坡實(shí)際情況及上下游區(qū)域概況，依經(jīng)驗(yàn)[2,16-17]假定該邊坡失事概率為1%，依前文建模方法，得到位移預(yù)警指標(biāo)估計(jì)值，并對(duì)預(yù)警指標(biāo)估計(jì)值序列進(jìn)行突變?cè)\斷。考慮到監(jiān)測(cè)數(shù)據(jù)具有非等時(shí)距性，為平緩數(shù)據(jù)，本文以該測(cè)點(diǎn)垂直位移日變化量(向下取正，單位：mm/d)為研究對(duì)象。

3.1監(jiān)測(cè)數(shù)據(jù)厚尾檢驗(yàn)

2009—2013年監(jiān)測(cè)數(shù)據(jù)厚尾檢驗(yàn)見(jiàn)圖1。由圖1可知，監(jiān)測(cè)數(shù)據(jù)中間部分近似直線(xiàn)，呈微凹型，下端翹起，具有厚尾分布，滿(mǎn)足極值POT模型成立的前提條件。

圖1　2009—2013年監(jiān)測(cè)數(shù)據(jù)厚尾檢驗(yàn)Fig.1　Heavy-tailed tests for monitoring data of 2009—2013

3.2Hill圖法擬定閾值u

根據(jù)1.4節(jié)中閾值的擬定方法，繪出2009—2013年Hill曲線(xiàn)圖(見(jiàn)圖2)，確定相應(yīng)年份閾值。

根據(jù)Hill圖中數(shù)據(jù)變化急緩情況，并結(jié)合已有學(xué)者文獻(xiàn)[10]的研究，擬定每年的k值，并由k值得到對(duì)應(yīng)的閾值u，假定邊坡失事概率為1%，由式(6)計(jì)算F(x)反函數(shù)，得到預(yù)警指標(biāo)估計(jì)值序列。結(jié)合年極值，將結(jié)果列于表1。為便于分析與比較，據(jù)表1繪制圖3。

圖2　2009—2013年Hill圖Fig.2　Hill figures of 2009—2013

圖3　預(yù)警指標(biāo)估計(jì)序列與極值序列對(duì)比Fig.3　Comparison between sequences of early 　warning indicator and extreme value

參數(shù)2009年2010年2011年2012年2013年k612163314u(mm·d-1)32.7631.5014.906.702.50位移預(yù)警指標(biāo)估計(jì)值(mm·d-1)41.3337.1437.2918.506.75年極值(mm·d-1)41.0737.4022.4013.004.60

3.3尖點(diǎn)突變模型診斷預(yù)警指標(biāo)估計(jì)值序列

3.4結(jié)論

據(jù)表1、圖3及Δi的計(jì)算結(jié)果可以看出：

(1)2009—2013年期間閾值、預(yù)警指標(biāo)估計(jì)值以及極值都呈下降趨勢(shì)，反應(yīng)邊坡向安全性態(tài)發(fā)展；

(2)2009—2010年間，該邊坡的位移預(yù)警值略小于年極值，該邊坡處于危險(xiǎn)狀態(tài)，2010年以后，預(yù)警值均高于年極值，該邊坡處于安全狀態(tài)；

(3)2011年時(shí)，狀態(tài)判別式Δ≈0，說(shuō)明位移預(yù)警指標(biāo)估計(jì)值序列處于臨界狀態(tài)，任之發(fā)展，將會(huì)發(fā)生突變，即Δ<0，之后卻發(fā)現(xiàn)Δ>0，顯而易見(jiàn)，邊坡位移預(yù)警指標(biāo)估計(jì)值向好的方向發(fā)生了轉(zhuǎn)變。2009年發(fā)現(xiàn)，該邊坡位移變化呈增大趨勢(shì)，便通過(guò)控制水庫(kù)水位升降、邊坡處理等工程措施加以控制，突變理論印證了采取控制措施后取得了實(shí)質(zhì)性效果。

4結(jié)語(yǔ)

極值理論P(yáng)OT模型以超閾值序列的GPD分布特征為研究對(duì)象，全面地考慮了所有較大測(cè)值，更好地刻畫(huà)了數(shù)據(jù)樣本的超閾值分布特征，因此得到的預(yù)警指標(biāo)估計(jì)值能更客觀地反映工程實(shí)際。

本文將基于極值理論的POT模型應(yīng)用于邊坡位移預(yù)警指標(biāo)估計(jì)中，通過(guò)分析研究，該邊坡處于安全狀態(tài)，并通過(guò)突變理論印證了控制措施起到了實(shí)質(zhì)性作用，這與邊坡實(shí)際情況較為吻合，與該邊坡實(shí)際監(jiān)測(cè)資料成果分析一致。研究分析表明極值POT模型與突變?cè)\斷模型相結(jié)合應(yīng)用于邊坡位移預(yù)警指標(biāo)估計(jì)與診斷上具有一定的科學(xué)性，較為可靠。

本文中邊坡失事概率尚無(wú)具體規(guī)范可依，依經(jīng)驗(yàn)對(duì)其假定，具有主觀性；監(jiān)測(cè)數(shù)據(jù)具有時(shí)間間隔性，為平緩數(shù)據(jù)使其滿(mǎn)足GPD分布，僅做了平均處理；建立POT模型的過(guò)程中，為簡(jiǎn)化分析，僅考慮超閾值序列數(shù)值大小而忽略超閾值序列的時(shí)間效應(yīng)。為更加客觀地估計(jì)預(yù)警指標(biāo)，存在的這些不足，需進(jìn)一步進(jìn)行研究與完善。

參考文獻(xiàn)：

[1]張振華,馮夏庭,周輝,等. 基于設(shè)計(jì)安全系數(shù)及破壞模式的邊坡開(kāi)挖過(guò)程動(dòng)態(tài)變形監(jiān)測(cè)預(yù)警方法研究[J]. 巖土力學(xué), 2009, 30(3): 603- 612. (ZHANG Zhen-hua, FENG Xia-ting, ZHOU Hui, et al. Research on dynamic early warning method of slope deformation monitoring during excavation based on design safety factor and failure mode[J]. Rock and Soil Mechanics, 2009, 30(3): 603- 612. (in Chinese))

[2]蘇懷智,王峰,劉紅萍. 基于POT模型建立大壩服役性態(tài)預(yù)警指標(biāo)[J]. 水利學(xué)報(bào), 2012, 43(8): 974- 978, 986. (SU Huai-zhi, WANG Feng, LIU Hong-ping. Early-warning index for dam service behavior based on POT model[J]. Journal of Hydraulic Engineering, 2012, 43(8): 974- 978, 986. (in Chinese))

[3]BALKEMA A A, DE HAAN L. Residual life time at great age[J]. Annals of Probability, 1974(2): 792- 804.

[4]PICKANDS J. Statistical inference using extreme order statistics[J]. Annals of Statistics, 1975(3): 119- 131.

[5]花擁軍,張宗益. 極值BMM與POT模型對(duì)滬深股市極端風(fēng)險(xiǎn)的比較研究[J]. 管理工程學(xué)報(bào), 2009, 23(4): 104- 108, 115. (HUA Yong-jun, ZHANG Zong-yi. Comparative research on extreme risk of stock market based on BMM and POT model[J]. Journal of Industrial Engineering and Engineering Management, 2009, 23(4): 104- 108, 115. (in Chinese))

[6]陳紅英. POT模型在巨災(zāi)保險(xiǎn)中的應(yīng)用[D]. 上海: 上海交通大學(xué), 2010.(CHEN Hong-ying. Application of POT model to catastrophe insurance[D]. Shanghai: Shanghai Jiaotong University, 2010. (in Chinese))

[7]錢(qián)小仕,王福昌,盛書(shū)中. 基于廣義帕累托分布的地震震級(jí)分布尾部特征分析[J]. 地震學(xué)報(bào), 2013, 35(3): 341- 350. (QIAN Xiao-shi, WANG Fu-chang, SHENG Shu-zhong. Characterization of tail distribution of earthquake magnitudes via generalized Pareto distribution[J]. Acta Seismologica Sinica, 2013, 35(3): 341- 350. (in Chinese))

[8]邊寬江,程波,王蕾蕾. 收益分布尖峰厚尾問(wèn)題的統(tǒng)計(jì)檢驗(yàn)[J]. 統(tǒng)計(jì)與決策, 2009(7): 83- 85. (BIAN Kuan-jiang, CHENG Bo, WANG Lei-lei. Problems about peak and fat-tail of income distribution based on statistics test[J]. Statistics & Decision, 2009(7): 83- 85. (in Chinese))

[9]HILL B M. A simple general approach to inference about the tail of a distribution[J]. Annals of Mathematical Statistics, 1975, 3(5): 1163- 1174.

[10]LORETAN M, PHILLIPS P C B. Testing the covariance stationarity of heavy-tailed time series: an overview of the theory with applications to several financial datasets[J]. Journal of Empirical Finance, 1994, 1(2): 211- 248.

[11]凌復(fù)華. 突變理論及其應(yīng)用[M]. 上海: 上海交通大學(xué)出版社, 1987. (LING Fu-hua. Catastrophe theory and its application[M]. Shanghai: Shanghai Jiaotong University Press, 1987.(in Chinese))

[12]阿諾德. 突變理論[M]. 周燕華, 譯. 北京: 高等教育出版社, 1990. (ARNOLD V I. Catastrophe theory[M]. Translated by ZHOU Yan-hua. Beijing: Higher Education Press, 1990. (in Chinese))

[13]姜璐,于連宇. 初等突變理論在社會(huì)科學(xué)中的應(yīng)用[J]. 系統(tǒng)工程理論與實(shí)踐, 2002(10): 113- 117. (JIANG Lu, YU Lian-yu. The application of primary catastrophe theory to social science[J]. Systems Engineering-Theory & Practice, 2002(10): 113- 117. (in Chinese))

[14]丁曉唐,顧沖時(shí),蔣勇. 基于有限元和突變理論的高拱壩壩體穩(wěn)定研究[J]. 河海大學(xué)學(xué)報(bào)：自然科學(xué)版, 2008, 36(2): 175- 178. (DING Xiao-tang, GU Chong-shi, JIANG Yong. Analysis of stability of high arch dams based on FEM and catastrophe theory[J]. Journal of Hohai University(Natural Sciences), 2008, 36(2): 175- 178. (in Chinese))

[15]鄭東健,雷霆. 基于突變理論的高拱壩失穩(wěn)判據(jù)研究[J]. 巖土工程學(xué)報(bào), 2011, 33(1): 23- 27. (ZHENG Dong-jian, LEI Ting. Instability criteria for high arch dams using catastrophe theory[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(1): 23- 27. (in Chinese))

[16]雷鵬,常曉林,肖峰,等. 高混凝土壩空間變形預(yù)警指標(biāo)研究[J]. 中國(guó)科學(xué):技術(shù)科學(xué), 2011, 41(7): 992- 999. (LEI Peng, CHANG Xiao-lin, XIAO Feng, et al. Study on early warning index of spatial deformation for high concrete dam[J]. Sci China Tech Sci,2011, 41(7): 992- 999. (in Chinese))

[17]華祖林,汪靚,顧莉,等. 基于門(mén)限極值理論的湖泊水質(zhì)參照狀態(tài)的確定[J]. 中國(guó)環(huán)境科學(xué), 2014, 34(12): 3215- 3222. (HUA Zu-lin, WANG Liang, GU Li, et al. Estimation of the lake quality reference condition based on the threshold extreme theory[J]. China Environmental Science, 2014, 34(12): 3215- 3222. (in Chinese))

Real-time estimation and diagnosis of early warning indicator of slope displacement based on POT-catastrophe theory

REN Jie1,2， SU Huai-zhi1,2， YANG Meng1,2， ZHOU Zhi-jie2

(1.StateKeyLaboratoryofHydrology-WaterResourcesandHydraulicEngineering,HohaiUniversity,Nanjing210098,China; 2.CollegeofWaterConservancyandHydropowerEngineering,HohaiUniversity,Nanjing210098,China)

Abstract：Real-time estimation and diagnosis of slope early warning indicators are important means to monitor the safety state of slope and prevent the situation of the evolution from slope to landslide. Safety early warning pays more attention to extreme events, which are not common. The key to achieve safety early warning lies in whether we can accurately analyze the quantile of the slope effect quantity sequences and describe the tail characteristics of the slope effect quantity sequences. First of all, the slope displacement monitoring sequence can be divided by the interval of time; then, the thresholds of the slope displacement monitoring sequences at each time interval are drafted using Hill figure method; the part of the displacement monitoring sequence at each time interval that is supra-threshold is fitted by Generalized Pareto Distribution (GPD), and the tail characteristics of the part that is supra-threshold of displacement monitoring sequences are described progressively. So, the ability of the slope to resist historical loads is analyzed based on the model of peaks over threshold (POT) using extreme value theory. If the failure probability of the slope is provided, the real-time estimation sequences of early warming indicators of displacement can be obtained and the ability of the slope to resist extreme loads that may happen in the future can be exploited and evaluated. Finally, the real-time estimation sequences of early warming indicators of displacement are diagnosed by using cusp catastrophic model from catastrophe theory. Based on specific slope projects, the rationality of the mixed model of POT and catastrophe theory is verified and the modeling method of the mixed model will be a new way to achieve the real-time estimation and diagnosis of engineering early warning indicators.

Key words：slope; early warning indicator of displacement; extreme value theory; POT model; cusp catastrophic theory

中圖分類(lèi)號(hào)：TV698.1

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1009-640X(2016)01-0030-07

作者簡(jiǎn)介：任杰(1992—)，男，山西臨汾人，碩士研究生，主要從事水工結(jié)構(gòu)工程安全監(jiān)控研究。

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(51179066，51409167，81279052，2014513311)；國(guó)家重點(diǎn)實(shí)驗(yàn)室研究資助項(xiàng)目(20145028312)； 江蘇省杰出青年基金資助項(xiàng)目(BK2012036)； 江蘇省“333高層次人才培養(yǎng)工程”科研項(xiàng)目(2016-B1307101)

收稿日期：2015-04-17

DOI:10.16198/j.cnki.1009-640X. 2016.01.005

任杰， 蘇懷智， 楊孟， 等. 邊坡位移預(yù)警指標(biāo)的實(shí)時(shí)估計(jì)與診斷[J]. 水利水運(yùn)工程學(xué)報(bào), 2016(1): 30-36. (REN Jie, SU Huai-zhi, YANG Meng, et al. Real-time estimation and diagnosis of early warning indicator of slope displacement based on POT-catastrophe theory[J]. Hydro-Science and Engineering, 2016(1): 30-36.)

E-mail: 1538406080@qq.com