線性工程沿線地下水超采引發(fā)不均勻地面沉降分析

賈 超，邊 超，丁朋朋，楊 霄

（山東大學海洋研究院，山東青島266232）

0 引 言

地面沉降是一種在人為因素和自然因素作用下地面高程緩慢降低的地質(zhì)災(zāi)害現(xiàn)象［1，2］。地下水的超量開采、石油和水溶性氣體的開采、自然固結(jié)和建筑荷載等是地面沉降常見的誘發(fā)機制［3］。其中，地下水超采而引發(fā)的地面沉降約占80%［4］。目前我國由于地下水超采引起的地面沉降主要出現(xiàn)在天津、北京、上海、江蘇、山東和山西等中東部地區(qū)，累計地面沉降量超過200 mm的國土面積大約有7.9 萬km2［5］。

地面沉降具有初期不易察覺、地域性明顯等特點，對鐵路、公路、管道等工程危害巨大。區(qū)域均勻沉降會導致路基的絕對高程降低，一旦路基的絕對高程降低，鐵路、公路等工程抵御極端惡劣災(zāi)害（如洪水）的能力就會降低。相較于均勻沉降，不均勻沉降的危害更為嚴重，它會破壞軌道線路、路基等［6］，具有高復(fù)雜性、局部化與快速發(fā)展等特征［7］。

因此，在地面沉降多發(fā)區(qū)的工程地質(zhì)選線中，必須充分考慮地面沉降的影響［8］。許多學者運用數(shù)值模擬、數(shù)學評價、監(jiān)測技術(shù)等手段對線性工程沿線地面沉降的分布特征及發(fā)展歷史開展了相關(guān)研究［9-12］。例如，李國和［9］采用水土耦合模型結(jié)合反分析方法對雅萬高鐵沿線地面沉降進行概要分析，并提出適宜的工程防治措施；祁彪［10］采用統(tǒng)計分析方法對某鐵路沿線地面沉降的幅度、速率及坡度進行了計算和預(yù)測，認為不均勻沉降對高鐵工程有較大影響并且鐵路運營會加劇不均勻沉降；謝海瀾等［11］分析了京津冀地面沉降年均沉降速率和高鐵線路坡度變化，將坡度變化作為地面沉降承載狀態(tài)評價的指標，建立地面沉降與高鐵安全運行之間的關(guān)系；Luo 等［12］探索了多時相合成孔徑雷達干涉測量技術(shù)（MT-InSAR）對京津城際鐵路地面沉降的監(jiān)測能力，并與高時空密度水準測量技術(shù)相結(jié)合，得到了更高精度的沉降數(shù)據(jù)，可以為重大線性工程地面沉降提供更精確的評估。

為了深入研究不均勻地面沉降對線性工程形變的作用程度，許多研究人員運用曲率分析的方法對其形變特征進行了定量分析［13，14］。例如，戴文彬［13］根據(jù)西安市地面沉降量等值線圖，按照地鐵走向繪制了沿線地面沉降曲線，得到了曲線的最小曲率半徑及累積地面沉降差和沉降曲線最小曲率半徑的關(guān)系，分析了西安市不均勻地面沉降的特征；Wu 等［14］分析了地面沉降下管道在不同位置的應(yīng)力變化規(guī)律，并采用最小曲率半徑和最小變形表征管道的變形程度。

本文根據(jù)地下水三維滲流和比奧固結(jié)理論，應(yīng)用多物理場有限元數(shù)值模擬軟件建立了地面沉降三維流固耦合數(shù)值模型，將魯南高鐵沿線區(qū)域作為研究區(qū)，定量評估研究區(qū)不同地下水開采方案下不均勻地面沉降發(fā)展趨勢，研究不均勻地面沉降對高鐵線路的影響。研究成果可為線性工程的選線、建設(shè)施工和后期運營提供技術(shù)支持，同時可為地下水超采誘發(fā)不均勻地面沉降的定量分析提供參考和借鑒。

1 理論與方法

1.1 地面沉降流固耦合數(shù)學模型

本文以比奧固結(jié)理論和達西定律為基礎(chǔ)［15，16］，假設(shè)土體飽和，微壓縮，且土骨架變形滿足線彈性理論，結(jié)合達西定律將三維應(yīng)力應(yīng)變關(guān)系結(jié)合三維地下水流方程，建立三維流固全耦合數(shù)學模型，其中，數(shù)值模擬的控制微分方程包括：運動方程、平衡方程、本構(gòu)方程［17］及邊界條件方程［18-20］。

1.1.1 運動方程：

式中：qi為通過單位面積的流量，m/s；kija為土中水的表觀流動性系數(shù)張量，m2/Pa·s，其中a為飽和度s的值，P為流體壓力，Pa；ρw為流體密度，kg/m3；gk為重力加速度矢量，m/s2。

1.1.2 平衡方程

對于小變形：

式中：ζ為單位體積的孔隙材料中流體體積的變化；qv為體積流源強度，L/s。

動量平衡的形式如下：

式中：σij是總應(yīng)力張量（Pa，以壓為正）；ρ=（1-n）ρs+n ρw是體積密度，ρs和ρw分別為固相和液相的密度。值得注意，（1-n）ρs對應(yīng)的是干密度ρd；ui為流速。

1.1.3 本構(gòu)方程

流體流動性的改變與孔隙水壓力p、飽和度s、體積應(yīng)變ε、溫度T等因素有關(guān)，孔隙流體的響應(yīng)方程可有如下表達：

式中：M是Biot 模量（Pa，M=Kwn，Kw為流體體積模量，n為孔隙率）；α為Biot系數(shù)；β為考慮流體和顆粒熱膨脹系數(shù)，℃-1。

在完全飽和狀態(tài)下，即s=1，kija=kij并且流體可以維持極限拉力的溫度為Tf。對此情況，響應(yīng)方程可簡化為：

結(jié)合式（1）、（2）和（5）得到滲流分析的微分控制方程為：

1.1.4 邊界條件方程

在計算中涉及到許多邊界條件，分別是：①地下水定流量邊界；②土體應(yīng)力邊界；③透水邊界。

對于定流量邊界［18］（設(shè)為S1或Γ1），邊界條件為：

或

式中：n為邊界S1或Γ1的外法線方向；q1和q2均為已知函數(shù)，表示S1上單位面積和Γ1上單位寬度的側(cè)向補給量。

設(shè)S2為應(yīng)力場的應(yīng)力邊界［19］，土體的應(yīng)力邊界條件可表示為：

式中：Fi為邊界應(yīng)力張量；nj為邊界表面的法外線矢量。

設(shè)S3為透水邊界［20］，透水邊界采用如下形式給出：

式中：qn是邊界外法線方向單位長度的流量，m2/s；p是邊界面處的孔隙水壓力；pe是滲流出口處的孔隙水壓力。

1.2 曲率分析

曲率是一種體現(xiàn)區(qū)域不均勻沉降的物理量［13］，根據(jù)對沿線沉降曲線圖的分析，求出曲線上某個點的切線方向角對弧長的轉(zhuǎn)動律。通過曲率的大小和地面一定范圍內(nèi)的沉降差的關(guān)系來判斷該區(qū)域不均勻沉降程度［21］。

若一條曲線L的函數(shù)為p=p（s），則該曲線的曲率k為：

對平面曲線p=p[x(s)，y(s)]，則其一階、二階導數(shù)分別是：

根據(jù)以上公式，計算得到地面沉降曲線的曲率k及曲率半徑R，其中曲率半徑R=1。

2 應(yīng)用算例

2.1 研究區(qū)概況

魯南高速鐵路菏澤至曲阜段位于山東省西南部，東起濟寧市管轄的曲阜市，西至菏澤市管轄的牡丹區(qū)，途經(jīng)濟寧市管轄的兗州區(qū)、任城區(qū)、汶上縣和嘉祥縣，菏澤市管轄的巨野縣、鄆城縣、定陶縣。研究區(qū)屬黃河沖積平原區(qū)，山前傾斜平原，地形平坦，地勢開闊。沿線以新生界第四系地層為主，厚度變化較大，從東往西厚度逐漸增加，濟寧附近約100～200 m，菏澤附近約300～400 m。在曲阜南部局部為剝蝕殘丘區(qū)，下伏基巖多為奧陶系、寒武系灰?guī)r等可溶巖地層。

以魯南高鐵施工圖階段明確的貫通方案為主線，拓寬線位兩側(cè)各10 km 為研究區(qū)范圍，面積約3 800 km2，如圖1所示。據(jù)現(xiàn)場調(diào)查，高鐵沿線多個城鎮(zhèn)（如菏澤及濟寧等）大量抽取地下水，形成多個以城鎮(zhèn)為中心的地下水降落漏斗，地下水位降低引起地面沉降災(zāi)害日益嚴重。

2.2 三維流固耦合模型

2.2.1 模型的建立

根據(jù)已有的鉆孔資料，確立地層埋深范圍為0～-500 m，選取高鐵線路兩側(cè)10 km 區(qū)域為模型的平面范圍，建立地面沉降三維數(shù)值模型，對地下水開采與地面沉降進行模擬研究。按地層巖性及土體參數(shù)將模型垂向上剖分為六層。第一層為表土層；第二層和第四層為隔水層；第三層和第五層為含水層，是主要采水層。具體地層劃分如圖2所示。

結(jié)合高鐵沿線地區(qū)水文地質(zhì)、工程地質(zhì)資料以及鉆孔土樣實驗結(jié)果，確定了各地層的巖土體物理力學參數(shù)，如表1所示。

表1 各地層土體的物理及力學參數(shù)Tab.1 Physical and mechanical parameters of various layers of soil

采用自由剖分四面體網(wǎng)格對模型進行離散化，剖分后該研究模型由522 836 個域單元，139 506 邊界單元和6 241 邊單元以及351 864 個自由度組成。模型三維結(jié)構(gòu)如圖3所示。

進行具體分析計算前，需要確定高鐵線路沿線區(qū)域各地下水開采井的開采量，在解析計算和數(shù)值實驗中，開采量一般采用流量或流速表示，為簡化計算，方便進一步進行三維數(shù)值模擬，將開采井的流量轉(zhuǎn)換為質(zhì)量通量表示各開采井的開采量［20］。研究區(qū)內(nèi)地下水開采井位置分布如圖4所示。

2.2.2 初始條件和邊界條件

模型邊界條件根據(jù)實測水文資料確定。以2014年為初始時刻，假設(shè)該時刻所有點的垂向位移為零，即為沉降初始時刻。側(cè)面為對稱約束，即水平方向位移為零；底面為固定邊界，即水平、豎直位移為零；上表面為自由邊界。上表面為對稱面，即流量為零；底面為不透水面，即流量為零；側(cè)面邊界按水頭變化值計算；井點均為流量邊界，流量大小按實際流量施加，本研究中按轉(zhuǎn)換后的質(zhì)量通量施加。

2.2.3 模型有效性及準確性驗證

從地面沉降時空分布趨勢和多年動態(tài)沉降值兩個方面對流固耦合三維數(shù)值模型進行校核和驗證。將數(shù)值模擬沉降結(jié)果與InSAR 逐年監(jiān)測數(shù)據(jù)對比分析，驗證三維流固耦合數(shù)值模擬沉降時空分布趨勢的準確性和有效性。

（1）地面沉降分布趨勢校準：基于研究區(qū)2014.01-2018.12期間InSAR 地面沉降速率實測數(shù)據(jù)（如圖5所示），將其與數(shù)值模擬計算結(jié)果對比分析，如圖6所示，兩者對比確定三維流固耦合數(shù)值模型地面沉降趨勢準確性和有效性。

從數(shù)值模擬結(jié)果和InSAR 數(shù)據(jù)對比發(fā)現(xiàn)，地面沉降數(shù)值模擬結(jié)果和InSAR 監(jiān)測結(jié)果總體上趨勢處于一致，但是在郭屯煤礦、趙樓煤礦、古城煤礦和興隆莊煤礦區(qū)域存在差別。調(diào)查發(fā)現(xiàn)這些區(qū)域的地下水開采井基本是農(nóng)田灌溉井和村鎮(zhèn)飲用井，年開采量僅為0.39～5.0 萬m3，地下水開采引起的最大沉降值約為76 mm，占監(jiān)測結(jié)果的33%。因此，煤礦開采等其他因素是這些區(qū)域發(fā)生地面沉降的主要影響因素。數(shù)值模擬結(jié)果在校核期內(nèi)的地面沉降模擬效果較為合理，能夠反映研究區(qū)由于地下水開采引發(fā)的實際地面沉降趨勢。

（2）地面沉降監(jiān)測點校準：高鐵沿線8個典型監(jiān)測點二等水準測量結(jié)果定量化校準研究區(qū)地面沉降三維數(shù)值模型，空間位置分布如圖1所示。分別將沉降監(jiān)測點2014-2018年的實際監(jiān)測數(shù)據(jù)與數(shù)值模擬結(jié)果進行對比分析，運用線性擬合的方法，驗證模型的有效性和準確度，如圖7所示。結(jié)果表明，2014-2018年監(jiān)測點處沉降模擬結(jié)果與監(jiān)測結(jié)果擬合的可決系數(shù)R均大于0.9，擬合效果較好，說明所建立的數(shù)值模型可用于研究區(qū)地面沉降的模擬分析。

2.3 不均勻地面沉降分析預(yù)測

采用已建立的地面沉降三維數(shù)值分析模型對不均勻地面沉降動態(tài)變化進行分析預(yù)測，重點分析研究了現(xiàn)狀開采和減采條件下未來10年高鐵沿線區(qū)域不均勻地面沉降量的動態(tài)特征。

本文采用曲率的大小來表征高鐵沿線不均勻沉降的程度。曲線的曲率是針對曲線上某個點的切線方向角對弧長的轉(zhuǎn)動率，曲率越大，在曲線上該點的彎曲程度越大，體現(xiàn)在高鐵沿線的沉降曲線上則說明地面沉降的不均勻性越嚴重。

2.3.1 地下水現(xiàn)狀開采

以2018年為分析預(yù)測初始時刻，2014年的沉降量作為初始條件，地下水開采量與2018年保持相同，保持模型位移邊界條件不變，對未來10年（2018-2028年）地面沉降趨勢及規(guī)律進行了模擬計算。計算得到研究區(qū)2028年地面沉降結(jié)果及高鐵沿線的沉降曲線剖面圖分別如圖8、圖9所示。

根據(jù)高鐵沿線地下水現(xiàn)狀開采條件下沉降預(yù)測結(jié)果及曲率計算公式得到高鐵沿線的沉降曲率，繪制模擬至2018 和2028的沉降曲率分布圖，如圖10所示。

經(jīng)分析發(fā)現(xiàn)，沉降預(yù)測剖面圖與沉降曲率分布圖在變化趨勢中呈現(xiàn)極高的相似度，例如DK305-DK320 段和DK360-DK420 段都是地面沉降較為復(fù)雜的區(qū)域，在曲率上表現(xiàn)的也是波動最厲害的，可以證明曲率可以較好地表征不均勻沉降程度。

2018年及2028年預(yù)測的最大沉降曲率的樁號點均為DK415，2018年的最大沉降曲率為7.221 95×10-5m-1，累積沉降量為102.021 mm；2028年的最大沉降曲率為7.647 54×10-5m-1，累積沉降量為273.142 mm。

2.3.2 地下水減采10%

為了有效地控制地面沉降對高鐵沿線的影響，對研究區(qū)在現(xiàn)狀開采條件的基礎(chǔ)上減采10%的地下水量，保持模型位移邊界條件不變，將2018年的地下水開采量減少10%作為預(yù)測期內(nèi)每一年的地下水開采量，模擬結(jié)果如圖11所示。根據(jù)沿線2018年和2028年累計地面沉降量繪制沿線地面沉降預(yù)測剖面圖，如圖12所示。

沿線的累積沉降量相較現(xiàn)狀開采條件下，有明顯地減少，10年的沉降速率也有減緩，表明減少地下水開采量可以會有效降低地面沉降的速率。

根據(jù)高鐵沿線地下水減采10%地下水條件下沉降預(yù)測結(jié)果及曲率計算公式得到高鐵沿線的沉降曲率，繪制模擬至2018年和2028年的沉降曲率分布圖，如圖13所示。沉降預(yù)測剖面圖與沉降曲率分布圖在變化趨勢中與現(xiàn)狀開采基本一致，且2018年及2028年預(yù)測的最大沉降曲率的樁號點均為DK415。2028年的最大沉降曲率為6.877 16×10-5m-1，累積沉降量為246.183 mm，與現(xiàn)狀開采條件下相對比可以看出限制開采地下水可以有效減緩地面沉降速率并且緩解不均勻沉降。

沉降預(yù)測剖面圖與沉降曲率分布圖在變化趨勢中與現(xiàn)狀開采基本一致，且2018年及2028年預(yù)測的最大沉降曲率的樁號點均為DK415。2028年的最大沉降曲率為6.877 16×10-5m-1，累積沉降量為246.183 mm，與現(xiàn)狀開采條件下相對比可以看出限制開采地下水可以有效減緩地面沉降速率并且緩解不均勻沉降。

2.3.3 綜合分析

兩種開采地下水方案的結(jié)果分析表明，在2018-2028年時間段內(nèi)，地面沉降量均逐年增加；隨著土層的固結(jié)壓密，年沉降速率逐年降低；減采條件下各年的沉降速率均低于現(xiàn)狀開采條件。

現(xiàn)狀開采條件下和減采10%條件下，鐵路沿線DK385處的2018-2028年曲率變化，如圖14所示。同一開采條件下，沉降點曲率經(jīng)歷了快速增長、緩慢增長和穩(wěn)定階段等3 個階段?，F(xiàn)狀開采和減采10%在2028年的曲率分別為1.129×10-5m-1和1.014×10-5m-1。

現(xiàn)狀開采和減采均會導致曲率的上升，但現(xiàn)狀開采下各個階段的曲率增長的速率比減采10%條件下更快，所以為了高鐵的安全及長期的運行，若不采取措施減少地下水的開采，則會導致不均勻地面沉降越發(fā)嚴重，曲率逐年增大，影響高鐵的安全運行。

綜上，地下水開采造成地下水位降低，導致地面沉降的變化是明顯的。減少開采量可以作為控制高鐵沿線區(qū)域不均勻地面沉降的措施之一，有效降低地面沉降量和沉降曲率，保證高鐵的安全運營。

3 結(jié) 論

本文以魯南高鐵菏澤至曲阜段為研究區(qū)，建立了地面沉降流固耦合三維數(shù)值模型，對研究區(qū)地下水開采誘發(fā)地面沉降進行模擬預(yù)測，并分析了高鐵沿線不均勻地面沉降特征及演化趨勢，結(jié)論如下。

（1）基于地下水滲流理論和土體比奧固結(jié)理論建立了魯南高鐵沿線地面沉降三維流固耦合數(shù)值模型，該模型可較為精確的考慮地下水開采與地層變形相互耦合，對抽水導致土體變形機理進行了較為真實的描述。

（2）根據(jù)現(xiàn)有的開采資料，通過對比未來10年地下水現(xiàn)狀開采及減少10%開采來預(yù)測不同工況下沉降量和速率，發(fā)現(xiàn)減少地下水開采可以有效減緩地面沉降的速率。因此，在鐵路沿線一定范圍內(nèi)設(shè)立控采區(qū)是十分必要的。

（3）本文針對不同工況的高鐵沿線進行了地面沉降曲線的曲率分析?，F(xiàn)狀開采條件下曲率雖然符合規(guī)范要求，但若不采取控采措施，沉降曲率將逐漸增加，將會危害鐵路安全。因此，減少地下水開采可以有效降低地面沉降和曲率，可以保證線路運行的安全性。 □