李　匡，胡宇豐，梁犁麗，徐海卿

PSO算法在馬斯京根法參數(shù)率定中的應(yīng)用

李匡，胡宇豐，梁犁麗，徐海卿

（北京中水科水電科技開發(fā)有限公司，北京 100038）

將PSO算法用于馬斯京根法參數(shù)的自動率定，開發(fā)了基于PSO算法的馬斯京根法參數(shù)自動率定系統(tǒng)，該系統(tǒng)可以對馬斯京根法參數(shù)進(jìn)行快速的率定和檢驗；是分析率定結(jié)果穩(wěn)定性的一種新方法；此系統(tǒng)具有較高的實用性，已在多個流域中得到了應(yīng)用。

馬斯京根法；參數(shù)率定；PSO算法；參數(shù)自動率定系統(tǒng)；數(shù)據(jù)穩(wěn)定性

1　引言

馬斯京根法是一種廣泛應(yīng)用于河道洪水流量演算的方法，是美國工程師G.T.Mccarthy于1938年提出的流量演算法，此法最早在美國馬斯京根河流域上使用，因而稱為馬斯京根法。該法建立馬斯京根槽蓄曲線，并與水量平衡方程聯(lián)解，進(jìn)行河段洪水演算[1]。該方法在實際應(yīng)用中的一個重要問題是模型的參數(shù)估計。在傳統(tǒng)方法中利用試算法對參數(shù)K、x進(jìn)行率定，再計算出決策變量C0、C1、C2，但是試算法具有盲目性、不確定性、不易程序化且精度不高等缺陷[1-2]。因此有學(xué)者選擇了直接優(yōu)化決策變量C0、C1、C2，得到較好的效果，主要方法有最小二乘法[3]，遺傳算法[4-5]，粒子群算法[6]等。

粒子群優(yōu)化算法（ParticleSwarmOptimization Algorithm，PSO）是進(jìn)化算法的一種，系統(tǒng)初始化為一組隨機(jī)解，通過迭代搜尋最優(yōu)值。傳統(tǒng)PSO算法在應(yīng)用中和其它全局優(yōu)化算法一樣可能會陷入局部最優(yōu)，導(dǎo)致收斂精度不高，后期收斂速度較慢，經(jīng)過學(xué)者的研究對粒子群算法進(jìn)行了改進(jìn)，使得粒子群算法的效率得到了較大的提高[7-8]。另外，同其它的進(jìn)化算法，如遺傳算法相比，PSO模型具有參數(shù)較少，原理簡單，易于實現(xiàn)，計算效率高等優(yōu)點(diǎn)。鑒于以上優(yōu)點(diǎn)，本文選擇PSO算法率定馬斯京根法參數(shù)。

馬斯京根法參數(shù)率定是一個非線性優(yōu)化問題，利用PSO算法進(jìn)行率定是合適的。在應(yīng)用PSO算法進(jìn)行馬斯京根法參數(shù)率定時，需要給出馬斯京根法各參數(shù)的取值范圍，然而在已有的研究文獻(xiàn)中，作者并沒有明確給出各參數(shù)的取值范圍。本文根據(jù)嚴(yán)格的數(shù)學(xué)證明給出了馬斯京根法參數(shù)的取值范圍，為廣大學(xué)者在應(yīng)用各類優(yōu)化算法進(jìn)行馬斯京根法參數(shù)率定提供參考。

作者開發(fā)了基于PSO算法的馬斯京根法參數(shù)自動率定系統(tǒng)，按照給定的數(shù)據(jù)庫表結(jié)構(gòu)輸入資料（流量，斷面信息），系統(tǒng)可以對流域中各河道的馬斯京根法參數(shù)進(jìn)行自動率定，并對率定出的馬斯京根法參數(shù)進(jìn)行檢驗。

2　馬斯京根法簡介

馬斯京根法基本原理如下：對水量平衡方程（1）和馬斯京根槽蓄方程（2）進(jìn)行求解：

式中：I上斷面的流量；I1，I2時段始末上斷面的流量；Q下斷面的流量；Q1，Q2時段始末下斷面的流量；△t計算時段長；S河段蓄水量；S1，S2時段始末河段蓄水量；K穩(wěn)定流情況下的河段洪水傳播時間；x槽蓄系數(shù)。C0，C1，C2是 K，△t，x的函數(shù)，并且C0+C1+C2=1

3　粒子群（PSO）算法及其計算流程

PSO算法是Kennedy和Eberhart于1995年提出的一種基于對鳥群捕食行為模擬的智能群集優(yōu)化算法。其基本思想可以簡單地表述為：在一個優(yōu)化問題的解空間中，每一個可行解被看做一個“粒子”，這些粒子在解空間內(nèi)不停地飛行，在飛行的過程中根據(jù)自身和種群中其他粒子積累的經(jīng)驗不斷調(diào)整自己的飛行策略，最終這些粒子都趨于解空間中的最優(yōu)區(qū)域，也即所謂的“食物”[10]。

PSO算法的數(shù)學(xué)描述如下：在一個n維的搜索空間中，由m個粒子組成一個種群，即X={x1，x2，x3，…，xi，…xm}T，第i個粒子的位置為xi=（xi1，xi2，xi3，…，xin），其速度為Vi=（Vi1，Vi2，Vi3，…，Vin），粒子的個體極值為：Pi=（Pi1，Pi2，Pi3，…，Pin），種群的全局極值為：Pg=（Pg1，Pg2，Pg3，…，Pgn），基本PSO算法的粒子在搜索過程中通過式（3）和式（4）不斷地進(jìn)行位置和飛行速度的更新：

其中：d=1，2，3，…，n；i=1，2，3，…，m。

式中：n為粒子的維數(shù)；m為種群的規(guī)模；t為當(dāng)前的迭代步數(shù)；r1，r2分別為0和1之間的隨機(jī)數(shù)，c1，c2為加速度常數(shù)。

粒子的維數(shù)n與具體的研究問題有關(guān)，例如在率定馬斯京根法參數(shù)時，有三個參數(shù)需要率定，則n=3。

種群規(guī)模m影響算法的穩(wěn)定性和效率，m取值越大，PSO算法的穩(wěn)定性越好，但同時計算量也較大，耗時長，計算效率較低。因此在對精度要求較高時可選擇較大的種群規(guī)模，對計算效率要求較高時選擇較小的種群規(guī)模，種群規(guī)模m的取值范圍一般在[50，100]之間。

從式（3）可以看出，粒子飛行速度更新公式包括三部分：第一部分是粒子的歷史飛行速度，說明了粒子目前的狀態(tài)，起到平衡全局和局部搜索能力的作用；第二部分是粒子自身的認(rèn)識，表示粒子在飛行過程中自身的思考；第三部分是粒子的社會認(rèn)識，表示粒子群中各粒子之間信息的交流。三個部分共同作用，決定了粒子的空間搜索能力。

傳統(tǒng)的粒子群算法容易陷入局部極小值，影響PSO算法找到全局最優(yōu)值。為了提高PSO算法的優(yōu)化性能，平衡算法的全局搜索能力和局部搜索能力，文獻(xiàn)[11]對基本PSO算法進(jìn)行了改進(jìn)，在基本PSO算法的基礎(chǔ)上引進(jìn)了慣性權(quán)重項ω，提出了標(biāo)準(zhǔn)PSO算法。式（3）也相應(yīng)地修正為：

式中ω的取值方式有兩種。一種是遞減策略，慣性權(quán)重ω滿足ω（t）=0.9-t×0.5，其中t為當(dāng)前迭T 代次數(shù)，T為最大迭代次數(shù)；另外一種是固定策略，通常取慣性權(quán)重為[0.4，0.6]中的某一固定值。

根據(jù)PSO算法基本思想，可以設(shè)計標(biāo)準(zhǔn)PSO算法的流程如下：

（1）將優(yōu)化的問題數(shù)學(xué)模型化，選定優(yōu)化問題的目標(biāo)函數(shù)（即適應(yīng)度函數(shù)），選定粒子的維數(shù)；

（2）初始化算法。對粒子群中的粒子位置和速度進(jìn)行初始化設(shè)定，即在一定的范圍內(nèi)隨機(jī)產(chǎn)生出每一個粒子的位置和速度；

（3）根據(jù)優(yōu)化問題的目標(biāo)函數(shù)計算每個粒子的適應(yīng)度值；

（4）對每個粒子，將其當(dāng)前適應(yīng)度值與其所經(jīng)歷的最優(yōu)適應(yīng)度值進(jìn)行比較，如果該粒子的當(dāng)前適應(yīng)度值更優(yōu)，那么將當(dāng)前位置記錄為該粒子的局部最優(yōu)位置；

（5）對每個粒子所經(jīng)歷的最優(yōu)適應(yīng)度值與全局最優(yōu)適應(yīng)度值進(jìn)行比較，如果個體粒子的最優(yōu)適應(yīng)度值較全局最優(yōu)適應(yīng)度值為優(yōu)，則將其作為全局最優(yōu)位置；

（6）對每個粒子的位置和飛行速度進(jìn)行更新；

（7）判斷是否達(dá)到優(yōu)化的終止條件。如果滿足終止條件，結(jié)束循環(huán)，否則返回第三步。

4　馬斯京根法參數(shù)率定系統(tǒng)

本文利用Vb.net開發(fā)了馬斯京根法參數(shù)率定系統(tǒng)，該系統(tǒng)利用PSO算法對馬斯京根法參數(shù)進(jìn)行率定。將各水文站的測站信息、流量資料等按照系統(tǒng)設(shè)計的數(shù)據(jù)庫表結(jié)構(gòu)輸入到數(shù)據(jù)庫中，系統(tǒng)即可按照默認(rèn)設(shè)置的PSO算法參數(shù)對馬斯京根法參數(shù)進(jìn)行率定，并且可以對率定的參數(shù)進(jìn)行檢驗及分析。

在用PSO算法進(jìn)行馬斯京根法參數(shù)率定時，還需要解決以下幾個問題：

4.1目標(biāo)函數(shù)

采用PSO算法率定模型參數(shù)時必須構(gòu)建目標(biāo)函數(shù)，一方面用以評價模擬效果的優(yōu)劣，另一方面作為粒子飛行策略調(diào)整的標(biāo)準(zhǔn)。在對馬斯京根法參數(shù)進(jìn)行率定時，一般有河槽蓄量誤差最小和流量誤差最小兩種目標(biāo)函數(shù)，經(jīng)學(xué)者研究表明，選用流量誤差最小作為目標(biāo)函數(shù)率定的馬斯京根法參數(shù)效果更好[12]。

本系統(tǒng)將計算流量和實測流量的誤差平方和最小作為判斷馬斯京根法參數(shù)尋優(yōu)的目標(biāo)函數(shù)：

式中：Qi為第i時段計算出流量，QiC為第i時段實測出流量，n為總時段數(shù)。

4.2參數(shù)率定結(jié)果的穩(wěn)定性分析

在利用PSO算法率定馬斯京根法參數(shù)時，由于PSO算法的初始解是隨機(jī)給定的，因此可能會造成率定的參數(shù)不穩(wěn)定。該問題可以通過多次計算，得到率定參數(shù)結(jié)果序列，然后檢驗參數(shù)結(jié)果序列是否穩(wěn)定，如果穩(wěn)定，則取參數(shù)結(jié)果的平均值作為最終率定參數(shù)結(jié)果；如果未通過穩(wěn)定性檢驗，則調(diào)整PSO參數(shù)重新率定，直到率定的參數(shù)結(jié)果穩(wěn)定為止。

可以通過以下算法檢驗率定的參數(shù)結(jié)果是否穩(wěn)定：

該算法的核心是計算聚集在平均值附近的數(shù)值的概率，當(dāng)該概率大于穩(wěn)定指標(biāo)時，則確定該數(shù)據(jù)序列是穩(wěn)定的，反之?dāng)?shù)據(jù)不穩(wěn)定。判斷過程可按照以下步驟進(jìn)行：

（1）分析數(shù)據(jù)的特點(diǎn)，確定合適的穩(wěn)定性指標(biāo)S以及穩(wěn)定性判別變化幅度Adec；

（2）計算該段數(shù)據(jù)序列的平均值Vavg；

（3）給定變化幅度A（0≤A≤1）圍繞平均值Vavg上、下波動，得到區(qū)間[Vavg-Vavg×A，Vavg+Vavg×A]，稱此區(qū)間為穩(wěn)定區(qū)域F。計算穩(wěn)定區(qū)域內(nèi)數(shù)據(jù)出現(xiàn)的概率P，即穩(wěn)定區(qū)域內(nèi)出現(xiàn)的數(shù)據(jù)個數(shù)與整個時間段內(nèi)數(shù)據(jù)的總個數(shù)的比值。如果P大于穩(wěn)定指標(biāo)F并且abs（P-S）≥ε，則轉(zhuǎn)至（4）繼續(xù)；如果P小于穩(wěn)定指標(biāo)S并且abs（P-S）≥ε，則轉(zhuǎn)至（5）繼續(xù)；如果abs（P-S）≥ε，則轉(zhuǎn)至（6）；

（4）減小變化幅度A，轉(zhuǎn)至（3）重新計算；

（5）增大變化幅度A，轉(zhuǎn)至（3）重新計算；

（6）如果A≤Adec，則認(rèn)為數(shù)據(jù)序列穩(wěn)定，否則認(rèn)為數(shù)據(jù)序列不穩(wěn)定。

此法的優(yōu)點(diǎn)是：不僅可以判斷數(shù)據(jù)是否穩(wěn)定，而且可以判別出數(shù)據(jù)在多大的范圍內(nèi)穩(wěn)定。

4.3有支流斷面的參數(shù)率定問題

對有支流匯入的斷面，系統(tǒng)中采用“先合后演”的方法進(jìn)行參數(shù)率定，即將干、支流所有上斷面相同時刻的流量相加，作為河段總?cè)肓?，按無支流河段方法進(jìn)行演算[2]。

4.4系統(tǒng)中PSO算法默認(rèn)參數(shù)

（1）計算次數(shù)：ComputeCount=50；

（2）種群規(guī)模：m=70；

（3）粒子維數(shù)：n=70；

（5）限制條件C2=1-C0-C1，當(dāng)C2超出其取值范圍，則調(diào)整至其取值范圍；

（6）加速因子：設(shè)定為C1=C2=2.0；

（7）慣性權(quán)重：從0.9～0.4遞減變化或者取固定值ω=0.5，默認(rèn)選擇后一種；

（8）速度限制系數(shù)：k=0.7；

（9）終止條件：最大迭代次數(shù)T=200，或者目標(biāo)函數(shù)之差小于給定精度ε=10-8。

表1　率定誤差與檢驗誤差比較

在實際應(yīng)用中用戶可以通過對以上參數(shù)進(jìn)行修改，找出率定效果較好的馬斯京根法參數(shù)。

5　實例

采用雅礱江流域上三個測站的資料進(jìn)行馬斯京根法參數(shù)率定和檢驗。三個測站的地理位置示意圖如圖1所示：

圖1　測站地理位置示意圖

將麥地龍站和列瓦站同時刻的資料相加作為入流量，將錦屏站的資料作為出流量。利用2006～2008年的資料進(jìn)行參數(shù)率定，再利用2009年的資料進(jìn)行檢驗，我們得到如下結(jié)果：

率定結(jié)果C0=0.056，C1=0.001，C2=0.942。

資料的率定誤差和檢驗誤差分別如表1。

先對每次計算所得的最優(yōu)目標(biāo)函數(shù)值序列以及參數(shù)序列分別進(jìn)行歸一化處理；再計算目標(biāo)函數(shù)值的方差為0.1466，率定參數(shù)的方差分別為0.0396，0.0218，0.0849，證明PSO算法具有較高的穩(wěn)定性。

每次計算當(dāng)目標(biāo)函數(shù)達(dá)到終止條件所需的平均迭代次數(shù)為27.66次，最大次數(shù)為54次，最小次數(shù)為13次，因此PSO算法也有較高的效率。

由于在參數(shù)率定時，選用的資料序列較長，計算量較大（因為計算50次需要的時間較長）。鑒于PSO算法的穩(wěn)定性和效率，在實際應(yīng)用中可以不用計算50次，只需計算5次左右即可，以期節(jié)省計算時間。另外，如果對精度要求不高，可以將終止條件中的目標(biāo)函數(shù)之差設(shè)為ε=10-5，此舉也可節(jié)省計算時間。

6　結(jié)束語

本文將PSO算法用于馬斯京根法參數(shù)的自動率定中，開發(fā)了基于PSO算法的馬斯京根法參數(shù)自動率定系統(tǒng)，此系統(tǒng)具有較高的實用價值，適合內(nèi)嵌入洪水預(yù)報系統(tǒng)進(jìn)行參數(shù)率定和檢驗；該系統(tǒng)在雅礱江流域和豐滿流域等的洪水預(yù)報過程中均得到了應(yīng)用，取得了較好的參數(shù)率定和檢驗效果。但是當(dāng)率定參數(shù)采用的資料較長時，需要的時間也較長，下一步我們將再對PSO算法進(jìn)行改進(jìn)，以提高系統(tǒng)的運(yùn)算效率。

[1]詹道江，葉守澤.工程水文學(xué)[M].北京:中國水利水電出版社，2000.

[2]包為民.水文預(yù)報[M].3版.北京:中國水利水電出版社，2006.

[3]翟國靜.馬斯京根模型參數(shù)估計方法探討[J].水文，1997，3: 40-43.

[4]魯帆，蔣云鐘，王浩，等.多智能體遺傳算法用于馬斯京根模型參數(shù)估計[J].水利學(xué)報，2007，38（3）：289-294.

[5]董純，康玲，侯國祥.基于遺傳模擬退火法的馬斯京根方程參數(shù)估計[J].華中科技大學(xué)學(xué)報（自然科學(xué)版），2002，30

（10）:20-22.

[6]李明明，李承軍，張銘.改進(jìn)PSO法在馬斯京根模型參數(shù)估計中的應(yīng)用[J].人民長江，2008，39（3）：60-62.

[7]李建勇.粒子群優(yōu)化算法研究[D].杭州:浙江大學(xué)，2004.

[8]熊鷹.微粒群算法的若干改進(jìn)及應(yīng)用[D].武漢:武漢理工大學(xué)，2006.

[9]李匡，付力，胡宇豐，等.馬斯京根法參數(shù)C0，C1，C2取值范圍的確定[J].南水北調(diào)與水利科技，2012，10（5）:43-45.

[10]劉蘇寧，甘泓，魏國孝.粒子群算法在新安江模型參數(shù)率定中的應(yīng)用[J].水利學(xué)報，2010，41（5）:537-544.

[11]ShiY，EberhartRC．EmpiricalStudyofParticleSwarm Optimization[C]//Proceedingofthe1999CongressonEvolutionaryComputation．Piscataway，NJIEEEServiceCenter，1999．

[12]董曉華，劉冀，徐剛，等.目標(biāo)函數(shù)對馬斯京根法河道流量演算精度的影響[J].水文，2009，29（6）：9-13.

TV734

B

1672-5387（2016）08-0088-04

10.13599/j.cnki.11-5130.2016.08.027

2016-06-29

李匡（1983-），男，高級工程師，從事洪水預(yù)報調(diào)度系統(tǒng)開發(fā)工作。