梁森，陳建華，李宏濤,羅威力，羅盈洲，艾姣姣，廖偉，

(1. 中建四局第一建設(shè)有限公司，廣東 廣州 510800； 2. 廣州大學(xué) 土木工程學(xué)院，廣東 廣州 510006)

0 引言

近年來(lái)，中國(guó)城市化進(jìn)程加快，基建工程項(xiàng)目越來(lái)越多。在巖溶地區(qū)修建建筑物易造成地基失穩(wěn)、地表塌陷等危害，且該危害往往具有突發(fā)性和不可預(yù)測(cè)性。因此，提前探測(cè)出地下溶洞并加以處理，對(duì)建筑物施工和運(yùn)營(yíng)安全具有重要意義。常用的地下溶洞探測(cè)方法有電阻率法[1-2]、瞬變電磁法[3-5]、地震波法[6-7]等。20世紀(jì)90年代后，日本學(xué)者提出了跨孔電阻率法[8-9]。跨孔電阻率法的反演算法常采用基于Tikhonov正則化的靈敏度迭代法，但研究發(fā)現(xiàn)該方法對(duì)初值、噪聲敏感，且容易陷入局部最優(yōu)，導(dǎo)致反演結(jié)果失真[10]。

群智能算法是通過(guò)模擬某種生物群體的某一行為而實(shí)現(xiàn)的一類基于群體搜索的智能優(yōu)化方法，它的優(yōu)點(diǎn)是全局搜索能力強(qiáng)、無(wú)需靈敏度信息、對(duì)初值不敏感。比如粒子群算法[11](PSO)、云模型果蠅算法(CMFOA)[12]、JAYA算法[13]、松鼠搜索算法(SSA)[14]等。本文通過(guò)數(shù)值算例對(duì)比分析，選出了4種算法中性能最優(yōu)異的算法——松鼠搜索算法，提出一種基于松鼠搜索算法的跨孔電阻率溶洞探測(cè)方法，以改善傳統(tǒng)算法的缺陷，提高反演精度。通過(guò)構(gòu)建小型、大型和串珠狀充水型溶洞這3種工況的地下溶洞模型，分析對(duì)比Tikhonov正則化的靈敏度迭代法和松鼠搜索算法在溶洞位置定位、溶洞輪廓、多解性等的反演效果；通過(guò)室內(nèi)物理模型和現(xiàn)場(chǎng)實(shí)驗(yàn)，驗(yàn)證松鼠搜索算法在跨孔電阻率溶洞探測(cè)領(lǐng)域中的實(shí)際應(yīng)用效果。

1 跨孔電阻率溶洞探測(cè)基本原理

1.1 基本原理

(1)

該優(yōu)化問(wèn)題為典型的非線性最小二乘優(yōu)化問(wèn)題，一般采用靈敏度迭代法進(jìn)行求解，即先對(duì)公式(1)的R(c)進(jìn)行Taylor展開(kāi)，從而得到線性化的優(yōu)化問(wèn)題，然后引進(jìn)Tikhonov正則化求解。該求解過(guò)程寫(xiě)成迭代形式為

(c-ck-1))‖2+λ‖c-ck-1‖2}=ck-1+

k=1,2,…,itermax。

(2)

其中u=[u1;u2;…;un]，包含了所有節(jié)點(diǎn)電勢(shì)；itermax為最大迭代次數(shù)。

由上可知，傳統(tǒng)的基于Tikhonov正則化的反演算法需要對(duì)初值集合c0進(jìn)行估計(jì)，也沒(méi)有脫離局部最優(yōu)解的數(shù)學(xué)手段。因此，該反演算法對(duì)測(cè)量噪聲敏感，容易陷入局部最優(yōu)解，最終導(dǎo)致反演結(jié)果與實(shí)際相差較大。

1.2 觀測(cè)模式

圖1 觀測(cè)模式示意Fig.1 A sketch of observation mode

2 松鼠搜索算法

2.1 方法原理

為了解決基于Tikhonov正則化的靈敏度迭代法存在的缺陷，本文引入松鼠搜索算法構(gòu)建跨孔電阻率溶洞探測(cè)反演方法以提高反演精度。松鼠搜索算法(SSA)是Mohit Jain于2018年提出的一種基于松鼠的動(dòng)態(tài)覓食行為的群體智能優(yōu)化算法[14]，該算法平衡了全局和局部搜索能力，在保證精度的情況下尋找全局最優(yōu)解。松鼠搜索算法的核心包括以下6個(gè)步驟。

1) 隨機(jī)初始化松鼠位置。假設(shè)森林中有NP個(gè)松鼠(FS)，其中，F(xiàn)Si,j表示第i只松鼠的第j維。森林中每只松鼠的初始位置由下式分配：

FSi=FSL+U(0,1)·(FSU-FSL) 。

(3)

式中：FSU和FSL是FSi的j個(gè)維度的上下限向量，U(0,1)是[0,1]范圍內(nèi)的隨機(jī)數(shù)。

2) 松鼠適應(yīng)性評(píng)估。由目標(biāo)函數(shù)計(jì)算出每只松鼠FS的適應(yīng)值，相應(yīng)的值存儲(chǔ)在數(shù)組f中：

(4)

式中：f(·)是目標(biāo)函數(shù)，用于計(jì)算每只松鼠的適應(yīng)值；n是松鼠FSi的維度；NP是種群數(shù)量。

4) 產(chǎn)生新松鼠位置。在松鼠動(dòng)態(tài)覓食過(guò)程中，可能會(huì)出現(xiàn)3種情況。

情況1：在橡樹(shù)上的松鼠可能會(huì)向山核桃樹(shù)移動(dòng)：

(5)

情況2：正常樹(shù)上的松鼠可能會(huì)向橡子樹(shù)移動(dòng)，以滿足它們的日常能量需求：

(6)

情況3：一些松鼠在正常的樹(shù)上已經(jīng)吃了橡子堅(jiān)果，它們可能會(huì)向山胡桃樹(shù)移動(dòng)，以便儲(chǔ)存山胡桃堅(jiān)果，以便在食物短缺時(shí)食用：

(7)

5) 判斷季節(jié)性監(jiān)測(cè)條件。季節(jié)變化顯著影響松鼠的覓食活動(dòng)，與秋天相比，氣候條件迫使它們?cè)诙觳荒敲椿钴S。在算法中引入了一個(gè)季節(jié)性的監(jiān)測(cè)條件，可以防止算法陷入局部最優(yōu)解。模擬行為步驟如下。

t=1,2,3,…,itermax。

(8)

c.如果季節(jié)監(jiān)測(cè)情況屬實(shí)，將無(wú)法在森林中尋找最佳冬季食物來(lái)源的松鼠隨機(jī)遷移。

6) 冬季末隨意搬遷。松鼠在冬季無(wú)法在森林中尋找最佳食物來(lái)源，但仍然存活下來(lái)，它們可能會(huì)向新的方向覓食。在建模中引入這種行為可以增強(qiáng)算法的探索能力，可同過(guò)方程式模擬：

(9)

式中的Lévy(n)是研究人員用來(lái)提高各種元啟發(fā)式算法的全局搜索能力的一種強(qiáng)大的數(shù)學(xué)工具[15]。

2.2 基于松鼠搜索算法的反演

(10)

FS=c=[c1,c2,…,cN]T。

式中：f為目標(biāo)函數(shù)值即松鼠的適應(yīng)值；s為電極方案數(shù)量；FS表示松鼠位置。在本文中，停止準(zhǔn)則設(shè)置為最大迭代次數(shù)iter≥itermax。圖2為采用松鼠搜索算法進(jìn)行溶洞反演的流程。

圖2 松鼠搜索算法反演流程Fig.2 An inversion flowchart of the squirrel search algorithm

3 數(shù)值算例

3.1 不同智能算法計(jì)算結(jié)果對(duì)比

為了檢驗(yàn)松鼠搜索算法在跨孔電阻率溶洞探測(cè)的算法優(yōu)勢(shì)，本文選取了3種常用的群智能優(yōu)化算法進(jìn)行對(duì)比，即：云模型果蠅算法(CMFOA)、粒子群算法(PSO)和JAYA算法。采用二維的四節(jié)點(diǎn)矩形單元進(jìn)行網(wǎng)格劃分，并采用Matlab編寫(xiě)程序和計(jì)算。圖3中紅色框線是溶洞真實(shí)位置輪廓線。模型參數(shù)為：孔間距5m，孔深度10m，單孔電極10個(gè)，圍巖電阻率263Ω·m，溶洞電阻率10Ω·m，溶洞坐標(biāo)為(1.75m，2.75m)。所有群智能算法的參數(shù)設(shè)置一致，最大迭代數(shù)量itermax=100，種群數(shù)量NP=250，最小電阻率值ρmin=10Ω·m，最大電阻率值ρmaxρmax=263Ω·m。

圖3、圖4對(duì)比了4種群智能優(yōu)化算法的反演計(jì)算結(jié)果。圖3顯示，只有松鼠搜索算法(SSA)可以精確地反演出溶洞的位置和大小，其他算法均不能反演出溶洞真實(shí)位置。從圖4可以看出，云模型果蠅算法(CMFOA)和粒子群算法(PSO)雖均已收斂，但都陷入局部最優(yōu)解；JAYA算法雖然未陷入局部解，但收斂速度弱于松鼠搜索算法；松鼠搜索算法的收斂速度快，最終的適應(yīng)值趨近于0(即趨近于精確解)。綜上，相對(duì)于其他3種常用的群智能優(yōu)化算法，松鼠搜索算法收斂速度快且可以跳出局部最優(yōu)解，精確度高。

圖3 不同智能優(yōu)化算法反演結(jié)果Fig.3 Inversion diagrams of different intelligent optimization algorithms

圖4 不同智能優(yōu)化算法的反演迭代結(jié)果Fig.4 Inversion iterations of different intelligent optimization algorithms

3.2 松鼠搜索算法與靈敏度迭代法計(jì)算結(jié)果對(duì)比

為了對(duì)比松鼠搜索算法與靈敏度迭代法的反演精確度，設(shè)計(jì)了3種數(shù)值算例，分別是小溶洞模型、大溶洞模型和串珠狀溶洞模型(小溶洞尺寸為2m×3m，大溶洞尺寸為3m×4m，串珠狀溶洞為小溶洞和大溶洞呈串珠狀排列的兩溶洞模型)。模擬區(qū)域范圍為20m×10m，單孔布置20個(gè)電極，圍巖電阻率為263Ω·m，溶洞電阻率為10Ω·m。松鼠搜索算法初始參數(shù)設(shè)置為：最大迭代數(shù)量itermax=100，種群數(shù)量NP=250，最小電阻率值ρmin=10Ω·m，最大電阻率值ρmax=263Ω·m。

圖5給出了3種溶洞(圖中紅色線框)的靈敏度迭代法反演結(jié)果?？梢钥闯觯`敏度迭代法可以清晰地反演出小溶洞和大溶洞的形態(tài)和位置，其中大溶洞反演結(jié)果在電極附近顯示較多假異常體；在反演串珠狀溶洞時(shí)，探測(cè)區(qū)域的中部出現(xiàn)了一個(gè)大塊假異常體，探測(cè)結(jié)果顯示了3個(gè)獨(dú)立的溶洞。靈敏度迭代法的反演結(jié)果呈現(xiàn)出了多解性，嚴(yán)重影響了溶洞的準(zhǔn)確定位和形態(tài)的辨別。

圖5 靈敏度迭代法反演結(jié)果Fig.5 The inversion results of the sensitivity iterative method

圖6為3種溶洞的松鼠搜索算法反演結(jié)果。從圖中可以看出：對(duì)3種溶洞，反演結(jié)果均可以精確地定位出溶洞位置，其形態(tài)輪廓基本對(duì)應(yīng)，并且?guī)缀鯖](méi)有大塊的假異常體。這表明反演結(jié)果沒(méi)有多解性，可以提供精確的溶洞識(shí)別結(jié)果。

圖6 松鼠搜索算法反演結(jié)果Fig.6 The inversion results of the squirrel search algorithm

圖7為松鼠搜索算法的迭代曲線，可以看出算法基本在第10步左右就已經(jīng)收斂，收斂速度極快且精確度高。

圖7 松鼠搜索算法迭代曲線Fig.7 Iterative graphs of the squirrel search algorithm

表1總結(jié)了松鼠搜索算法和靈敏度迭代法的反演結(jié)果。整體上看，靈敏度迭代法雖可以很好地定位小溶洞和大溶洞的位置,卻無(wú)法準(zhǔn)確定位串珠狀溶洞，而松鼠搜索算法可以精準(zhǔn)地定位出3種溶洞的位置。在溶洞輪廓反演方面，靈敏度迭代法雖可以較好地反演小和大溶洞的輪廓卻無(wú)法很好地反演出串珠狀溶洞的真實(shí)輪廓，而松鼠搜索算法可以精準(zhǔn)地反演出3種情況下溶洞的真實(shí)輪廓。在反演多解性方面，靈敏度迭代法的反演結(jié)果在測(cè)量電極附近顯示出較多假異常體，尤其在反演串珠狀溶洞時(shí)，探測(cè)區(qū)域的中部出現(xiàn)了大塊假異常體，嚴(yán)重影響溶洞的準(zhǔn)確定位；與之相比，松鼠搜索算法可以有效抑制反演的多解性，反演結(jié)果精確，假異常體少。

表1 反演結(jié)果對(duì)比Table 1 Contrast of inversion results

綜上所述，無(wú)論是溶洞位置定位還是溶洞輪廓反演的表現(xiàn)，松鼠搜索算法均優(yōu)于靈敏度迭代法。同時(shí)靈敏度迭代法的反演結(jié)果呈現(xiàn)多解性，而松鼠搜索算法則能有效抑制多解性，使得反演結(jié)果的精確度顯著提升。

4 實(shí)驗(yàn)驗(yàn)證

室內(nèi)模型實(shí)驗(yàn)采用國(guó)產(chǎn)集中式高密度電法探測(cè)儀，主要由多功能直流電法儀、電極轉(zhuǎn)換器、高密度電纜、電源箱和電源組成(圖8)。

圖8 高密度電法探測(cè)儀Fig.8 High-density electric detector

圖9為制作的室內(nèi)溶洞物理模型，原型與模型的幾何相似比G為10。表2為各幾何因素原型尺寸與模型尺寸的對(duì)照。基于相似性原理，模型需要滿足以下兩個(gè)方面的要求：① 幾何因素比值統(tǒng)一為10，幾何因素包括地質(zhì)體的大小、埋深、位置、電極位置等；② 各電性不均勻體的比值參數(shù)須與實(shí)際地質(zhì)條件一致。采用粉質(zhì)黏土模擬灰?guī)r，其電阻率約為263Ω·m；采用邊長(zhǎng)為0.1m的六面體空心導(dǎo)電鐵塊模擬充水溶洞，20 ℃(常溫)時(shí)的電阻率為78×10-8Ω·m；采用邊長(zhǎng)為0.2m 的六面體混凝土塊模擬充氣溶洞，20 ℃(常溫)干燥時(shí)的電阻率約為 1000Ω·m。

圖9 溶洞模型Fig.9 Cave models

表2 各幾何因素原型尺寸與模型尺寸對(duì)照Table 2 The comparison of the prototype dimensions of each geometric factor with the size of the model

本次實(shí)驗(yàn)主要采集電流數(shù)據(jù)和電壓數(shù)據(jù)。按新型的四極觀測(cè)模式跑極并記錄電勢(shì)差數(shù)據(jù)，將每種A-B組合下測(cè)量到的10次電流值加和后取平均值，該平均值即為同種輸入輸出A-B組合下的電流值。由于儀器的測(cè)量結(jié)果只能反映出一個(gè)正值電流，實(shí)際上輸入輸出的電流值應(yīng)互為相反數(shù)，故將上述的平均值取相反數(shù)，其相反數(shù)即為輸出的電流；另一方面，也是為了確保邊值問(wèn)題上解的存在性，所有輸入輸出的電流總和必須為0。電勢(shì)差是每個(gè)測(cè)點(diǎn)相對(duì)于中間測(cè)點(diǎn)N的電勢(shì)差，在進(jìn)行反演計(jì)算時(shí)，需要在Matlab上模擬出中間測(cè)點(diǎn)N的電勢(shì)值，然后用測(cè)量的電勢(shì)差值加上N點(diǎn)的電勢(shì)值即可近似為每個(gè)測(cè)點(diǎn)的電勢(shì)。有限元方法下邊界條件的方程為

(11)

對(duì)電勢(shì)值具體的處理方式為：同種輸入輸出組合下的電勢(shì)值加和并求平均值，然后用每個(gè)電勢(shì)值減去所求的平均值。

模擬區(qū)域設(shè)置為1.5m×1.0m；松鼠搜索算法參數(shù)設(shè)置為itermax=500；充水溶洞模型參數(shù)設(shè)置為：NP=250，ρmin=0.1Ω·m，ρmax=263Ω·m；充氣溶洞模型參數(shù)設(shè)置為：NP=250，ρmin=263Ω·m，ρmax=1000Ω·m。將處理好的電流和電勢(shì)差數(shù)據(jù)導(dǎo)入松鼠搜索算法的跨孔電阻率探測(cè)反演程序進(jìn)行反演計(jì)算，通過(guò)迭代反演得到最終的反演成像結(jié)果。

圖10為充水溶洞模型的反演結(jié)果和迭代曲線，圖中紅色方框表示導(dǎo)電鐵塊的實(shí)際埋設(shè)位置。反演結(jié)果顯示：在框線內(nèi)有一個(gè)大小為0.05m×0.1m、電阻率為0.1Ω·m的低阻體，此外，在框線外有少量的假異常體。從迭代曲線圖可以看出，適應(yīng)值在第200次迭代反演計(jì)算中趨于0，說(shuō)明反演結(jié)果已收斂，收斂速度較快。反演結(jié)果基本與實(shí)際情況相符。

圖10 充水溶洞模型反演結(jié)果(a)與迭代曲線(b)Fig.10 Inversion results of water-filled karst cave model (a) and iterative curve (b)

圖11為充氣溶洞模型的反演結(jié)果和迭代曲線，圖中紅色方框表示混凝土塊的實(shí)際埋設(shè)位置。圖11a顯示框線內(nèi)有一個(gè)0.05m×0.05m、電阻率值為1000Ω·m的高阻體，框線外有少量的假異常體。圖11b與圖10b同樣，適應(yīng)值也在第200次迭代反演計(jì)算中趨于0。

圖11 充氣溶洞模型反演結(jié)果(a)與迭代曲線(b)Fig.11 Inversion results of the gas-filled cave model (a) and iterative curve (b)

綜上，不管是充水還是充氣的溶洞模型，在考慮有測(cè)量誤差、模型誤差和背景土體非均質(zhì)性的影響下，反演結(jié)果都基本反映了室內(nèi)模型的真實(shí)情況。

5 實(shí)際工程應(yīng)用

為了進(jìn)一步檢驗(yàn)改進(jìn)觀測(cè)模式和反演方法的有效性，在廣州市花都區(qū)花山鎮(zhèn)平西村、新和村美華航空電子研發(fā)項(xiàng)目工程建設(shè)工地開(kāi)展了溶洞探測(cè)現(xiàn)場(chǎng)實(shí)驗(yàn)。根據(jù)鉆探揭露，施工場(chǎng)地下伏石炭系灰?guī)r，大部分地段溶洞極發(fā)育，勘察中鉆遇溶洞和土洞的鉆孔有47個(gè)，見(jiàn)洞率達(dá)32.6 %，見(jiàn)溶洞層數(shù)為1～2層，為巖溶強(qiáng)發(fā)育區(qū)；大部分溶洞頂板較薄，多為串珠狀，半充填—全充填，部分為空洞，充填物為軟塑狀或流塑狀黏性土及少量風(fēng)化巖塊。

根據(jù)地質(zhì)調(diào)繪、物探和鉆探資料，測(cè)區(qū)內(nèi)大致有6類不同介質(zhì)的地層：素填土、粉質(zhì)黏土、礫砂、溶洞、中風(fēng)化石灰?guī)r和微風(fēng)化石灰?guī)r，地電斷面大致分為5層，詳見(jiàn)表3。溶洞與圍巖存在較大的電性差異，具備直流電法勘探的物性前提。溶洞、軟弱夾層、溶溝、富水帶等強(qiáng)溶蝕帶為此次探測(cè)的目標(biāo)體。

表3 測(cè)區(qū)視電阻率分布Table 3 The visual resistivity distribution table of the measured area

現(xiàn)場(chǎng)實(shí)驗(yàn)選取了2個(gè)鉆探好的孔位，兩孔相距約10m，井口高程相差0.25 m,鉆探深度分別為23 m和20 m(表4)。分析鉆探結(jié)果可以推知，在鉆孔深度13.2～18m處存在溶洞(充填粉質(zhì)黏土)，圍巖為微/中風(fēng)化石灰?guī)r。兩孔之間的具體溶洞分布情況未知。

表4 鉆孔柱狀圖部分結(jié)果Table 4 Summary of the results of the drilling bar chart section

根據(jù)新型的四極觀測(cè)模式測(cè)量電勢(shì)差數(shù)據(jù)，將得到的電勢(shì)差數(shù)據(jù)導(dǎo)入反演程序中，模擬區(qū)域設(shè)置為20m×12m，單孔電極數(shù)量為16個(gè)。松鼠搜索算法參數(shù)設(shè)置為:itermax=1000，NP=250，ρmin=1Ω·m，ρmax=3000Ω·m。

鉆孔剖面的基于松鼠搜索算法的反演成像結(jié)果見(jiàn)圖12。從迭代曲線圖(圖12d)可以看出，適應(yīng)值在第200次迭代反演計(jì)算中趨于0，說(shuō)明反演結(jié)果已收斂，收斂速度較快。圖12a顯示：在深度0～3m范圍電阻率較高，可以推斷屬于素填土層，含水率較低所以電阻率較高；在深度3～12m范圍電阻率值大幅降低，推測(cè)此深度范圍位于地下水位或發(fā)生土層電性改變，引起電阻率升高，部分區(qū)域的電阻率較高推測(cè)可能是孤石引起；在深度12m以下的探測(cè)范圍內(nèi)電阻率普遍升高，可以推測(cè)到達(dá)了石灰?guī)r層；在深度14～18m范圍內(nèi)，出現(xiàn)了3個(gè)電阻率值與上層土相近的孔洞，推測(cè)石灰?guī)r層出現(xiàn)了串珠狀的溶洞，充填物為上層土(粉質(zhì)黏土)，這與2個(gè)孔位的鉆探結(jié)果顯示的深度13.2～18m處有串珠狀溶洞基本相符。

由于剖面中部的溶洞是否存在未知，為進(jìn)一步驗(yàn)證算法反演結(jié)果，在兩孔中間打一個(gè)驗(yàn)證孔，其柱狀圖如圖12b所示。將反演圖與驗(yàn)證柱狀圖對(duì)比發(fā)現(xiàn)，在剖面中部位置確實(shí)存在溶洞，且溶洞位置與反演結(jié)果預(yù)報(bào)位置基本吻合。

圖12 現(xiàn)場(chǎng)實(shí)驗(yàn)反演結(jié)果Fig.12 Inversion results of field experiment

綜上所述，基于松鼠搜索算法的反演結(jié)果與鉆探結(jié)果基本吻合。

6 結(jié)論

1) 數(shù)值算例的反演結(jié)果表明：對(duì)比粒子群算法(PSO)、云模型果蠅算法(CMFOA)和JAYA算法，基于松鼠搜索算法的跨孔電阻率溶洞探測(cè)方法具有收斂速度快、精確度高、全局優(yōu)化能力強(qiáng)的優(yōu)勢(shì)。

2) 數(shù)值模擬小溶洞、大溶洞和串珠狀溶洞的反演結(jié)果表明：在溶洞位置定位、溶洞輪廓識(shí)別、多解性抑制方面，松鼠搜索算法均優(yōu)于靈敏度迭代法。

3) 通過(guò)室內(nèi)溶洞模型試驗(yàn)可知，基于松鼠搜索算法的跨孔電阻率溶洞探測(cè)反演結(jié)果基本反映了室內(nèi)模型的真實(shí)地下情況?，F(xiàn)場(chǎng)實(shí)驗(yàn)最終的反演結(jié)果與鉆探結(jié)果對(duì)比分析，反演結(jié)果基本可以準(zhǔn)確地反演出地下溶洞的真實(shí)位置和形態(tài)，進(jìn)一步驗(yàn)證了基于松鼠搜索算法的反演方法實(shí)際工程應(yīng)用效果。