王文川，薛一銘，徐 雷

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

0 引 言

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

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

1 基本原理

1.1 算法原理

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

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

對所有位置的變化控制的E的計算公式如下：

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

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

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

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

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

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

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

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

1.2 反演模型設(shè)計

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

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

當(dāng)u ≤1時，

當(dāng)u ≥1時，

式中：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ù)過程中，可按照各時刻實測降深值與模擬降深值的離差平方和的均值構(gòu)建目標(biāo)適應(yīng)度函數(shù)：

式中：si為第i 個記錄點的實測水位降深，m；N 表示抽水試驗的記錄點總數(shù)。

2 實例應(yīng)用

2.1 數(shù)據(jù)來源

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

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

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

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

2.2 評價指標(biāo)

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

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

2.3 結(jié)果分析

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

圖1 實例1實測降深與模擬降深對比圖Fig.1 Comparison of measured and simulated depth drawdown in Example 1

圖2 實例2實測降深與模擬降深對比圖Fig.2 Comparison of measured and simulated depth drawdown in Example 2

從表1 和表2 中可以得到：①實例1 中，Gold-SA、AO、WOA和AOA 算法獲得的最小值分別低于最大值90.47%，92.36%，52.14%，99.16%，實例2 中，Gold-SA、AO、WOA 和AOA 算法獲得的最小值分別低于最大值98.21%，59.84%，36.99%，91.51%，HGS 算法運行5 次，其標(biāo)準(zhǔn)差為0，說明HGS 穩(wěn)健性強；②HGS的迭代時間相較于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ù)問題的有效性和可行性。

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

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

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

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

3 結(jié) 論

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

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

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

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

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