張國(guó)新，陳 亮，李小雨，顏書(shū)法

（1.中電投蒙東能源集團(tuán)有限責(zé)任公司扎哈淖爾工業(yè)供水分公司，內(nèi)蒙古通遼 028000；2.河海大學(xué)巖土力學(xué)與堤埋工程教育部重點(diǎn)試驗(yàn)室，江蘇南京 210098；3.河海大學(xué)巖土工程科學(xué)研究所，江蘇南京 210098）

0 引言

據(jù)調(diào)查研究顯示，由滲透破壞造成堤防和大壩破壞的事故占大壩失事事故的40%以上[1]。滲漏會(huì)造成土體顆粒流失，如果任由滲漏發(fā)展會(huì)進(jìn)一步形成集中滲漏通道，將可能導(dǎo)致壩體塌陷或滑坡，甚至形成決堤，造成嚴(yán)重后果[2]。受堤壩規(guī)模、材料和地形、地質(zhì)條件的影響，堤壩往往存在多個(gè)集中滲漏通道，如果不能探測(cè)出全部滲漏通道，在主要通道進(jìn)行治理后，次要滲漏通道會(huì)逐漸發(fā)展成主要通道，繼續(xù)對(duì)堤壩造成危害。因此，有必要對(duì)堤壩進(jìn)行全面摸排，找到所有滲漏通道，進(jìn)行徹底加固。

目前，用于堤壩滲漏探測(cè)的技術(shù)主要有物探法[3]、同位素法[4]和溫度示蹤法[5-7]。物探法包括高密度電法[8]、瞬變電磁法[9]和地質(zhì)雷達(dá)法[10]。由于滲漏一般處在堤壩深層，在含水層中物探信號(hào)會(huì)受到強(qiáng)烈干擾，導(dǎo)致探測(cè)效果不明顯。同位素測(cè)定則成本昂貴，對(duì)環(huán)境也會(huì)造成一定的污染。溫度場(chǎng)屬于地層的天然物理場(chǎng)，通過(guò)測(cè)量地層不同部位的溫度值，可以得到整個(gè)地層連續(xù)的溫度場(chǎng)。因此，利用溫度場(chǎng)對(duì)滲漏通道進(jìn)行探測(cè)，具有成本低、效率高等優(yōu)點(diǎn)。

20世紀(jì)50年代，國(guó)外將地溫運(yùn)用于堤壩滲漏探測(cè)[11]，目前，溫度探測(cè)主要用于對(duì)大壩滲漏的定性分析，缺乏對(duì)滲漏通道的定量計(jì)算。本文基于已有的多個(gè)集中滲漏通道溫度場(chǎng)探測(cè)模型，建立反分析目標(biāo)函數(shù)，通過(guò)粒子群優(yōu)化算法[12]優(yōu)化目標(biāo)函數(shù)中的參數(shù)得到集中滲漏通道中心位置，并與同位素示蹤結(jié)果進(jìn)行對(duì)比，證明此優(yōu)化算法可以運(yùn)用在溫度示蹤確定大壩集中滲漏的中心位置。

1 多個(gè)集中滲漏通道反分析

1.1 溫度示蹤原理

滲流是土體中普遍存在的現(xiàn)象。滲流過(guò)程中流速緩慢，土和水之間有充分的時(shí)間進(jìn)行熱量交換，因此，在滲流過(guò)程中，巖土體和水的溫度是相同的。發(fā)生滲漏時(shí)，水流流速變快，巖土和滲漏水沒(méi)有時(shí)間進(jìn)行充分的熱量交換，滲漏水僅僅與附近土體發(fā)生部分熱量的交換。因此，滲漏附近的溫度場(chǎng)與滲流部分的溫度場(chǎng)存在明顯差異，離滲漏中心越近，這種差異越明顯，整個(gè)地層形成了具有一定特征的溫度場(chǎng)。

1.2 數(shù)學(xué)物理模型

正常滲流情況下，基于單個(gè)滲漏通道的一般熱傳導(dǎo)物理問(wèn)題的數(shù)學(xué)模型：

式中：M為滲漏通道到某點(diǎn)x的溫度函數(shù)；L0為滲漏通道的中心位置；T（x，t）為堤壩x測(cè)點(diǎn)t時(shí)刻的溫度值。

由于所采用的數(shù)學(xué)模型是線(xiàn)性的，溫度又是標(biāo)量值，運(yùn)用疊加原理可得n個(gè)集中滲漏通道下的一般熱傳導(dǎo)物理方程數(shù)學(xué)模型：

式中：n為集中滲漏通道的數(shù)量；Loi為第i個(gè)集中滲漏通道的位置；T(xj，t)為堤壩xj測(cè)點(diǎn)在n個(gè)集中滲漏通道作用下t時(shí)刻的溫度值。

在構(gòu)建數(shù)學(xué)模型時(shí)，將滲漏通道的中心作為坐標(biāo)原點(diǎn)，構(gòu)建了熱傳導(dǎo)方程。在實(shí)際測(cè)量中，并不知道滲漏的位置。因此，往往將某個(gè)探測(cè)孔的地面高程作為坐標(biāo)原點(diǎn)。想要運(yùn)用上述方程，就必須要把探測(cè)坐標(biāo)系統(tǒng)下測(cè)點(diǎn)的坐標(biāo)轉(zhuǎn)換為以滲漏通道中心為原點(diǎn)的坐標(biāo)。為此，采用轉(zhuǎn)換尺度無(wú)變化的Brusa-Wolf六參數(shù)模型，圖1是以滲漏通道中心為原點(diǎn)的坐標(biāo)系統(tǒng)O和探測(cè)坐標(biāo)系統(tǒng)O′之間相互關(guān)系示意圖。轉(zhuǎn)換公式如下：

圖1 模型坐標(biāo)O與探測(cè)坐標(biāo)O′轉(zhuǎn)換示意圖

式中：x，y，z為測(cè)點(diǎn)在熱傳導(dǎo)模型中的坐標(biāo)值；αx，αy，αz為歐勒角；x′，y′，z′為測(cè)點(diǎn)在探測(cè)坐標(biāo)系中的坐標(biāo)值；X0，Y0，Z0為熱傳導(dǎo)模型坐標(biāo)原點(diǎn)O在探測(cè)坐標(biāo)系中的坐標(biāo)值；矩陣R如下：

因此，可得在熱傳導(dǎo)模型坐標(biāo)系下，測(cè)點(diǎn)j在n個(gè)集中滲漏通道作用下的溫度：

式中：Tj為測(cè)點(diǎn)理論溫度值；Tji為測(cè)點(diǎn)j在集中滲漏通道i單獨(dú)作用下的溫度值；C1i，C2i為第i個(gè)滲漏通道產(chǎn)生的溫度場(chǎng)解析方程的未知參數(shù)；lji為測(cè)點(diǎn)j到集中滲漏通道i中心的距離

1.3 目標(biāo)函數(shù)的構(gòu)建

在實(shí)際工程中，通常會(huì)通過(guò)測(cè)量部分位置的溫度來(lái)計(jì)算式（2）中的參數(shù)。通過(guò)觀測(cè)數(shù)據(jù)求解數(shù)學(xué)物理模型中的參數(shù)問(wèn)題稱(chēng)為反問(wèn)題。但是，觀測(cè)過(guò)程中，由于人為偶然誤差和測(cè)量?jī)x器的系統(tǒng)誤差，觀測(cè)量一般可表示：

式中：為j點(diǎn)溫度觀測(cè)值；εj為誤差。

此處主要研究多個(gè)滲漏通道在同一平面同一深度的情況，因此，式（5）可以簡(jiǎn)化成：

式中：Xoi為第i個(gè)滲漏通道中心的橫軸坐標(biāo)值。

為此，最小二乘法最優(yōu)化目標(biāo)函數(shù)為

式中：m為測(cè)點(diǎn)個(gè)數(shù)。

2 粒子群優(yōu)化算法

粒子群優(yōu)化算法（Particle Swarm Optimization）是由KENNEDY等于1995年提出的一種模擬鳥(niǎo)類(lèi)覓食行為的智能算法。PSO利用鳥(niǎo)群中個(gè)體對(duì)信息的共享機(jī)制，使得整個(gè)鳥(niǎo)群朝食物最多的方向運(yùn)動(dòng)。PSO算法中有2個(gè)參數(shù)：速度和位置。通過(guò)每一次迭代后對(duì)適應(yīng)度函數(shù)進(jìn)行評(píng)價(jià)，不斷調(diào)整速度和位置，最終達(dá)到全局最優(yōu)值。具體表現(xiàn)見(jiàn)圖2。

圖2 粒子群優(yōu)化算法流程圖

每次迭代后速度更新公式：

式中：t為當(dāng)前迭代次數(shù)；w為慣性權(quán)重；vi（t）表為當(dāng)前粒子i的速度；c1，c2為學(xué)習(xí)因子；rand1和rand2為取值0～1之間的隨機(jī)數(shù)；pi（t）為當(dāng)前粒子i的個(gè)體極值；xi（t）為當(dāng)前粒子i的位置。

位置更新公式：

式中：xi（t+1）為第t+1次迭代開(kāi)始時(shí)粒子i的位置。

優(yōu)化過(guò)程中，為了避免陷入局部最小，動(dòng)態(tài)調(diào)整加速度因子和慣性權(quán)重。每次迭代后，w，c1和c2可表示：

式中：wmax為迭代過(guò)程中慣性權(quán)重的最大值，取wmax=1.4；wmin為迭代過(guò)程中慣性權(quán)重的最小值，取wmin=0.7；tmax為最大迭代次數(shù)，取tmax=500。

式中：c1max為迭代過(guò)程中學(xué)習(xí)因子的最大值，取c1max=2.5；c1min為迭代過(guò)程中慣性權(quán)重的最小值，取c1min=0.7。

文中每次迭代過(guò)程中c2取值與c1相同。

3 工程實(shí)例應(yīng)用

3.1 工程概況

某壩水庫(kù)工程等級(jí)為Ⅲ等，主要建筑物級(jí)別為3級(jí)，次要建筑物級(jí)別為4級(jí)。壩頂長(zhǎng)1 230.00 m，壩頂寬3.50 m，壩頂高程953.80 m，最大壩高21.60 m。上游壩坡采用模袋混凝土護(hù)坡，坡度為1∶3.5。下游壩坡為碎石護(hù)坡，在940.00 m處設(shè)有馬道，馬道寬2.50 m，馬道以上坡度為1∶3.0，以下為1∶3.5。河床覆蓋層截滲采用懸掛式混凝土心墻，心墻分布樁號(hào)為0-150～0+520和0+730～0+770，墻深25.00 m，墻厚1.20 m。

該水庫(kù)自2008年蓄水以來(lái)，出現(xiàn)滲漏主要分布在樁號(hào)0+200，0+250及樁號(hào)0+700～0+900附近，大部分區(qū)域壩后坡腳潮濕，后期對(duì)大壩進(jìn)行加固，加固后發(fā)現(xiàn)在下游排水溝仍然存在一股滲漏量較大的穩(wěn)定滲漏水。為此，在馬道布置16個(gè)勘探孔，從左到右孔號(hào)分別為ZK8到ZK23，孔深18.00 m，入基巖5.00 m，16個(gè)勘探孔一字型排開(kāi)與壩軸線(xiàn)平行。

3.2 溫度反演

鉆孔完成后，靜置鉆孔2 d，等鉆孔中的流速場(chǎng)和溫度場(chǎng)恢復(fù)到本底之后開(kāi)始對(duì)這16個(gè)鉆孔進(jìn)行溫度探測(cè)，進(jìn)行了一周探測(cè)后發(fā)現(xiàn)每個(gè)鉆孔的最低溫度在16.00 m的位置，這16個(gè)鉆孔的16.00 m溫度見(jiàn)圖3。

圖3 各孔最低溫度分布情況

從圖3中可以發(fā)現(xiàn)，在ZK17附近存在明顯的異常低溫區(qū)，因此，在ZK17附近肯定存在一條明顯集中滲漏通道。首先，式（7）中n=1，利用單集中滲漏通道熱傳導(dǎo)模型式進(jìn)行數(shù)據(jù)擬合，擬合的溫度曲線(xiàn)如圖4單個(gè)集中滲漏通道溫度擬合曲線(xiàn)所示，C1，C2，Xo分別為0.08，9.0和600。各點(diǎn)擬合殘差結(jié)果見(jiàn)表1。

表1 各孔溫度殘差分布情況

圖4 單個(gè)集中滲漏通道擬合曲線(xiàn)

從圖4中發(fā)現(xiàn)，實(shí)測(cè)數(shù)據(jù)點(diǎn)離散分布在擬合曲線(xiàn)的上下側(cè)，為了使各測(cè)點(diǎn)溫度對(duì)單個(gè)集中滲漏堤壩溫度曲線(xiàn)異常程度具有可比性，把擬合曲線(xiàn)整體向上移動(dòng)，移動(dòng)量為0.36。此時(shí)各點(diǎn)的殘差即為表1中修正后的殘差，殘差值越大，代表測(cè)點(diǎn)處越有可能存在其他滲漏通道。利用3次樣條曲線(xiàn)對(duì)修正后的殘差值作插值曲線(xiàn)，如圖5所示，由曲線(xiàn)的波谷數(shù)量確定所有滲漏通道的位置。圖5中共發(fā)現(xiàn)了5個(gè)波谷，其中波谷1經(jīng)過(guò)分析認(rèn)為是由于勘探范圍的限制，未能完全揭示可能存在的滲漏通道，這種情況求得的滲漏通道位置誤差較大，因此不考慮這部分的異常。其余4個(gè)波谷認(rèn)為是由于滲漏通道造成的，因此，除了在17號(hào)孔附近還存在另外3個(gè)滲漏通道。

圖5 3次樣條插值曲線(xiàn)

通過(guò)3次樣條插值曲線(xiàn)確定了4條滲漏通道，因此，在目標(biāo)函數(shù)中一共有12個(gè)參數(shù)：C11，C21，Xo1，C12，C22，Xo2，C13，C23，Xo3，C14，C24，Xo4。采 用 粒子群優(yōu)化算法初始化10個(gè)12維粒子，利用Python編制優(yōu)化算法進(jìn)行500次迭代后，目標(biāo)函數(shù)收斂，收斂過(guò)程見(jiàn)圖6，計(jì)算結(jié)果見(jiàn)表2。

表2 計(jì)算參數(shù)結(jié)果

圖6 優(yōu)化過(guò)程

將優(yōu)化后的參數(shù)帶入式（7），目標(biāo)函數(shù)方差為0.838，優(yōu)化曲線(xiàn)見(jiàn)圖7，各點(diǎn)殘差值見(jiàn)表1中優(yōu)化后的殘差一欄。

在圖7中明顯可以看到2條主要集中滲漏通道，其中心位置為600.00 m和910.00 m。此外,還有2條次要集中滲漏通道，它們的中心位置為517.47 m和664.60 m。

圖7 多通道溫度優(yōu)化曲線(xiàn)

3.3 同位素示蹤試驗(yàn)

為了進(jìn)一步說(shuō)明溫度示蹤法和粒子群優(yōu)化算法在滲漏探測(cè)中的可靠性，選擇同位素示蹤試驗(yàn)，測(cè)試8～23孔的流速，經(jīng)過(guò)測(cè)試，得到各孔16.00 m位置的流速見(jiàn)圖8。從圖8中可以發(fā)現(xiàn)：在ZK17和ZK22附近存在明顯較大流速，流速超過(guò)1.00×10-4cm/s，而在次要集中滲漏通道附近流速較小。ZK17的位置正好是600.00 m，所以，此處流速和溫度相符合，流速大，溫度低。ZK22離主要滲漏通道1還有10.00 m的距離，ZK22存在大流速，從側(cè)面說(shuō)明主要集中滲漏通道已經(jīng)發(fā)展到了此處。從圖7也可以看出主要滲漏通道1已經(jīng)包含了ZK22，而ZK15和ZK18流速正常，從側(cè)面說(shuō)明次要滲漏通道規(guī)模較小，還沒(méi)有形成大規(guī)模的集中滲漏通道。

圖8 各孔流速分布

4 結(jié)語(yǔ)

運(yùn)用粒子群優(yōu)化算法成功地反演出了內(nèi)蒙古某大壩的集中滲漏通道，豐富了溫度探測(cè)多滲漏通道位置方法的理論化和定量化；同位素示蹤法只能反應(yīng)局部滲漏情況，不能連續(xù)分析大壩滲流場(chǎng)，而溫度示蹤可以對(duì)大壩溫度場(chǎng)連續(xù)分析，為徹底探測(cè)堤壩滲漏和有效、便捷、經(jīng)濟(jì)整治滲漏提供了更加可行的辦法。