饑餓搜索算法在地下水參數(shù)反演中的應(yīng)用

王文川，薛一銘，徐 雷

（1.華北水利水電大學(xué)水資源學(xué)院，鄭州 450046；2.河海大學(xué)水文水資源學(xué)院，南京 210024）

0 引 言

區(qū)域內(nèi)某一時(shí)段內(nèi)的地下水資源通常具有隨機(jī)性和非平穩(wěn)性特征，對(duì)區(qū)域水文地質(zhì)參數(shù)（導(dǎo)水系數(shù)T 和儲(chǔ)水系數(shù)S）的精確估計(jì)，是進(jìn)行模擬地下水運(yùn)動(dòng)數(shù)學(xué)模型、科學(xué)利用地下水資源的關(guān)鍵。因此，需要研究探討簡(jiǎn)便實(shí)用的參數(shù)反演方法來(lái)輔助地下水資源開(kāi)發(fā)利用［1］。一般而言，為獲得較為可靠的參數(shù)，首先需要進(jìn)行非穩(wěn)定流抽水試驗(yàn)，再基于泰斯公式［2］利用試驗(yàn)數(shù)據(jù)反演分析。目前，傳統(tǒng)反演方法有泰斯配線法［3］、改進(jìn)直線圖解法［4］和線性回歸法［5］、Jacob 直線圖解法［6］和簡(jiǎn)算法［7］等，但這些方法求解過(guò)程復(fù)雜，得到的參數(shù)可能不是全局最優(yōu)解。為了更精確地求解參數(shù)，近些年來(lái)，許多學(xué)者已開(kāi)始將智能優(yōu)化算法應(yīng)用于求解參數(shù)反演的問(wèn)題。崔東文［8］將飛蛾火焰優(yōu)化（Moth-Flame Optimization，MFO）算法應(yīng)用于含水層參數(shù)反演，然而MFO 算法不可避免地存在早熟收斂和全局搜索精度低的缺陷。邱淑偉等［9］利用遍歷搜索算法進(jìn)行求解水文地質(zhì)參數(shù)，避免了陷入局部最優(yōu)解的情況，求參過(guò)程受人主觀因素影響較?。坏枰啻嗡阉?，當(dāng)搜索空間范圍較大時(shí)效率不高。張娟娟［10］采用改進(jìn)模擬退火算法雖較好地處理了地下水參數(shù)反演問(wèn)題，但其所求得的參數(shù)精度不高。陳童［11］等采用自適應(yīng)的變異因子F和交叉因子CR改進(jìn)差分進(jìn)化算法，將其應(yīng)用于泰斯井函數(shù)，較傳統(tǒng)差分進(jìn)化算法尋優(yōu)率與計(jì)算結(jié)果精度更高。譚璟［12］設(shè)計(jì)了一種模擬退火—人工魚(yú)群的混合算法將其應(yīng)用到地下水模型參數(shù)反演問(wèn)題中。該方法相較原始人工魚(yú)群算法而言，極大地節(jié)省了參數(shù)計(jì)算時(shí)間。秋云翔等［13］提出改進(jìn)的入侵雜草優(yōu)化算法應(yīng)用于泰斯模型中，然而當(dāng)導(dǎo)水系數(shù)和儲(chǔ)水系數(shù)較小時(shí)，所得結(jié)果適應(yīng)度和誤差值與常用方法接近。

要想得到理想的反演結(jié)果，要求智能優(yōu)化算法具有較好的全局尋優(yōu)能力和收斂性。隨著算法研究的不斷發(fā)展，近幾年涌現(xiàn)出新一批如天鷹優(yōu)化算法［14］、鯨魚(yú)優(yōu)化算法［15］、黃金正弦算法［16］以及阿基米德算法［17］等新型算法，它們?cè)诨鶞?zhǔn)函數(shù)測(cè)試上已表現(xiàn)出較強(qiáng)優(yōu)越性。但“沒(méi)有免費(fèi)的午餐”定理在邏輯上證明了沒(méi)有一種算法可以作為最好的方法解決所有的優(yōu)化問(wèn)題?；诖?，本文引入Yutao Yang 等人［18］提出的一種新型的饑餓搜索算法，基于兩例抽水試驗(yàn)數(shù)據(jù)［10，19］，利用HGS 算法求解相關(guān)水文地質(zhì)參數(shù)，并與以上多種算法反演結(jié)果進(jìn)行比較。本項(xiàng)工作旨在驗(yàn)證HGS 算法在地下水參數(shù)反演問(wèn)題中的可行性和高效性，為其他相關(guān)地下水模型研究工作提供參考。

1 基本原理

1.1 算法原理

饑餓是動(dòng)物生活中決定行為和行動(dòng)最重要的穩(wěn)態(tài)動(dòng)機(jī)和原因之一。盡管自然界各種各樣的刺激和相互競(jìng)爭(zhēng)的需求總是影響動(dòng)物的生活質(zhì)量，但當(dāng)它們面臨熱量不足時(shí)，它們必須尋找食物。健康的動(dòng)物比弱小的更容易找到食物來(lái)源，有更大的生存機(jī)會(huì)。饑餓越強(qiáng)烈，對(duì)食物的渴望越強(qiáng)烈，機(jī)體就越積極地在短時(shí)間內(nèi)尋找食物，當(dāng)食物來(lái)源有限時(shí)，為避免因饑餓導(dǎo)致死亡，饑餓的動(dòng)物之間會(huì)有一個(gè)尋找食物來(lái)源并取得勝利的博弈?；诖耍琘ang 等人［18］于2021年提出了一種基于群居動(dòng)物因饑餓進(jìn)行食物搜索行為的群智能算法——饑餓搜索算法。其數(shù)學(xué)模型描述如下：

（1）接近食物。社會(huì)動(dòng)物在覓食過(guò)程中經(jīng)常相互合作，但不排除少數(shù)個(gè)體不參與合作的可能性。以下公式代表了HGS算法在個(gè)體合作交流和覓食行為中的核心方程：

對(duì)所有位置的變化控制的E的計(jì)算公式如下：

式中：i∈1，2，…，n，F(xiàn)（i）代表每個(gè)個(gè)體的適應(yīng)度值；BF是在當(dāng)前迭代過(guò)程中獲得的最佳適應(yīng)度值；sech 是雙曲函數(shù)，sech(x) =

式中：rand∈?［0，1］；T代表最大迭代數(shù)。

式（1）允許個(gè)體探索接近最優(yōu)解的可能位置和遠(yuǎn)離最優(yōu)解的位置，在一定程度上保證了對(duì)解空間邊界內(nèi)所有位置的搜索多樣化。HGS 在搜索過(guò)程中，執(zhí)行式中的Game1可以完成全局搜索，無(wú)論是Game2還是Game3都可以完成局部搜索。同樣的概念也可以應(yīng)用于高維搜索空間。

（2）饑餓角色。在這一部分，個(gè)體在搜索中的饑餓特征通過(guò)一個(gè)提出的式（1）進(jìn)行數(shù)學(xué)模擬：

式中：hungry 代表每個(gè)個(gè)體的饑餓值；N 代表個(gè)體的數(shù)量；SHungry 所有個(gè)體饑餓感的總和，即sum（hungry）；r3，r4，r5∈［0，1］；hungry（i）的公式如下：

式中：AllFitness（i）是當(dāng)前迭代過(guò)程中第i 個(gè)個(gè)體的適應(yīng)度。在每次迭代中，最佳個(gè)體的饑餓被設(shè)置為0，對(duì)于其他個(gè)體，新的饑餓（H）是在原有饑餓的基礎(chǔ)上增加的，不同個(gè)體對(duì)應(yīng)的H 是不同的。H的計(jì)算公式為：

式中：r6∈［0，1］；F（i）代表每個(gè)個(gè)體的適應(yīng)值；BF 是當(dāng)前迭代過(guò)程中群體獲得的最佳適應(yīng)度值；而WF 則代表相應(yīng)的最差適應(yīng)度值；UB和LB分別為搜索空間的上限和下限；為了使算法得到更好的性能，通過(guò)控制饑餓的上下限，LH 為饑餓感H 下限。公式（8）中的UB 和LB 用來(lái)反映個(gè)體在不同環(huán)境中最大的饑餓值；F（i）-BF 代表個(gè)體不再需要挨餓的食物量；WF-BF 得出個(gè)體在當(dāng)前過(guò)程中的總覓食能力；代表饑餓比率；r6×2代表消除環(huán)境因素對(duì)饑餓的正面或負(fù)面影響。

由于以下獨(dú)特優(yōu)勢(shì)，饑餓搜索算法表現(xiàn)出高效的性能：①參數(shù)l 和E 協(xié)助HGS 演化初始位置和改變搜索模式，以確保能盡可能地探索整個(gè)解空間，并在很大程度上增強(qiáng)算法的多樣化能力；②饑餓權(quán)重在搜索過(guò)程中增加對(duì)HGS 的擾動(dòng)，這兩個(gè)動(dòng)態(tài)自適應(yīng)參數(shù)可以防止算法陷入局部最優(yōu)；③參數(shù)R→確保以特定的速率減少HGS的搜索步驟，從而滿(mǎn)足前期大范圍探索目標(biāo)解空間和后期開(kāi)發(fā)目標(biāo)搜索深度的需要，在探索和開(kāi)發(fā)兩項(xiàng)相反功能之間具有較好的平衡性能；④HGS 可以根據(jù)最佳解（和正常解來(lái)進(jìn)化個(gè)體位置，以確保更多的探索模式和對(duì)解空間隱藏區(qū)域的更多覆蓋。

1.2 反演模型設(shè)計(jì)

對(duì)于含水層滿(mǎn)足泰斯模型的假定條件，在平面上無(wú)限延展、均質(zhì)各向同性的承壓含水層中的完整井進(jìn)行抽水試驗(yàn)，觀測(cè)井中的水位降深可用泰斯模型的解析解來(lái)表示。其模型為［10］：

式中：s 為水位降深，m；Q 為抽水井的流量，m3/d；T 為導(dǎo)水系數(shù)，m2/d；t為從抽水開(kāi)始持續(xù)的時(shí)間，d；W(u)為泰斯井函數(shù)；r為計(jì)算點(diǎn)到抽水井的距離，m；S 為含水層的儲(chǔ)水系數(shù)。計(jì)算W(u)時(shí)需要用廣義積分，計(jì)算過(guò)程將極大增加了算法的運(yùn)行時(shí)間。本文采用Rajesh Srivastava給出的W(u)近似表達(dá)式［20］：

當(dāng)u ≤1時(shí)，

當(dāng)u ≥1時(shí)，

式中：a0=-0.577 22；a1=-0.999 99；a2=-0.249 91；a3=0.055 19；a4=-0.009 76；a5=0.001 08；b0=0.267 77；b1=8.634 76；b2=18.059 02；b3=8.573 33；c0=3.958 50；c1=21.099 65；c2=25.632 96；c3=9.573 32。

因此，在優(yōu)化算法篩選參數(shù)過(guò)程中，可按照各時(shí)刻實(shí)測(cè)降深值與模擬降深值的離差平方和的均值構(gòu)建目標(biāo)適應(yīng)度函數(shù)：

式中：si為第i 個(gè)記錄點(diǎn)的實(shí)測(cè)水位降深，m；N 表示抽水試驗(yàn)的記錄點(diǎn)總數(shù)。

2 實(shí)例應(yīng)用

2.1 數(shù)據(jù)來(lái)源

為驗(yàn)證本文提出的方法的可靠性，采用2 組實(shí)例抽水試驗(yàn)數(shù)據(jù)進(jìn)行模型數(shù)值驗(yàn)證。實(shí)例1數(shù)據(jù)來(lái)源于山東省萊蕪縣西尚莊地區(qū)抽水試驗(yàn)［19］，選擇112孔為觀測(cè)孔，已知測(cè)量點(diǎn)距離抽水主井r=104 m，抽水量Q=5.221 9 m3/min，抽水時(shí)間t=840 min。實(shí)例2數(shù)據(jù)來(lái)源于文獻(xiàn)［10］，有一承壓含水層做流量抽水試驗(yàn)，假設(shè)該含水層的水文地質(zhì)條件滿(mǎn)足泰斯井流模型的假定條件。已知觀測(cè)孔距抽水井的距離r=100 m，主井做定流量抽水，抽水量Q=162.9 m3/min。利用HGS 算法反演參數(shù)（T，S），此外，試驗(yàn)同時(shí)采用近年來(lái)4種新型優(yōu)化算法黃金正弦算法（Gold-SA）、天鷹優(yōu)化算法（AO）、鯨魚(yú)優(yōu)化算法（WOA）以及阿基米德算法（AOA）對(duì)本例進(jìn)行含水層參數(shù)求解［14-17］。

根據(jù)大量實(shí)驗(yàn)分析，最大迭代次數(shù)設(shè)置太大導(dǎo)致迭代太長(zhǎng)消耗較多計(jì)算時(shí)間，設(shè)置太小迭代太短無(wú)法出現(xiàn)最優(yōu)解。種群取值過(guò)大會(huì)浪費(fèi)計(jì)算資源，取值太小不能夠提供足夠采樣點(diǎn)。當(dāng)搜索精度達(dá)到一定水平后，迭代次數(shù)和種群大小對(duì)結(jié)果的影響很小，值設(shè)置太大會(huì)增加計(jì)算時(shí)間。因此，測(cè)試參數(shù)設(shè)置標(biāo)準(zhǔn)如下：HGS 算法的允許最大迭代次數(shù)設(shè)為1000，種群大小設(shè)為30，根據(jù)文獻(xiàn)［18］將饑餓感下限LH設(shè)置為100，參數(shù)l設(shè)置為0.03。其他比較算法（包括Gold-SA、AO、WOA 和AOA）的種群大小與允許迭代次數(shù)設(shè)置均與HGS相同，其他算法特征參數(shù)設(shè)置參考原始文獻(xiàn)［14-17］。根據(jù)文獻(xiàn)［9］，實(shí)例1 中T，S 的搜索空間分別為（4.592×10-4，4.592×106）和（1×10-7，0.5），根據(jù)文獻(xiàn)［10］，實(shí)例2 中T，S 的搜索空間分別為（3，300）和（0.006，0.6）。為了測(cè)試各優(yōu)化算法在求解過(guò)程中的穩(wěn)定性，各方法獨(dú)立隨機(jī)運(yùn)行5次，統(tǒng)計(jì)包括最大值、最小值、平均值、中值和標(biāo)準(zhǔn)差在內(nèi)的五個(gè)項(xiàng)目結(jié)果，并記錄迭代完成時(shí)間，相關(guān)統(tǒng)計(jì)結(jié)果見(jiàn)表1和表2。

表1 實(shí)例1各方法獨(dú)立運(yùn)行5次的統(tǒng)計(jì)結(jié)果Tab.1 Statistical results of each method running independently for 5 times in Example 1

表2 實(shí)例2各方法獨(dú)立運(yùn)行5次的統(tǒng)計(jì)結(jié)果Tab.2 Statistical results of each method running independently for 5 times in Example 2

2.2 評(píng)價(jià)指標(biāo)

評(píng)價(jià)指標(biāo)采用平均相對(duì)誤差（MRE）、均方根誤差（RMSE）、平均絕對(duì)誤差（MAE）和納什效率系數(shù)（NSE），計(jì)算公式如下所示。以上4 種評(píng)價(jià)指標(biāo)可以給出不同方法反演過(guò)程的具體表現(xiàn)，例如納什效率系數(shù)越接近于1，表示模型總體質(zhì)量更好，可信度更高。

式中：s（i）是實(shí)測(cè)值；s*（i）是模型的模擬值；n是樣本數(shù)是實(shí)測(cè)水位的算數(shù)平均值。

2.3 結(jié)果分析

各方法的實(shí)測(cè)降深與理論降深對(duì)比圖見(jiàn)圖1和圖2，表3和表4記錄了5種方法所得的反演結(jié)果和評(píng)價(jià)指標(biāo)值。

圖1 實(shí)例1實(shí)測(cè)降深與模擬降深對(duì)比圖Fig.1 Comparison of measured and simulated depth drawdown in Example 1

圖2 實(shí)例2實(shí)測(cè)降深與模擬降深對(duì)比圖Fig.2 Comparison of measured and simulated depth drawdown in Example 2

從表1 和表2 中可以得到：①實(shí)例1 中，Gold-SA、AO、WOA和AOA 算法獲得的最小值分別低于最大值90.47%，92.36%，52.14%，99.16%，實(shí)例2 中，Gold-SA、AO、WOA 和AOA 算法獲得的最小值分別低于最大值98.21%，59.84%，36.99%，91.51%，HGS 算法運(yùn)行5 次，其標(biāo)準(zhǔn)差為0，說(shuō)明HGS 穩(wěn)健性強(qiáng)；②HGS的迭代時(shí)間相較于Gold-SA、AO、WOA 和AOA 算法分別縮短了94.43%，55.91%，27.07%，3.85% 和95.05%，63.16%，35.64%，16.29%，可以看出HGS具有較快的收斂速度。以上結(jié)果充分表明HGS 明顯優(yōu)于其余4 種算法，證明HGS 在求解地下水參數(shù)問(wèn)題的有效性和可行性。

根據(jù)圖1、2、表3 和表4 可以得到：①HGS 算法和其他4 種方法得到的含水層參數(shù)值十分接近，可認(rèn)為利用HGS算法求參的方法是有效可信的；②2 個(gè)實(shí)例結(jié)果均顯示HGS 算法的反演精度最高，適應(yīng)度分別為6.978×10-6和4.347×10-5，實(shí)例1中HGS算法較Gold-SA 算法、AO 算法、WOA 算法以及AOA 算法分別提高了99.85%，59.24%，39.16%，99.42%，對(duì)于實(shí)例2，HGS 算法較其余4 種算法分別提高了96.95%，65.91%，32.29%，99.10%。表明HGS算法具有較為可靠的收斂速度和全局尋優(yōu)效果；③基于平均相對(duì)誤差MRE、均方根誤差RMSE 和平均絕對(duì)誤差MAE指標(biāo)體系下各算法的評(píng)可以得價(jià)結(jié)果，將5 種方法的3 個(gè)誤差指標(biāo)值分別從小到大排列為：HGS＜WOA＜AO＜AOA＜Gold-SA以及HGS＜WOA＜AO＜Gold-SA＜AOA。在5 種方法中HGS 算法誤差指標(biāo)數(shù)值最小，排名第一，本文采取的HGS 算法具有一定的可行性并處于先進(jìn)水平；④用HGS 方法求得的Nash-Sutcliffe效率系數(shù)值在2組實(shí)例中都最接近1.000 0，從圖1和圖2中可以看出HGS算法預(yù)測(cè)曲線與實(shí)測(cè)曲線基本吻合，表明文中提出的方法反演得到的水文地質(zhì)參數(shù)符合相應(yīng)實(shí)際含水層的水力特性。

表3 實(shí)例1不同優(yōu)化方法參數(shù)求解結(jié)果的比較Tab.3 Comparison of parameter solution results of different optimization methods in Example 1

表4 實(shí)例2不同優(yōu)化方法參數(shù)求解結(jié)果的比較Tab.4 Comparison of parameter solution results of different optimization methods in Example 2

5種方法得到的結(jié)果存在差異，對(duì)于群智能算法而言，隨著迭代次數(shù)的增加，種群會(huì)不斷接近一個(gè)最優(yōu)個(gè)體區(qū)域，沒(méi)有確保整個(gè)解空間得到搜索，導(dǎo)致算法的多樣性喪失。例如，在WOA 中，最優(yōu)解是獵物的位置，而在位置更新公式中獵物的位置未被充分利用。并且在整個(gè)迭代過(guò)程中，最優(yōu)位置只有在出現(xiàn)優(yōu)于它的位置時(shí)才會(huì)更新，所以總的更新次數(shù)不多，導(dǎo)致算法搜尋效率不高。而在HGS算法中，因?yàn)橛袃蓚€(gè)饑餓權(quán)重的擾動(dòng)，在每次迭代后都對(duì)種群進(jìn)行一個(gè)干擾，使得其能很好地跳出局部最優(yōu)。在搜索階段，可以搜索盡可能完整的解空間，使算法能夠在早期達(dá)到快速收斂的效果。在開(kāi)發(fā)階段，可以在附近找到最優(yōu)解，保證了解的準(zhǔn)確性。

3 結(jié) 論

本文利用2 組實(shí)際抽水試驗(yàn)數(shù)據(jù)，測(cè)試一種新型群體智能優(yōu)化算法——饑餓搜索算法（HGS）來(lái)反演地下水含水層參數(shù)。將HGS 算法的計(jì)算結(jié)果與黃金正弦算法（Gold-SA）、天鷹優(yōu)化算法（AO）、鯨魚(yú)優(yōu)化算法（WOA）以及阿基米德算法（AOA）的計(jì)算結(jié)果對(duì)比，可得出以下結(jié)論：

（1）在利用算法進(jìn)行地下水參數(shù)反演的過(guò)程中，HGS 在較少有限次數(shù)內(nèi)即能快速找到全局最優(yōu)，具有較快的收斂速度，節(jié)省了計(jì)算資源。

（2）Gold-SA、AO、WOA 和AOA 算法的實(shí)驗(yàn)測(cè)試結(jié)果均劣于HGS 算法，且HGS 算法多次運(yùn)行結(jié)果的方差為零，彌補(bǔ)了以往采用群體智能優(yōu)化算法反演地下水參數(shù)結(jié)果不穩(wěn)定的不足。

（3）2 個(gè)實(shí)例中各方法的綜合性能從優(yōu)到劣排列分別為：HGS＞W(wǎng)OA＞AO＞AOA＞Gold-SA 和 HGS＞W(wǎng)OA＞AO＞Gold-SA＞AOA，驗(yàn)證了HGS算法在含水層參數(shù)反演的可行性并處于先進(jìn)水平，HGS可為水文地質(zhì)參數(shù)反演提供新的參考途徑。

（4）在后續(xù)開(kāi)展改進(jìn)群智能優(yōu)化算法工作中，算法開(kāi)發(fā)階段可以嘗試與饑餓驅(qū)動(dòng)搜索策略相結(jié)合，引入兩個(gè)基于動(dòng)態(tài)自適應(yīng)參數(shù)和可以平衡算法的全局和局部搜索能力，提高算法的尋優(yōu)性能，以解決更多優(yōu)化問(wèn)題。