倪紹虎，田 偉，趙毅鋒，魏盛夏子，趙瑞存，鄭海圣

（1.中國電建集團(tuán)華東勘測設(shè)計研究院有限公司,浙江杭州，311122；2.浙江省抽水蓄能工程技術(shù)研究中心，浙江杭州，311122；3.華東天荒坪抽水蓄能有限責(zé)任公司，浙江湖州，313300）

0 引言

地下水的影響是地下工程設(shè)計中必須考慮的關(guān)鍵環(huán)節(jié)之一。一方面地下水通過物理力學(xué)作用軟化巖體，改變其物理力學(xué)特性，另一方面也是結(jié)構(gòu)主要荷載，直接關(guān)系到工程穩(wěn)定及結(jié)構(gòu)安全，滲流作用下隧洞突水突泥、內(nèi)水外滲發(fā)生滲漏、外水滲漏等諸多問題，均為地下水滲流引起。因此，當(dāng)?shù)叵滤^豐富時，為全面了解地下水分布及其對工程的影響，通常需要進(jìn)行三維滲流場分析[1-2]。

地下水滲流分析最為關(guān)鍵的環(huán)節(jié)是通過大量的工程地質(zhì)、水文地質(zhì)探查、勘測，查明地下水的發(fā)育程度、空間分布特征及區(qū)域巖體的滲透性。通常采用現(xiàn)場地質(zhì)鉆孔進(jìn)行原位試驗，這些現(xiàn)場試驗獲得的結(jié)果能較好地反映試驗區(qū)滲透特性。然而巖土介質(zhì)是十分復(fù)雜的非均一介質(zhì)，試驗成果往往離散性較強(qiáng)，代表性較差。同時原位試驗也需要耗費(fèi)大量人力和物力，因此，利用原位觀測、監(jiān)測獲取的地下水?dāng)?shù)據(jù)信息進(jìn)行滲流場反分析，反演獲得各滲透分區(qū)滲透系數(shù)、滲流場分布、涌水量、滲漏量等指標(biāo)，具有很好的工程實(shí)用性和應(yīng)用價值。

反演分析方法已廣泛應(yīng)用于工程滲流分析中。朱岳明等[3]通過優(yōu)化算法反演分析了巖體滲透系數(shù)張量。劉先珊等[4]提出模擬退火算法神經(jīng)網(wǎng)絡(luò)，對大壩滲流場進(jìn)行了反分析。盛金昌等[5]運(yùn)用序列二次優(yōu)化方法，反演了工程區(qū)域的天然滲流場。謝紅強(qiáng)等[6]結(jié)合工程區(qū)域水文地質(zhì)資料，反演分析了巖體各向異性滲流場。張璇等[7]采用正交設(shè)計參數(shù)樣本，利用遺傳神經(jīng)網(wǎng)絡(luò)反演了工程區(qū)域的初始滲流場。段斌等[8]采用二次規(guī)劃反演，對天然滲流場進(jìn)行了反演分析。魏進(jìn)兵等[9]采用BP網(wǎng)絡(luò)對滑坡體進(jìn)行了非飽和滲流反演。

然而巖土工程反演多維多極值問題的求解是一個極其復(fù)雜的問題。為了得到理想可靠的滲流反演結(jié)果，對軟、硬件兩方面都有相當(dāng)高的要求，需要一種高效適用的數(shù)學(xué)或仿生智能優(yōu)化方法，能夠?qū)崿F(xiàn)對多參數(shù)多維問題的多目標(biāo)優(yōu)化快速搜索?；诖?，運(yùn)用有限單元法進(jìn)行三維滲流場分析，提出采用粒子群智能優(yōu)化算法進(jìn)行復(fù)雜三維滲流場優(yōu)化反演分析，顯著提高反演計算效率和可靠性，可有效應(yīng)用于地下巖土工程滲流分析、滲流控制及水文地質(zhì)條件評價中。

1 三維滲流場優(yōu)化反演分析

1.1 等效連續(xù)介質(zhì)滲流分析的基本前提

嚴(yán)格來講，巖體是非連續(xù)介質(zhì)，其滲流問題多屬于裂隙滲流問題，完整的巖塊滲透性極低，巖體滲流主要受其內(nèi)部結(jié)構(gòu)面和微裂隙控制。但裂隙網(wǎng)絡(luò)滲流比較復(fù)雜，裂隙網(wǎng)絡(luò)分布及參數(shù)基本不可能準(zhǔn)確獲得。因此，通常在滿足一定條件時，可采用等效連續(xù)介質(zhì)的方法將裂隙巖體滲流概化為連續(xù)介質(zhì)模型。對于巖體滲透性表征單元體存在且相對計算模型很小時，采用等效連續(xù)介質(zhì)滲流模型是可行的、高效的[10-11]，同時也完全能滿足工程需要。

滲流問題能否視為等效連續(xù)介質(zhì)滲流問題的關(guān)鍵在于分析域的滲流表征單元體積，即尺寸效應(yīng)。尺寸效應(yīng)是裂隙巖體的典型特征，如圖1所示。巖體中的微裂紋、節(jié)理、裂隙等不連續(xù)結(jié)構(gòu)面的存在，使不同試件尺寸以及不同結(jié)構(gòu)面發(fā)育程度的巖體的特性各異。不同研究范疇的裂隙巖體的表征單元體積不同，其物理力學(xué)特性和滲透特性有不同的表征單元體積，分別代表巖體力學(xué)特性的尺寸效應(yīng)和滲透特性的尺寸效應(yīng)，稱滲透特性的表征單元體積為滲透表征單元體積（permeability representative element volume，PREV），如圖2所示。土體等孔隙介質(zhì)中固體顆粒間的孔隙很小，其PREV非常小，滲流分析時通常可作為連續(xù)介質(zhì)考慮。而裂隙巖體不同，其空隙結(jié)構(gòu)主要為孔隙較大的裂隙網(wǎng)絡(luò)，其PREV一般較大甚至不存在，因此，在多數(shù)情況下不能等效為連續(xù)介質(zhì)，而應(yīng)按照非連續(xù)介質(zhì)考慮，即使在某些情況下可以按等效連續(xù)介質(zhì)考慮，也必須考慮其滲透各向異性。

WEI Z Q等人[12-13]提出了裂隙巖體PREV估算方法，當(dāng)模型尺寸大于10倍任意節(jié)理的隙間距St時，裂隙巖體的滲透性趨于穩(wěn)定；當(dāng)模型尺寸大于50倍任意節(jié)理的隙間距St時，巖體的滲透性變化很小。單組裂隙PREV可表示為：

圖1 裂隙巖體尺寸效應(yīng)示意圖Fig.1 Size effect of the fractured rock

圖2 裂隙巖體滲透表征單元體積示意圖Fig.2 Volume of the permeability representative element of the fractured rock mass

式中，St為單組裂隙平均間距。對于幾何特性和分布不同的多組裂隙，可用多組裂隙的平均間距Sa0代替St，并可表示為：

在裂隙巖體PREV存在且遠(yuǎn)小于研究區(qū)域尺寸時，可以等效為連續(xù)介質(zhì)考慮，通過裂隙巖體的滲透張量來描述裂隙巖體的滲透各向異性及裂隙網(wǎng)絡(luò)的水力學(xué)行為。若表征單元體不存在或過大，可考慮將研究區(qū)域進(jìn)行滲透性分區(qū)，再確定分區(qū)巖體的滲透性表征單元體。若表征單元體確實(shí)不存在或太大，則不能采用等效連續(xù)介質(zhì)模型。

1.2 滲流有限單元法基本原理

采用等效連續(xù)介質(zhì)滲流模型進(jìn)行三維穩(wěn)定滲流場反演分析，計算方法采用有限單元法。對于穩(wěn)態(tài)滲流問題，其有限單元法基本方程可表示為：

式中，[K]為滲透傳導(dǎo)矩陣；{h}為節(jié)點(diǎn)勢列陣；{F}為節(jié)點(diǎn)流量列陣。

1.3 滲流優(yōu)化反演分析方法

巖體滲流反分析通常分為直接法和間接法兩大類。直接法要求實(shí)測數(shù)據(jù)較為完備且準(zhǔn)確性較高，應(yīng)用較少。間接法則利用正問題解的適定性，通過數(shù)學(xué)方法將計算值逐步逼近觀測值，運(yùn)用廣泛。優(yōu)化反演時，研究區(qū)域內(nèi)測點(diǎn)的計算值與實(shí)測值的誤差，表示為：

式中：n為測點(diǎn)總數(shù)；xi(i=1,2,...,n)為一組反映巖體滲流場的水文地質(zhì)參數(shù)，如導(dǎo)水系數(shù)、給水度等多種初始參數(shù)；ωj為第 j測點(diǎn)的權(quán)重因子為第 j測點(diǎn)處的水頭（流量）計算值；為第 j測點(diǎn)處的現(xiàn)場水頭（流量）監(jiān)測值。

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

優(yōu)化反演需借助數(shù)學(xué)方法，即優(yōu)化算法，采用改進(jìn)的粒子群優(yōu)化算法進(jìn)行滲流場優(yōu)化反演分析。算法基于種群的全局搜索策略，根據(jù)粒子自身速度指導(dǎo)優(yōu)化搜索，粒子群根據(jù)記憶調(diào)整搜索策略，以實(shí)現(xiàn)路徑最優(yōu)。

粒子群算法中，各粒子均是搜索空間中的一個解。采用“群體”和“進(jìn)化”的概念，依據(jù)粒子個體的適應(yīng)值大小進(jìn)行操作。根據(jù)粒子自身的飛行經(jīng)驗和同伴粒子的飛行經(jīng)驗，調(diào)整粒子飛行策略和路徑。每個粒子通過飛行過程的最好位置以及粒子群體經(jīng)歷過的最好位置，不斷更新自己調(diào)整飛行軌跡，從而產(chǎn)生新一代群體，再進(jìn)入下一次飛行。對于第k次迭代，每一個粒子按下式更新：

式中：i=1,2,…,M，M為群體中粒子的總數(shù)；w為慣性權(quán)重系數(shù)；為第k次迭代粒子i飛行速度矢量的第d維分量；為第k次迭代粒子i位置矢量的第d維分量；pid為粒子i個體最好位置的第d維分量；pgd為群體最好位置的第d維分量；c1、c2為學(xué)習(xí)因子，通常取[0,4]的兩個定常數(shù)，推薦取為c1=c2=2.0；r1、r2為[0,1]的相互獨(dú)立的隨機(jī)數(shù)。

3 工程運(yùn)用

將上文所述的粒子群優(yōu)化算法運(yùn)用于浙江省天荒坪一級抽水蓄能電站長期運(yùn)行期的三維滲流場優(yōu)化反演分析中。

3.1 工程概況及水文地質(zhì)條件

天荒坪抽水蓄能電站位于浙江省安吉縣天荒坪鎮(zhèn)境內(nèi)，距離杭州57 km，距離上海175 km，接近華東電網(wǎng)的負(fù)荷中心。電站樞紐主要由上水庫、下水庫、輸水系統(tǒng)、地下廠房洞室群和開關(guān)站組成，安裝6臺300 MW的抽水蓄能機(jī)組，總裝機(jī)容量1 800 MW。

輸水系統(tǒng)由上水庫進(jìn)/出水口、斜井、混凝土岔管、鋼支管、尾水隧洞及下水庫進(jìn)/出水口等組成，引水隧洞采用一洞三機(jī)布置型式，尾水隧洞采用一洞一機(jī)布置型式。

輸水系統(tǒng)皆處于地下水位線以下。由于沿線地形陡峻，地下水補(bǔ)給條件差～較差，巖體完整～較完整，賦水性差，開挖過程中初始地下水出露點(diǎn)不多見，沿結(jié)構(gòu)面滲出，多呈滴水，局部可見小股線流，流量小，初始洞內(nèi)除幾處出水點(diǎn)外，部分洞段潮濕，大部分洞段干燥。

輸水系統(tǒng)沿線基巖為凝灰?guī)r，后期侵入花崗斑巖及煌斑巖脈。圍巖完整性好，巖質(zhì)堅硬。根據(jù)工程建設(shè)期地下巖體灌前壓水試驗數(shù)據(jù)，圍巖滲透系數(shù)反演范圍為K=1.0×10-9～1.0×10-7m/s。根據(jù)地下巖體灌前壓水試驗數(shù)據(jù)，固結(jié)灌漿和灌漿帷幕滲透系數(shù)反演范圍為K=1.0×10-10～1.0×10-8m/s。

電站自1998年充水發(fā)電以來，電站運(yùn)行安全。運(yùn)行近20年，獲得了大量地下水、滲漏量等實(shí)測數(shù)據(jù)和資料。利用這些實(shí)測數(shù)據(jù)，反演工程區(qū)域三維滲流場，對地下水分布、滲控效應(yīng)、帷幕防滲效果評價、修復(fù)方案制訂等有重要參考意義。

3.2 三維滲流場優(yōu)化反演分析模型

三維滲流場反演計算區(qū)域選取上庫到下庫范圍，包含引水發(fā)電系統(tǒng)地下洞室群所有建筑物。X軸與引水鋼管軸線平行，指向下游為正；Y軸與引水鋼管軸線垂直，從1號鋼管指向6號鋼管為正；Z軸與大地坐標(biāo)重合，指向上為正。整體模型范圍：X方向從-615.61到472.14共1 087.75 m；Y方向從-266.70到347.20共613.90 m；Z方向從64.94到1 049.59共984.65 m。模型采用空間8節(jié)點(diǎn)等參單元，共劃分等參單元230 732個。計算模型如圖3和圖4所示。

圖3 三維滲流反演分析有限元計算模型Fig.3 FEM model of 3D seepage back analysis

圖4 三維滲流反演分析地下洞室有限元計算模型Fig.4 FEM model of the underground caverns and tunnels in 3D seepage back analysis

3.3 三維滲流場優(yōu)化反演分析結(jié)果

以地下水位測孔、平硐量水堰作為反演主要依據(jù)。通過100代例子反演，達(dá)到收斂標(biāo)準(zhǔn)。反演所得上下部圍巖和固結(jié)灌漿滲透系數(shù)如表1所示，反演結(jié)果均在鉆孔實(shí)測范圍之內(nèi)。從滲透系數(shù)量值大小上看，上部圍巖最大，下部圍巖次之，固結(jié)灌漿最小。結(jié)合地質(zhì)勘察資料，在斜井以上部位發(fā)育多條斷層，從宏觀上來看，上部圍巖滲透系數(shù)要比下部圍巖大。固結(jié)灌漿提高圍巖整體性的同時也降低了透水性，因而其滲透系數(shù)理論上要比圍巖小。因此，反演的參數(shù)大小關(guān)系合理。

表1 反演滲透系數(shù)表Table 1 Permeability coefficients in the back analysis

反演各水位測孔相對誤差均在5%以內(nèi)，反演結(jié)果合理。個別水位測孔誤差較大，主要原因是復(fù)雜的地質(zhì)條件，各種斷層、節(jié)理等的分布對滲流場都會產(chǎn)生影響，是導(dǎo)致誤差的一個重要原因。

工程區(qū)大范圍滲流等水頭線光滑，從上庫到下庫呈現(xiàn)逐漸減小的趨勢。在上下部圍巖交界面上，滲流等水頭線出現(xiàn)微小折線，反映了兩種圍巖滲透系數(shù)的差異。在輸水系統(tǒng)高壓管道附近，滲流等水頭線出現(xiàn)了上壅的趨勢，反映了斜井內(nèi)水外滲導(dǎo)致上體地下水位壅高的效應(yīng)。由于山體陡峭，上體地下水位和下庫水位相對都較高，廠房區(qū)位于自由水面之下，這與實(shí)際觀測中廠房頂拱滲水相吻合。在岔管和鋼襯接頭處，滲流等水頭線密集，說明該處滲流坡降大，滲透荷載消散較快。而在鋼管段和廠房段，滲流等水頭線稀疏，說明該處滲流坡降較小，這說明灌漿帷幕和排水孔幕很好地阻止了斜井的地下水滲漏到高壓岔管段和廠房段，其阻排水效果顯著。在廠區(qū)完善的排水廊道、排水孔和灌漿帷幕的排水、阻水作用下，整個地下廠房周圍滲透水頭很小。

總體來看，反演結(jié)果與實(shí)測資料相差不大，反演所得滲流場可以反映工程區(qū)大范圍的實(shí)際地下水分布和滲流狀態(tài)，此滲流反演結(jié)果可用于指導(dǎo)工程運(yùn)行、維護(hù)和檢修。

4 結(jié)語

基于工程區(qū)域長期實(shí)測水位、流量數(shù)據(jù)，采用粒子群智能優(yōu)化算法進(jìn)行滲流場優(yōu)化反演，顯著提高了復(fù)雜三維滲流場的反演分析效率和計算精度，應(yīng)用于實(shí)際工程的三維滲流場反演分析中，較好地解決了復(fù)雜三維滲流場反演難題。

（1）對于滲透性表征單元體存在且相對計算模型較小的巖土工程滲流問題，可按等效連續(xù)介質(zhì)模型進(jìn)行滲流分析。

圖5 運(yùn)行期輸水發(fā)電系統(tǒng)反演三維滲流場壓力水頭分布Fig.5 Distribution of water head along water transmission and power system in 3D seepage field back analysis

圖6 運(yùn)行期高壓管道區(qū)域反演三維滲流場壓力水頭分布Fig.6 Distribution of water head around the high pressure bifurcated tunnel during operation in 3D seepage field back analysis

（2）基于連續(xù)介質(zhì)的有限單元法，采用粒子群優(yōu)化算法進(jìn)行滲流場反演，對于邊界條件能相對準(zhǔn)確確定的問題，通過文中所述優(yōu)化反演分析方法，反演獲得的參數(shù)及滲流場分布是比較準(zhǔn)確和可靠的，反演精度基本能滿足工程需要。

（3）三維滲流場反演分析方法可有效應(yīng)用在巖土工程地下水分析中，但對實(shí)測數(shù)據(jù)準(zhǔn)確性、合理性要求較高。實(shí)測數(shù)據(jù)信息越準(zhǔn)確，反演結(jié)果越精確，同時實(shí)測數(shù)據(jù)信息越多，反演結(jié)果也越可信。

（4）反演分析方法除應(yīng)用于地下巖土工程滲流反演分析外，也可應(yīng)用于大壩、邊坡等工程滲流反演分析中。 ■

[1]杜延齡,許國安.滲流分析的有限單元法和電網(wǎng)絡(luò)法[M].北京：水利電力出版社,1992：1-2.

[2]張有天.巖石水力學(xué)與工程[M].北京：中國水利水電出版社,2005：8-25.

[3]朱岳明,張燎軍,吳愔.裂隙巖體滲透系數(shù)張量的反演分析[J].巖石力學(xué)與工程學(xué)報,1997,16(5)：461-470.

[4]劉先珊,周創(chuàng)兵,張立君.基于模擬退火的Gauss-Newton算法的神經(jīng)網(wǎng)絡(luò)在滲流反分析中的應(yīng)用[J].巖土力學(xué),2005,26(3)：404-414.

[5]盛金昌,速寶玉,詹美禮.三維天然滲流場反演分析及工程應(yīng)用[J].巖石力學(xué)與工程學(xué)報,2003,22(2)：203-207.

[6]謝紅強(qiáng),何江達(dá),符文熹,等.大崗山水電站廠壩區(qū)各向異性滲流場反演分析[J].四川大學(xué)學(xué)報（工程科學(xué)版）,2007,39(增)：132-135.

[7]張璇,李毅,羅琳,等.基于遺傳神經(jīng)網(wǎng)絡(luò)的壩區(qū)初始滲流場反分析[J].水電能源科學(xué),2012,30(12)：74-77.

[8]段斌,張林,何江達(dá),等.復(fù)雜裂隙巖體天然滲流場反演分析[J].水力發(fā)電學(xué)報,2012,31(3)：188-198.

[9]魏進(jìn)兵,鄧建輝,高春玉,等.三峽庫區(qū)泄灘滑坡非飽和滲流分析及滲透系數(shù)反演[J].巖土力學(xué),2008,29(8)：2262-2266.

[10]Duncan C Wyllie,Christopher W Mah.Rock Slope Engineering(4th edition)[M].London：Spon Press,2005：74-78.

[11]倪紹虎,何世海,汪小剛,等.裂隙巖體水力學(xué)特性研究[J].巖石力學(xué)與工程學(xué)報,2012,31(3),488-498.

[12]WEI Z Q，HUDSON J A.Permeability of jointed rock masses[C].Proceeding of ISRM Symposium：Rock Mechanics and Power Plants,Madrid：Romana，1988：613-625.

[13]WEI Z Q，EGGER P，DECOEUDRES F.Permeability predictions for jointed rock masses[J].International Journal of Rock Mechanics and Mining Sciences&Geomechanics Abstract，1995，32(3)：251-261.

[14]周志芳,王錦國.裂隙介質(zhì)水動力學(xué)[M].北京：水利電力出版社,2004：241-245.

[15]Kennedy J,Eberhart R C.Particle swarm optimization[C].Conf.on Neural NetworksⅣ,1995：1942-1948.

[16]Shi Y H,Eberhat R C.Parameter selection in particle swarm optimization[J].In Evolutionary ProgrammingⅦ：1998：591-600.

[17]Shi Y H,Eberhat R C.A modified particle swarm optimizer[C].Proceedings of the IEEE International Conference on Evolutionary Computation,1998：591-600.

[18]高瑋.基于粒子群優(yōu)化的巖土工程反分析研究[J].巖土力學(xué)，2006,27(5)：795-798.

[19]倪紹虎,肖明.基于圍巖松動圈的地下工程參數(shù)場位移反分析[J].巖石力學(xué)與工程學(xué)報,2009,28(7)：1439-1446.

[20]倪紹虎,肖明,王繼偉.改進(jìn)粒子群算法在地下工程反分析中的運(yùn)用[J].武漢大學(xué)學(xué)報（工學(xué)版）,2009,42(6)：326-330.

[21]常曉林,喻勝春,馬剛,等.基于粒子遷徙的粒子群優(yōu)化算法及其在巖土工程中的應(yīng)用[J].巖土力學(xué),2011,32(4)：1077-1082.

[22]倪紹虎,肖明,何世海,等.地下工程并行優(yōu)化反演分析及算例驗證[J].巖石力學(xué)與工程學(xué)報,2013,32(3)：501-511.