劉云潺,張慧寧
(黃河水利職業(yè)技術(shù)學(xué)院,河南 開封 475004)
在節(jié)水灌溉中,需要對農(nóng)田土壤水分情況進行測量[1]。人工監(jiān)測,既耗費人力,又不能實時監(jiān)控,因此通過土壤水分特征曲線描述土壤含水量與吸力之間的關(guān)系,對于研究土壤水分特征十分重要[2]。
土壤水分測量與土壤本身的機理無關(guān),最著名的是Van Genuchten模型,模型中參數(shù)的物理意義明確且與實測曲線擬合較好[3],但是Van Genuchten模型屬于非線性復(fù)雜方程并且參數(shù)較多,擬合參數(shù)過程屬于非線性問題,若使用線性方法和線性擬合會出現(xiàn)較大誤差。目前解決方法有:最小二乘法(Least Square Method,LSM)擬合Van Genuchten模型方程參數(shù)[4],能夠?qū)?shù)的求解提供較好的方法,但是會遇到求解停止或參數(shù)為負(fù)的問題;粒子群算法(Particle Swarm Optimization,PSO),編程簡單[5],但是在Van Genuchten模型參數(shù)尋優(yōu)后期容易產(chǎn)生停滯現(xiàn)象;差異進化粒子群算法(Differential Evolution PSO,DEPSO),考慮了多參數(shù)間的相關(guān)性,能夠在一定程度上簡化方程參數(shù)估計[6];模擬退火算法(Simulated Annealing,SA)、遺傳算法(Genetic Algorithm,GA),均能夠獲得參數(shù)解[7,8],但全局收斂性相對較弱。
本文引入平行、競爭策略,提出平行競爭PSO算法(Parallel Competition PSO,PCPSO),競爭機制通過更新部分粒子來改善群體的多樣性同時保持最好的粒子,最差粒子被重置;平行策略通過把粒子群分成若干子群,引入共享因子對子群之間信息共享進行動態(tài)調(diào)節(jié),實驗仿真顯示本文算法對測試函數(shù)求解具有收斂速度較快、解精度較高的特點,測量粉壤土水分的相對誤差相比其他算法都較小。
粒子群算法(Particle Swarm Optimization,PSO)通過學(xué)習(xí)因子c1和c2及慣性權(quán)重ω進行更新速度vq,h和位置xq,h[9]:
(1)
式中:q=1,2,…,h=1,2,…,t為迭代次數(shù);xq,h(t)為第q個粒子在第t次迭代時第 維的位置;vq,h(t)為第q個粒子在第t次迭代時第 維的速度;pg,h為歷史最優(yōu)值;pq,h為當(dāng)前最優(yōu)值;r1∈(0,1)、r2∈(0,1)為相互獨立的隨機函數(shù)。
經(jīng)典粒子群優(yōu)化算法存在早熟收斂和停滯現(xiàn)象,因為當(dāng)評價函數(shù)高維或者形狀復(fù)雜,粒子最終被困在一個局部極小值范圍內(nèi)。
1.2.1 競爭機制
通過更新部分粒子來改善群體的多樣性同時保持最好的粒子,最差粒子被重置,以便更好地探索Van Genuchten模型參數(shù)[10]。新生成的粒子將試圖尋找更好的最優(yōu)值,這樣最好的粒子圍繞當(dāng)前的全局最佳位置保持搜索。如果最新的粒子找到更好的最優(yōu)值,粒子群就會聚集在全局最佳周圍,否則它們將回到當(dāng)前全局最佳狀態(tài),直到再次被重置。只有粒子在檢測到早熟收斂時才被觸發(fā)重置,這樣粒子群獲得了內(nèi)部競爭。粒子狀態(tài)檢測到需要重置的方法如下:在第t次迭代時,如果粒子群最大半徑值δt小于閾值ε時,競爭才被觸發(fā)。ε表達式根據(jù)經(jīng)驗確定,不但對于控制競爭觸發(fā)效果要求較高,而且要阻止群體中任何粒子在局部最小值中停滯太久。最終δt與ε設(shè)計為:
(2)
式中:‖‖為歐幾里得范數(shù);s為粒子總數(shù)。
粒子被重置的比例應(yīng)該隨著迭代次數(shù)增加而減少,不應(yīng)該一直以一定的規(guī)模被重置,因此通過非線性策略,粒子在迭代t上由Logistic函數(shù)σ(t)重新設(shè)置比例,σ(t)減少與迭代次數(shù)的線性相關(guān)性。通過引入一個競爭性參數(shù)γ∈[0,2]控制:
(3)
Logistic函數(shù) 與不同的 值關(guān)系如圖1所示。
圖1 σ(t)與不同的γ值Fig.1 σ(t)and different γ values
通過圖1可以發(fā)現(xiàn),當(dāng)γ值較大時,有助于增加粒子之間的競爭,但是收斂較慢;當(dāng)γ值較小時,粒子之間的競爭減弱,但是收斂較快,γ=0時,不存在競爭性。本文令γ=1,這樣可以保持競爭性同時加快收斂速度。
1.2.2 平行策略
粒子之間的搜索信息交流采用平行子群方式[11,12],把粒子群分成若干子群,子群的群平均適應(yīng)度與原始粒子群平均適應(yīng)度相差不能小于設(shè)定的閾值ζ,每個子群由若干個粒子構(gòu)成。ζ計算為:
(4)
(5)
(6)
式中:fn為各個子群的適應(yīng)度值;f為歸一化因子,調(diào)節(jié)ζ大小。
在搜索過程中,為了控制粒子個體之間信息共享的程度,實現(xiàn)子群中較差個體的更新,引入共享因子對子群之間信息共享進行動態(tài)調(diào)節(jié),任一個子群中的粒子與周圍子群信息共享次數(shù)為e1∈[1,E],E為最大共享總次數(shù),共享因子ψ1計算為:
(7)
子群中的粒子與自身子群信息共享次數(shù)為e2∈[1,E′],E′為最大共享總次數(shù),共享因子ψ2計算為:
(8)
(9)
式中:s為sigmoid 型隸屬函數(shù)。
sigmoid 型隸屬函數(shù)控制共享過程如圖2所示。
圖2 sigmoid 型隸屬函數(shù)控制共享過程Fig.2 Sigmoid type membership function controls sharing process
從圖2可以看出,e1/e2較小時尋優(yōu)前期共享主要在自身子群中進行;當(dāng)e1/e2達到一定程度時,s迅速增大,共享尋優(yōu)主要與周圍子群進行,擴大搜索范圍,避免局部解。
則粒子更新速度vq,h和位置xq,h為:
(10)
土壤水分特征曲線Van Genuchten模型為:
(11)
式中:θ為土壤體積含水量,cm3/cm3;h為土壤基質(zhì)勢,cm;θr和θs分別為土壤的剩余體積含水量和飽和體積含水量,cm3/cm3;a和n均為經(jīng)驗擬合參數(shù),m=1-1/n。
參數(shù)θs、θr、a、n相互獨立,確定需要優(yōu)化變量的取值范圍θs∈[0,0.6]、θr∈[0,0.3]、a∈[0,1]、n∈[1,10]。但是由于方程的參數(shù)較多,常規(guī)的線性處理與線性擬合無法真實實現(xiàn)數(shù)據(jù)曲線,需要通過非線性擬合方法實現(xiàn)[13-15]。
把均方誤差作為Van Genuchten方程曲線擬合度評價標(biāo)準(zhǔn):
(12)
則粒子群適應(yīng)度函數(shù)f(h)選擇為:
(13)
當(dāng)f(h)獲得最小的時候則Van Genuchten方程曲線擬合度較高。
算法流程:①粒子群初始化,產(chǎn)生粒子個體;②計算粒子的適應(yīng)度;③按照公式(2)、(3)進行競爭操作;④按照公式(7)、(8)進行平行子群操作;⑤按照公式(10)更新粒子群;⑥若滿足最大迭代次數(shù)或Van Genuchten方程曲線擬合度,精度控制在10-6,滿足進行步驟⑦,否則進行步驟②;⑦輸出θs、θr、a、n參數(shù)值。
實驗PC配置為CPU3.6GHz、內(nèi)存8GB,集成顯卡,Matlab2014編程實現(xiàn)仿真,粒子群設(shè)置為400個,最大迭代次數(shù)為300,子群數(shù)量最大為20個,c1=c2=2,ω∈(0.2,0.8)。測量土壤水分傳感器采用時域反射型,探針長度:8.5 cm,探針直徑:3 mm,密封材料:環(huán)氧樹脂,探針材料:304不銹鋼。
為了檢測PCPSO算法的性能,選擇4個復(fù)雜函數(shù)Sphere function、Rosenbrock function、Rastrigin function、Quartic function進行測試,分別如下:
圖3至圖6給出了4個測試函數(shù)求解過程中平均最優(yōu)值的變化曲線。
圖3 Sphere function 變化曲線Fig.3 Sphere function variation curve
圖4 Rosenbrock function變化曲線Fig.4 Rosenbrock function variation curve
圖5 Rastrigin function 變化曲線Fig.5 Rastrigin function variation curve
從圖3至圖6可以看出PCPSO算法仿真實驗結(jié)果較好,PCPSO算法對測試函數(shù)都能搜索到比PSO更優(yōu)的結(jié)果,收斂速度較快,同時解的精度較高。這是因為PCPSO算法在迭代過程中,通過閾值 觸發(fā)競爭機制保持適應(yīng)度值較大的粒子優(yōu)勢,避免早熟,平行策略向全局最優(yōu)位置收斂,進而使得方程解優(yōu)化達到最優(yōu)。
圖6 Quartic function 變化曲線Fig.6 Quartic function variation curve
在求解粉壤土Van Genuchten模型參數(shù)過程中,使用了LSM、PSO、DEPSO、SA 、GA、PCPSO方法,實驗中每個算法采用蒙特卡羅方法法獨立運行30次,結(jié)果取平均值獲得參數(shù)解如表1所示。
對粉壤土進行不同算法模擬數(shù)據(jù)與實測數(shù)據(jù)對比如表2所示。
根據(jù)表1獲得的θs、θr、a、n參數(shù)值,從表2可以看出,本文算法對粉壤土模擬數(shù)據(jù)接近實測數(shù)據(jù)。不同方法對粉壤土脫濕、吸濕θ~h實測數(shù)據(jù)的相對誤差比較如圖7和圖8所示。
表1 不同方法確定粉壤土Van Genuchten模型參數(shù)Tab.1 Determination of Van Genuchten model parameters for silt loam by different methods
圖7 不同方法對粉壤土脫濕θ~h實測數(shù)據(jù)的相對誤差比較Fig.7 Comparison of relative error of θ~h measured by different methods for silt loam dehumidification soil
從圖7和圖8可以看出,本文PCPSO算法對粉壤土脫濕、吸濕θ~h實測數(shù)據(jù)相對誤差比較小,脫濕實測數(shù)據(jù)相對誤差最大為5%,吸濕實測數(shù)據(jù)相對誤差最大為4%,其他算法都比PCPSO大,從而在實際測試中會導(dǎo)致較大的誤差。這是因為PCPSO算法在可行解空間內(nèi)通過多個子群從多個初始解開始搜索,而競爭策略防止搜索過程收斂于局部最優(yōu)解,從而獲得全局最優(yōu)解。
表2 不同方法模擬粉壤土θ~h數(shù)據(jù)與實測數(shù)據(jù)對比 cmTab.2 Comparison of θ~h data and measured data for different methods of simulating loam soil
圖8 不同方法對粉壤土吸濕θ~h實測數(shù)據(jù)的相對誤差比較Fig.8 Comparison of relative error of θ~h measured by different methods for silt loam absorption soil
本文提出平行競爭策略PSO算法,競爭機制更新最差的粒子,平行策略引入共享因子對子群之間信息共享進行動態(tài)調(diào)節(jié),實驗仿真顯示本文算法對測試函數(shù)求解具有收斂速度較快、解精度較高的特點,模擬測量粉壤土水分的相對誤差相比其他算法都較小,因此為測量土壤水分提供了一種新方法。