基于鹽巖腔體收斂過程的地表變形動態預測

童 凱，杜圣賢，賈 超，狄勝同，馬玉驍

(1.山東大學土建與水利學院，山東 濟南 250061；2.山東大學海洋地質與工程研究所，山東 青島 266237；3.山東省地質科學研究院，山東 濟南 250013)

0 引 言

鹽巖因其易溶于水，所以水溶開采鹽礦是一種十分便捷經濟的方案。又因具有極低的滲透率、良好的損傷自愈能力和塑性變形大等優良特性，被全球各國廣泛的運用在石油、天然氣等資源的儲備上。

在鹽巖力學特性研究方面，楊春和等[1]、YANG等[2]、WANG等[3]、姜德義等[4]、HOU等[5]學者開展了諸多研究，這對掌握鹽巖特性十分重要。由于鹽巖的流變特性，使得鹽腔在運營過程中會發生收斂，進而引發地面沉降，威脅到地表建筑物的安全[6]。由此可見研究引起地表沉降的機制對于控制由于鹽礦開采引發的次生地質災害十分必要。然而我國在地下礦產開采引起地表沉降的研究方面起步較晚，劉寶琛等[7]最先引進并發展了概率積分法，并在煤礦地表變形預測中進行了運用。任松等[8]又將用于煤礦開采沉降計算中的概率積分法引入到鹽巖的沉降預測中，并根據該理論發展了新概率積分三維預測模型。李銀平等[9]利用半無限域內球形空洞受力收縮位移解，利用疊加原理得到考慮內壓作用下的地表沉降彈性積分解析解。陳雨等[10-11]采用Gaussian曲線表示沉降分布，結合腔體收斂函數，建立了傳遞函數法預測地表變形的理論在此基礎上綜合考慮影響腔體收斂的多種時效因素(腔體水溶速率、體積收斂系數以及收斂傳遞率等)，根據SCHOBER等[12]提出的兩階段收斂模型，建立了一套較為完整的變形動態預測方法。

基于傳遞函數建立的動態預測模型是一套較完整的用于求解鹽礦開采誘發地表沉降的計算模型，但模型中收斂系數和傳遞率對于地表沉降的動態預測影響較大，其值估計的準確性直接決定了地表沉降結果的準確性。因此，本文基于傳遞函數方法和冪指函數模型，在考慮腔體收斂過程的情況下，對腔體收斂形式進行合理的簡化，建立起一套適用于預測鹽巖開采誘發地面變形的動態預測模型。

1 腔體收斂變形理論模型

1.1 傳遞函數模型

傳遞函數模型考慮地下腔體收斂與地表變形之間的聯系，認為腔體臨空邊界所受應力以水平方向為主，假設其受力過程中水平方向的變形遠大于垂直方向的變形，從而建立起地下腔體收斂與地表變形之間的函數[10]，見式(1)。

式中：S(x,y)為地表沉降分布函數；a為考慮收斂體積經上覆巖層傳遞折減后到達地面形成的沉降體積與溶腔體積之比；G(x,y)為聯系地下腔體收斂和地表變形的關系函數；q(z)腔體收斂函數；zt、zb分別為腔體頂部和底部埋深。

1.2 收斂函數公式推導

基于傳遞函數法建立的動態預測模型可以運用于求解地面沉降問題[10,12]，詳見式(2)～(4)。

(2)

F(ζ,η,t′)=

式中：S(d,t)為距離腔體中心d處在t時刻的沉降量；Sm為腔體中心處的最大沉降量；d為計算點到腔體中心的距離；r～2=rt×rb，rt、rb其中分別是溶腔頂部和底部的影響半徑；T為建腔時間；ζ為收斂系數，η為傳遞率。

該模型在進行計算時，由于ζ、η參數的取值決定了計算結果的精度，而該參數又不易獲得，為減小這一影響，本文在此基礎上對腔體收斂過程進行一定合理簡化，建立一套新的計算模型。

假設腔體建成后運行t時刻時，腔體體積滿足式(5)關系。

V(t)=ψV0(5)

式中：V0為腔體建腔時體積；V(t)為t時刻收斂后的腔體體積；ψ為體積收斂率。

假設腔體收斂時主要發生水平方向的變形，而忽略垂直方向的變形。以圓柱形、橢圓形和梨形腔體為例，建立單腔變形預計模型，見圖1。利用圖1中的模型參數，根據腔體收斂的假設，便可以建立腔體的收斂函數q(z)。

圖1 腔體收斂模型Fig.1 Convergence model of salt cavern

橢圓形腔體收斂前后的橢圓可描述為式(6)～式(7)。

(6)

式中：Δr為腔體收斂半徑；B為橢圓形短半軸；A為橢圓形長半軸；zm為橢圓形腔體中心埋深。

將橢圓形腔體體積計算公式帶入式(5)可得式(8)。對于整個橢圓形，收斂體積可表達為式(9)。聯立式(6)～(9)，可得式(10)。

(8)

同理可得圓柱形腔體和梨形腔體的收斂函數,分別表達為式(11)和式(12)。

q(z)=π(1-ψ)×

式中：Zm為梨形腔體上下腔分界線埋深；h1、h2分別為梨形腔體上下腔高度；Rp為梨形腔體上下腔分界線處半徑。

1.3 沉降函數推導

利用高斯分布函數描述地下收斂與地表變形的傳遞函數[10]，見式(13)。

(13)

式中：R為腔體開采在的地表影響半徑，R=zcotγ；d為沉降計算點到腔體中心在地表的水平投影距離；γ為收斂影響角，一般小于45°[12]。

定義式(14)，并聯立式(10)～(12)和式(14),可以求得三種不同形狀腔體的地表沉降預測公式。圓柱形腔體地表沉降預測表達式見式(15)；橢圓形腔體體地表沉降預測表達式見式(16)；梨形腔體體地表沉降預測表達式見式(17)。

(14)

1.4 地表沉降與水平變形關系

地表沉降和地表水平變形之間的關系可用式(18)表示[6]。

(18)

式中，n為水平變形影響系數，可由現場監測數據所得，一般簡化計算取n=1。

考慮腔體收斂過程，并由此可以得到三種不同腔形開采下地表的水平變形計算公式。定義式(19)，圓柱形腔體的地表的水平變形計算見式(20)；橢圓形腔體的地表的水平變形計算見式(21)；梨形腔體的地表的水平變形計算見式(22)。

(19)

Uo(d)=

當式(15)～(17)和式(20)～(22)中的ψ=0時，即此時儲氣庫完全收斂，地表變形達到最大值。故式(15)～(17)和式(20)～(22)可統一表達為式(23)。

S=S0(1-ψ)(23)

式中：S0為腔體完全收斂的地表變形值；S為當腔體收斂率為ψ時的地表變形值。

至此，針對不同形狀單腔開采考慮腔體收斂過程的地表沉降變形預測公式已經全部給出。

1.5 地表沉降動態預測模型

地下采空區的存在勢必會對地表的建筑物和居民造成安全隱患[13]，而鹽巖的蠕變特性又使得地表變形隨著儲庫的運行持續變化。在煤田開采條件下，地表移動過程可從6個月延續到數年，而在鹽礦開采條件下，移動持續時間甚至達到100年以上[14]。在儲庫運行初期還很難看出由于鹽礦開采對地表建筑物的影響，但隨著時間的推移，地表變形部分受壓區可能會在移動期間遭受拉伸，反之亦然。因此在儲庫建設過程中除了考慮穩定后的地表變形狀態，還須考慮地表變形隨時間的發展過程[15]。

地表變形預測公式中，體積收斂率ψ隨著儲庫運行時間的不同取值也不同，若知道不同時刻腔體的體積收斂率即可求得不同時間的地表變形分布。目前確定腔體收斂體積主要有以下三種方法：現場聲吶測腔、數值模擬和工程經驗預測。考慮到現場聲吶測試的花費較高，又具有一定的滯后性，而經驗預估又具有很大的不確定性，因此在條件有限的情況下采用數值模擬進行求解腔體收斂體積是一個較為可行的辦法。

2 地表變形算例分析

2.1 計算模型及參數說明

為了驗證本模型的預測效果，筆者在Matlab中編寫了相應的計算程序，用以求解地面變形分布，同時運用Ansys建立相應的有限元模型，在FLAC3D中進行計算求解，將二者的計算結果進行對比。

本次僅針對橢球體單腔開采引起的地面變形進行預測，根據楊春和等[16]的研究成果，腔體采用長短軸為7∶3的橢球體模型。腔體取長半軸為70 m，短半軸比為30 m，腔體中心埋深1 000 m，各地層覆存深度及有限元模型見圖2。數值模擬的瞬時穩態性計算采用Mohr-Coulomb彈塑性模型，蠕變計算采用FLAC3D自帶的Cpower蠕變模型。參考文獻[17]和文獻[18]取各地層巖性及蠕變參數，詳見表1和表2。

圖2 巖層分布及有限元模型Fig.2 Rock strata distribution and finite element model

表1 各巖層力學參數Table 1 Mechanical parameters of each rock layer

巖性彈性模量/GPa泊松比凝聚力/MPa內摩擦角/(°)抗拉強度/MPa泥巖21.780.2322.0241.5鹽巖5.150.3001.5351.0砂巖20.000.2500.5350.5

表2 鹽巖蠕變參數Table 2 Creep parameters of salt rock

模型計算中相關參數參考之前學者研究和經驗進行取值，沉降預測公式中的a取0.8，影響角取37°。數值模擬直接從腔體運營期開始計算，暫不考慮腔體運行壓力和建腔過程對于地表沉降的影響。假設腔體運行壓力為0，雖然這將導致計算結果偏大，但對于儲庫設計考慮是偏安全的。

2.2 數值模擬與預測模型結果對比

通過數值模擬可以求得不同運行時間的腔體體積收縮量，再帶入式(23)即可求得地表沉降變形的動態過程。

2.2.1 地表沉降動態曲線對比

根據數值模擬求出不同蠕變時間的腔體體積收縮量，再利用預測模型求取地表中心沉降值。以數值模擬的中心沉降值為標準，將預測模型和數值模擬的誤差統計在表3中。

表3 鹽腔地表中心沉降值計算表Table 3 Central surface subsidence of salt cavern

將數值模擬和預測模型的計算的不同蠕變時間的地表中心沉降值繪制在圖3中，可以更直觀地比較二者之間的關系。

從圖3和表3中可以看出，預測模型和數值模擬的發展趨勢基本相同。蠕變前10年二者相對誤差較大，由于該階段沉降值較小，即便絕對差值只有0.84 mm的第一年，相對誤差也達到了18.8%，隨著蠕變時間的增加，誤差逐漸減小，蠕變50年時，相對誤差僅有1.3%，整個計算周期內，最大沉降差值為1.4 mm，說明本預測模型具有較高的精度。

2.2.2 地表變形空間分布對比

為了進一步檢驗預測模型的地表變形在空間中的分布情況，計算儲氣庫運行50年時，以數值模擬的地表變形值為標準，將預測模型和數值模擬的誤差統計在表4中，將其與數值模擬的計算結果繪制在圖4中。

圖3 數值模擬和預測模型地表中心沉降對比Fig.3 Comparison of surface subsidence between numerical simulation and prediction model

表4 地表變形分布曲線對比Table 4 Comparison of surface deformation value

計算結果水平距離/m預測模型/mm數值模擬/mm相對誤差/%0-19.89-19.641.3100-19.54-19.142.1200-18.50-17.972.9300-16.93-16.403.2400-14.93-14.453.3地表沉降/mm500-12.70-12.342.9600-10.43-10.172.6700-8.26-8.735.4800-6.32-7.019.8900-4.66-5.6617.71 000-3.32-4.2622.100.000.20-1001.961.970.52003.713.904.93005.095.466.84005.996.457.1水平變形/mm5006.366.877.46006.276.676.07005.796.338.58005.065.7912.69004.205.0316.51 0003.324.1019.0

從圖4和表4中可以看出，數值模擬和預測模型計算的地表中心沉降最大值計算結果分別為19.64 mm和19.89 mm，相對誤差僅為1.3%。在整個計算區域內，當徑向距離小于600 m時，數值模擬的計算結果稍小于預測模型，但隨著距離的增大，數值模擬的結果會偏大一些。當距離達到1 000 m時，由于地表沉降值已經很小，即使只有0.94 mm的差值，相對誤差也達到了22.1%。

分析水平變形分布曲線和誤差統計可以得出，地表水平變形最大值發生在徑向距離為500 m的位置，二者計算的最大值分別為6.36 mm和6.87 mm，相對誤差7.4%，計算結果很接近。徑向距離小于700 m時，二者的相對誤差均小于10%，但隨著徑向距離的增大，很明顯數值模擬的計算結果較預測模型稍大，最大誤差達到了19%。

綜合地表沉降和水平變形分布曲線的對比結果，可以得出以下結論，本文所提出的預測模型計算的地表變形空間分布規律和數值模擬一致，但在距離腔體中心較遠時，預測結果存在一定誤差。但主要關心的變形較大的集中區域擬合程度較好，因此在一定條件下能滿足工程精度需求。

圖4 儲氣庫運行50年地表變形對比Fig.4 Comparison of surface deformation in operation for 50 years of gas reservoir

2.3 群腔地面變形動態預測

為了保證地表變形在可接受范圍內，盡可能多地開采鹽巖，充分利用鹽礦資源，使有限空間內的經濟效益達到最大化，實行群腔開采。運用本文提出的地表變形動態預測模型，再利用疊加原理進行群腔的地面沉降動態預測，假設各腔體沉降計算的參數相互獨立[12]，則地表任一點p(xi,yi)的變形值可由式(24)表達。

(24)

式中：Sp(xi,yi)(t)為點p在t時刻的地表沉降值；di為點p到i號腔體的徑向距離。

各地層分布情況及腔體形狀和圖2相同，開采形式如圖5所示，鄰腔夾角60°，腔距為2.0D[19]。基于疊加原理的假設可知，各腔體之間的變形互不影響，故以表3中的腔體收斂率作為該算例中各腔體的收斂動態過程，將地表沉降結果繪制在圖6中。

圖5 群腔布置平面圖Fig.5 Plan position of salt caverns field

圖6 群腔開采地面沉降曲線Fig.6 The subsidence curve of salt caverns

采用圖5所示鹽群布置時，地面沉降最大值發生在腔群中心地表位置，在50年時達到了132.5 mm，是相同單腔開采條件下的6.66倍。由此可見，群腔開采對地面沉降的影響很大，因此需要對鹽巖開采規模和腔群布置方式進行有效的控制。

3 結 論

1) 本文提出的預測模型所含未知參數較少，計算方便，在確定某工程的埋深及腔體形狀的前提下，可以通過預估或者參考類似工程儲庫運行年份的腔體收斂情況，求得地表變形的動態過程。根據不同的腔群布置形式可以很方便的求解地表沉降分布，從而進行優化設計。

2) 通過和數值模擬的計算結果對比，本模型在變形動態預測和空間分布上具有較高的準確性，地表中心最大沉降值的動態預測差值最大為1.4 mm。當徑向距離小于700 m時，地表變形的空間分布誤差均小于10%。腔群開采的沉降計算表明，多腔儲庫建設對地表沉降的影響顯著，所以在工程設計之初有必要針對腔群布置形式和建設規模進行進一步的研究和優化。

由于缺少實際監測資料，計算結果的準確性有待驗證，同時本模型未考慮運行壓力以及造腔過程的影響，這也必然導致計算結果與實際情況存在一定偏差。模型計算中的參數還需根據實際儲庫運行情況和腔體收斂結果進行進一步校核，以提高預測精度。