王 琰,呂 航,谷復光*
(1.吉林建筑大學測繪與勘查工程學院,吉林 長春 130118;2.吉林大學新能源與環境學院,吉林 長春 130021)
隨著工業發展和城市化進程的逐步推進,地面沉降在全世界范圍內普遍發生,引發了建筑物開裂破壞、城市內澇、防潮抗洪能力降低等各種災害。經過眾多學者多年來的監測、分析與模擬研究,認為地下水開采是引發大面積地面沉降的主要原因,特別是在第四紀松散沉積層分布的地區,由地下水水位急劇、廣泛下降引起的地面沉降問題日趨嚴重。如何基于描述地下水流場的變化,確立準確、高效的地面沉降計算和預測方法,對于控制地面沉降的發展至關重要。
國內外許多學者基于各地區水文地質條件以及地下水水位和地面沉降量的動態監測結果,提出了一些地下水開采引起地面沉降的預測方法,主要包括數值計算法和數理統計法兩大類。數值計算法主要有兩步計算模型(如Gambolati等最早提出兩步計算模型)、部分耦合模型(如Fokker等在兩步計算模型的基礎上,結合地下水流方程和Terzaghi固結理論,建立了模擬多層含水層系統的地下水流和總沉降量的數學模型)、考慮三維滲流與一維固結模型(如薛禹群等建立了三維變系數滲流模型和垂向一維沉降模型,分析了長江三角洲南部區域的地面沉降問題)和三維滲流與三維固結完全耦合模型(如Luo等、徐成華等在黏彈塑性比奧固結理論的基礎上建立的地下水開采與地面沉降全耦合模型),該方法相對來講具有計算復雜、占用計算空間大、耗時多、許多參數難以確定等缺陷,采用數值計算方法進行地面沉降量的計算和預測時需花費大量的時間和人力進行模型建立和參數識別,計算速度相對較慢。相對于數值計算法,數理統計法在統計資料完整的情況下,擁有計算速度快、應用性較強等優勢。
本文在分析地下水滲流與地面沉降間相關原理的基礎上,采用數理統計方法,依據工程實例分析了地面沉降量與地下水水位變化之間的相關關系,建立了兩者之間的多元線性回歸模型,用來預測由地下水開采引起的地面沉降量的發展變化趨勢,并通過對預測模型的檢驗,證明了該方法在地面沉降預測方面具有一定的可靠性和實用性。
天然土體是由礦物顆粒骨架、孔隙水和氣體構成的固、氣、液三相體系。根據有效應力原理,在含水介質骨架壓縮的初始時段,土體中某一截面的上部荷載由有效應力(作用在礦物顆粒骨架上的應力)和孔隙水壓力共同承擔;隨著地下水的持續開采,土體中地下水水位逐漸下降,孔隙水壓力逐漸降低,作用在顆粒骨架上的有效應力隨之增加;由于有效應力逐漸增加,土體骨架發生變化,從而導致土體體積改變,土體產生壓縮變形。宏觀上,各層土體壓縮變形積累在地面即表現為地面沉降。
依據上述理論,分析并建立地下水開采過程中地下水水位變化值與土體變形量之間的關系,即可得到含水層的土體變形方程。
根據太沙基有效應力原理,考慮土體變形過程中側向受限,土體變形主要產生在垂直方向上,且一段時間內土體承受的總應力不隨時間變化等假設條件,有如下計算公式:
σ
′=σ
-u
(1)
式中:σ
′為含水層土體承受的有效應力(MPa);σ
為含水層土體承受的總應力(MPa);μ為含水層土體承受的孔隙水壓力(MPa),
Δu=ρ
g
ΔH[
其中ρ
為水的密度(
kg/m)
,g
為重力加速度(
m/s);
ΔH
為孔隙水頭變化值(
m)]
。無壓含水層系統中,含水層土體承受的有效應力的變化可表示如下:

(2)

土體中某一截面上,含水層土體承受的有效應力的變化可表示如下:
Δσ
′=
-ρ
g
ΔH
(3)
由于土體有效應力的變化會導致土體骨架發生壓縮,因此土體壓縮系數α
可表示為
(4)
式中:V
為所取單元體積(
m);σ
′為含水層土體承受的有效應力(MPa)。

(5)
式中:b
為可壓縮含水層的初始厚度(
m)
;b
為含水層厚度(
m)
。根據一維太沙基有效應力原理,將公式(5)代入公式(2)、(3)中,可得:

(6)

(7)
式中:μ
為含水層土體的骨架儲水率(
1/m)
,它是土體物理參數的函數。因此,公式(6)、(7)可表示如下:
Δb=
-ΔH(
1-n
+n
)μ
b
(8)
Δb=
-ΔHμ
b
(9)
式中:
Δb
為含水層彈性壓縮量與非彈壓縮量(
m)
,正為壓縮,負為擴張。公式(8)、(9)即為計算潛水含水層(潛水含水層或因抽水含水層由飽和狀態轉為非飽和狀態)土體垂向變形與承壓含水層土體垂向變形的基本公式。對應土體變形的不同階段,含水層土體骨架儲水率的變化分別如下:

(10)

根據公式(8)、(9),由地下水水位下降引起的含水層土體壓縮量的計算模型如下18:
潛水含水層的彈性變形量為
Δb=
-ΔH(
1-n
+n
)μ
b
=
-ΔHμ
(11)
潛水含水層的非彈性變形量為
Δb
=
-ΔH(
1-n
+n
)μ
b
=
-ΔHμ
(12)
承壓含水層的彈性變形量為
Δb=
-ΔHμ
b
=
-ΔHμ
(13)
承壓含水層的非彈性變形量為
Δb
=
-ΔHμ
b
=
-ΔHμ
(14)
式中:μ
為含水層土體骨架成分的彈性儲水因子(
無量綱)
;μ
為含水層土體骨架成分的非彈性儲水因子(
無量綱)
;其他符號意義同上。實際地下水開發利用中,上述建立的含水層土體壓縮量計算模型將被推廣至各黏性土弱透水層,含水層及包含在其中的弱含水層被視為一個整體,整個土體骨架都可發生壓縮,參數μ
、μ
是在整個土層規模上的等效參數。本研究選擇的地面沉降多層位監測場地位于浙江省東部、長江三角洲東南角,平原區第四系厚度為50~110 m,最厚達120余米。監測目的地層由下至上分別為第四系上更新統下組礫石和中細砂組成的第Ⅰ承壓含水層、上更新統上組下部細砂和含礫砂組成的第Ⅰ承壓含水層、上更新統上組中部亞黏土組成的第Ⅱ黏性土弱含水層、全新統下組亞砂土構成的淺部孔隙承壓含水層、全新統中組淤泥或淤泥質亞黏土構成的第Ⅰ黏性土弱含水層和全新統上組亞砂土構成的潛水含水層。場地所在區域自1982年開始少量開發地下水,開采量逐年增加,由于大規模開采地下水引起了較嚴重的地面沉降問題,為了保護環境,自2008年開始當地對地下水禁采,目前只預留少部分地下水井應急使用。
根據地面沉降機理,地面沉降與地下水水位的變化存在密切的聯系,故本文對研究場地內地面沉降及地下水水位監測點1985—2008年各含水層壓縮沉降量與相應的地下水水位變幅數據進行了整理與分析,詳見圖1至圖6。

圖1 潛水含水層沉降逐年變化量與地下水水位逐年 變化量對比Fig.1 Comparison of annual variation of land subsidence and groundwater level of phreatic aquifer

圖2 第Ⅰ黏性土弱含水層沉降逐年變化量與地下水 水位逐年變化量對比Fig.2 Comparison of annual variation of land subsidence and groundwater level of the first aquitard

圖3 淺部孔隙承壓含水層沉降逐年變化量與地下水 水位逐年變化量對比Fig.3 Comparison of annual variation of land subsidence and groundwater level of shallow confined aquifer

圖4 第Ⅱ黏性土弱含水層沉降逐年變化量與地下水 水位逐年變化量對比Fig.4 Comparison of annual variation of land subsidence and groundwater level of the second aquitard

圖5 第Ⅰ1承壓含水層沉降逐年變化量與地下水 水位逐年變化量對比Fig.5 Comparison of annual variation of land subsidence and groundwater level of the Ⅰ1 confined aquifer

圖6 第Ⅰ2承壓含水層沉降逐年變化量與地下水水位 逐年變化量對比Fig.6 Comparison of annual variation of land subsidence and groundwater level of the Ⅰ2 confined aquifer
從研究場地內各含水層沉降逐年變化量與地下水水位逐年變化量的對比結果來看,地面沉降變化量與地下水水位標高變化量數值在各監測年份存在著明顯的對應關系,且后續的地面沉降與地下水水位之間的回歸分析在數理上也充分證明了地面沉降量與地下水水位變化之間存在密切的相關性。
由前述理論分析和地下水水位監測點數據分析可知,地面沉降量與地下水水位的變化之間存在線性關系,可以用線性回歸方程來描述。本文基于數理統計方法,建立時間序列上地面沉降量與各含水層地下水水位變幅之間的多元線性回歸模型,以此為基礎,可以實現地下水開采引起的地面沉降量預測。
通過對計算點多年統計資料的分析,本文建立了研究場地內地面沉降量與地下水水位變化之間的多元線性回歸模型:
Δb=α
+α
ΔH
+α
ΔH
+α
ΔH
+α
ΔH
+α
ΔH
+α
ΔH
(15)
式中:Δb
為月地面沉降量(mm);ΔH
為潛水含水層月地下水水位變幅(m);ΔH
為潛水含水層與淺部孔隙承壓含水層間黏性土弱含水層月地下水水位變幅(m); ΔH
為淺部孔隙承壓含水層月地下水水位變幅(m); ΔH
為淺部孔隙承壓含水層與第Ⅰ承壓含水層間黏性土弱含水層月地下水水位變幅(m);ΔH
為第Ⅰ承壓含水層月地下水水位變幅(m);ΔH
為第Ⅰ承壓含水層月地下水水位變幅(m);α
為常數項,與各含水層和黏性土弱含水層的監測控制程度和監測精度以及非地下水開采引起的地面沉降有關;α
、α
、α
、α
、α
、α
為系數項,與各含水層和黏性土弱含水層的土體骨架釋水系數有關。根據已有266組因變量(月地面沉降量)與自變量(各含水層地下水水位變幅)的觀測統計數據,經計算得到了研究場地內月地面沉降量與各含水層地下水水位變幅之間的多元線性回歸方程:
Δb=
-0.87+1.22ΔH
+0.46ΔH
+2.61ΔH
+0.32ΔH
+0.63ΔH
+ 0.83ΔH
(
16)
本文所建立的多元線性回歸模型的擬合度檢驗結果,詳見表1。

表1 多元線性回歸模型擬合度檢驗結果Table 1 List of goodness of fit test of the multiple linearregression model
由表1可知,建立的多元線性回歸模型的復相關系數R
為0.917,表明模型自變量與因變量之間的密切程度很高;模型的決定系數(也稱判定系數)R
值為0.841,為避免隨自變量個數增加其可靠性降低,經調整后R
值為0.837(統計學上一般R
值大于0.8時說明其擬合程度較高),表明利用該多元線性回歸方程來描述和說明因變量的變化所解釋的程度較高;而模型R
的估計標準誤差值為2.682 26,其值較小,表明R
的估值具有一定的可靠性。本文所建立的多元線性回歸模型的方差分析結果,見表2。

表2 多元線性回歸模型方差分析結果Table 2 List of ANOVA of the multiple linear regressionmodel
由表2可知,建立的多元回歸模型的顯著性指標Sig.值為0,F
分布統計量值為195.534,說明該模型具有高度顯著的統計意義。本文所建立的多元線性回歸模型的共線性檢驗結果,見表3。

表3 多元線性回歸模型共線性檢驗結果Table 3 List of collinearity tests of the multiple linearregression model
由表3可知,建立的多元線性回歸模型的特征值均不為0(若多個維度的特征值等于0,則可能存在共線性問題),條件指數值均小于30(若某個維度的條件指數大于30,則存在共線性問題),說明該模型各自變量間存在共線性問題的可能性很小。
本文所建立的多元線性回歸模型的殘差分析結果,見表4。

表4 多元線性回歸模型殘差分析結果Table 4 List of residual statistics of the multiple linearregression model
由表4可知,建立的多元線性回歸模型的殘差統計量數據中無離群值,且數據的標準殘差也較小,可以認為該模型是健康的。
通過對本文所建立的多元線性回歸模型的檢驗,繪制了標準化殘差的直方圖和標準化殘差的正態P
-P
圖,見圖7和圖8。
圖7 標準化殘差直方圖Fig.7 Histogram of standardized residuals

圖8 標準化殘差正態P-P圖Fig.8 Normalized residual normal P-P diagram
由圖7和圖8可見,建立的多元線性回歸模型的標準化殘差分布具有正態分布的趨勢,說明模型是恰當、可用的。
綜上所述,本文所建立的多元線性回歸模型正確、可信,可以用來預測某處地面沉降總量月變化量的發展趨勢。
(1) 由地下水開采引起的地面沉降,地面沉降量與地下水水位變化之間存在時間上的對應,通過對地面沉降機理的研究證實兩者之間確實存在線性相關關系。本文應用數理統計方法,通過工程實例對地面沉降與地下水水位進行了回歸分析及預測,模型檢驗結果表明該預測模型是可用的。
(2) 在未來地下水開發利用規劃情況下,場地所在地區可應用該方法開展由地下水開采引起的地面沉降量預測,其計算速度快、實用性較高。
(3) 隨著全國范圍內地下水與地面沉降監測網的逐年完善,大量精度較高的監測數據將被獲取與積累。本文的研究方法可在第四紀松散沉積層較厚的沖、洪積平原地區進行廣泛的推廣與應用。