地層形變與地下水位分層監(jiān)測數(shù)據(jù)的交叉小波分析

孔祥如

(北京市地質(zhì)環(huán)境監(jiān)測所，北京 100195)

地面沉降是指由于自然或人為因素引發(fā)的地下松散巖層固結(jié)壓縮，并導(dǎo)致一定區(qū)域范圍內(nèi)地面高程降低的地質(zhì)現(xiàn)象，是一種緩變型地質(zhì)災(zāi)害[1]，地下水超量開采是城市中發(fā)生地面沉降的主要原因[2]。L.G.Devin等[3]研究記錄了中國、意大利、美國和墨西哥等近90 個國家因超采地下水造成不同程度的地面沉降。我國也有超100 座城市遭受地面沉降威脅，嚴(yán)重沉降區(qū)面積為7.9 萬km2(累積沉降量＞200 mm)[4-5]。北京平原區(qū)自20 世紀(jì)50 年代開始出現(xiàn)地下水過度開采現(xiàn)象，導(dǎo)致地下水位持續(xù)下降，逐漸使北京發(fā)展成為地面沉降嚴(yán)重地區(qū)[6]。根據(jù)研究資料顯示，北京地面沉降仍處于快速發(fā)展階段[7]，截至2016 年，北京平原區(qū)最大累計沉降量達(dá)1 864 mm，最大年沉降速率超過110 mm/a[8]，地面沉降已經(jīng)成為北京市經(jīng)濟(jì)社會發(fā)展的主要制約因素之一。

長期超量開采地下水導(dǎo)致地下水位大幅下降，鄰近黏性土層向含水層釋水，孔隙水壓力消散，導(dǎo)致黏性土層壓縮，形成地面沉降[9]。受壓縮層的巖性、厚度和固結(jié)程度等因素控制，從水位下降到地面沉降發(fā)生這一過程具有一定滯后性[10]。對于滯后性的研究一直是地面沉降研究領(lǐng)域的重點(diǎn)。張建偉等[11]分析了上海城區(qū)地下水采灌量、開采層次及地面沉降變化特征，通過數(shù)學(xué)統(tǒng)計方法得出地面沉降與地下水開采量存在明顯正相關(guān)，并具有時滯效應(yīng)。施小清等[12]通過分析常州市分層標(biāo)及水位觀測孔數(shù)據(jù)并結(jié)合壓縮試驗，發(fā)現(xiàn)常州地區(qū)含水砂層在地面沉降過程中同樣存在變形滯后現(xiàn)象，并揭示了其蠕變特性。徐進(jìn)等[13]基于Biot 理論的地面沉降耦合模型，利用半解析數(shù)值原理和黏彈性流變理論，推導(dǎo)了可壓縮土層黏彈性耦合變形的求解格式，反映土體黏滯性所導(dǎo)致的變形滯后效應(yīng)。Zhou Zhifang 等[14]提出了釋水滯后率的概念用來描述弱透水層中的滯后釋水規(guī)律，其研究表明弱透水層的水動態(tài)響應(yīng)滯后效應(yīng)不僅與弱透水層的性質(zhì)有關(guān)，也與弱透水層厚度的平方成正比關(guān)系。已有研究主要側(cè)重于土層的滯后機(jī)理分析，對滯后時間的研究也主要是通過室內(nèi)試驗或數(shù)值模擬得出理論或預(yù)測結(jié)果，對實(shí)際監(jiān)測數(shù)據(jù)中的準(zhǔn)確滯后時間缺乏有效的識別方法。

小波變換是一種數(shù)據(jù)時頻局部化分析工具，由其衍生的交叉小波變換能有效識別出兩組序列的周期相關(guān)性，在天文、氣象、地學(xué)等領(lǐng)域已取得一些應(yīng)用成果[15-17]，在地下水學(xué)科方面，主要用于分析大氣降水與地下水位的周期性變化特征[18-20]，而對于地面沉降尤其是分層沉降與地下水位的時滯響應(yīng)方面缺乏相關(guān)研究應(yīng)用。筆者選取北京平原區(qū)地面沉降監(jiān)測站地下水位與分層標(biāo)的長時間序列監(jiān)測數(shù)據(jù)，采用Mann-Kendall趨勢檢驗、連續(xù)小波變換、交叉小波變換等方法，對不同深度的地層垂向形變量和水位變化數(shù)據(jù)進(jìn)行處理和對比分析，獲取地層形變和地下水位動態(tài)的周期規(guī)律，繼而進(jìn)行交叉小波變換定量分析不同埋藏深度下的地層壓縮-反彈對地下水位動態(tài)響應(yīng)的滯后特征，以期進(jìn)一步認(rèn)識地面沉降成因機(jī)理，提高地面沉降預(yù)測模型精度。

1 研究區(qū)概況

北京平原位于華北平原北部，北依軍都山，西靠太行山余脈西山。在大地構(gòu)造上屬于中朝準(zhǔn)地臺(Ⅰ)華北斷陷拗(Ⅱ)，中生代以來，受印支、燕山、喜馬拉雅多期構(gòu)造活動作用，進(jìn)一步形成隆起-凹陷相間的次級構(gòu)造單元。北京平原第四系屬典型山前沖洪積地層，整體上由永定河、潮白河、溫榆河、泃河和拒馬河五大水系的沖洪積扇群和多個沉積凹陷構(gòu)成，第四紀(jì)沉積物巖性橫向變化迅速，自山前向平原沉積物顆粒直徑逐漸減小，在垂向上，則呈礫石層、砂層、黏土層交互出現(xiàn)的互層結(jié)構(gòu)。第四系孔隙水是研究區(qū)地下水主要開采層。根據(jù)地層年代、巖性、埋藏條件及地下水補(bǔ)徑排條件，通常將北京平原第四系含水層概化為3 組：第一組含水層(Ⅰ)為潛水和淺層承壓水，對應(yīng)地層為全新統(tǒng)(Q4)和上更新統(tǒng)(Q3)；第二組含水層(Ⅱ)為中深層承壓水，對應(yīng)地層為中更新統(tǒng)(Q2)；第三含水層(Ⅲ)為深層承壓水，對應(yīng)地層為下更新統(tǒng)(Q1)。每一組含水層多夾雜數(shù)層以粉質(zhì)黏土、黏土或粉土為主的弱透水層(圖1)。由于地下水超量開采導(dǎo)致含水層水位下降，弱透水層向含水層釋水使得弱透水層孔隙水壓力降低，進(jìn)而引發(fā)黏性土層壓縮，是該地區(qū)發(fā)生地面沉降的主要原因。

圖1 北京平原區(qū)典型第四系地質(zhì)剖面Fig.1 Typical Quaternary geological section in the Beijing Plain

2 研究方法

2.1 Mann-Kendall 趨勢檢驗

Mann-Kendall 趨勢檢驗是一種非參數(shù)統(tǒng)計檢驗方法，具有樣本不必遵從某種特定分布、結(jié)果也不受少數(shù)異常值干擾的特點(diǎn)，適用于長時間序列的趨勢分析，被世界氣象組織推薦并廣泛應(yīng)用，在水文學(xué)領(lǐng)域也積累了大量研究成果[21-24]。

設(shè)xn(n=1,2,···,N，N為樣本數(shù)量)為平穩(wěn)獨(dú)立序列，統(tǒng)計量U定義為：

當(dāng)n≥10，統(tǒng)計量U近似服從正態(tài)分布，不考慮序列中等值數(shù)據(jù)點(diǎn)情況，其均值E(U)=0。

標(biāo)準(zhǔn)化的檢驗統(tǒng)計量ZMK采用下式計算：

式中：g、h為任意樣本序號，且g＞h。當(dāng)β為正值，表示上升趨勢；當(dāng)β為負(fù)值，表示下降趨勢[19]。

2.2 連續(xù)小波變換

小波變換將波頻顯示在不同時間尺度上，來反映長時間序列中不同尺度顯著周期及其所在的時間段。本文選用在時間和頻率的局部化上具有良好平衡性的Morlet 小波作為母小波φ0，函數(shù)為：

式中：ω0為 無量綱頻率；η為無量綱時間。

時間序列的連續(xù)小波變換可定義為對序列xn基于小波標(biāo)準(zhǔn)化算法的卷積，函數(shù)為：

式中：δt為均一時間步長；s為小波尺度。

背景功率譜采用紅噪聲檢驗，紅噪聲檢驗過程用一階自回歸方程。背景紅噪聲功率譜Pk定義為：

式中：r為紅噪聲功率譜中自回歸方程的相關(guān)系數(shù)；k為傅里葉頻率系數(shù)。

2.3 交叉小波變換

交叉小波變換能夠分析兩組序列在時頻空間中能量共振和協(xié)方差分布規(guī)律，可以在不同時間尺度上揭示兩組序列頻率周期的一致性和相關(guān)性，并能通過小波相位角分析，計算時頻空間中的相位關(guān)系。

式中：σX、σY分別為時間序列Xn和Yn的標(biāo)準(zhǔn)差；當(dāng)背景功率譜為實(shí)小波時，v=1，為復(fù)小波時v=2；Zv(p)是概率p的置信度水平。

交叉小波的相位角am定義為：

式中：ai(i=1,2,···,N)為小波變換時間域上的單個相位角。

小波相位角的標(biāo)準(zhǔn)偏差為：

本文所采用的小波分析計算方法和程序代碼主要參考C.Torrence 等[25]和A.Grinsted 等[26]的研究成果。

3 數(shù)據(jù)來源與預(yù)處理

3.1 數(shù)據(jù)來源

本文使用的地面沉降和地下水水位數(shù)據(jù)來自北京平各莊地面沉降監(jiān)測站，該站建于2008 年，站內(nèi)埋設(shè)有分層標(biāo)和地下水位監(jiān)測井，監(jiān)測層位包含北京平原區(qū)第四系主要含水層(砂層)和弱透水層(黏性土層)，地層結(jié)構(gòu)具有代表性。本文使用的數(shù)據(jù)時間段為2008 年6 月-2018 年12 月，監(jiān)測頻率為1 次/d。地下水位監(jiān)測井(D-1-D-5)主要監(jiān)測不同含水層水位，地下水位時間序列(D1-D5)為單井監(jiān)測數(shù)據(jù)；分層標(biāo)監(jiān)測井(F-1-F-5)監(jiān)測2 個分層標(biāo)之間的地層，包括弱透水層和含水層，其地層形變時間序列(F1-F5)為地層頂?shù)? 個分層標(biāo)監(jiān)測數(shù)據(jù)之差，即監(jiān)測層位的垂向形變量。地下水位和沉降監(jiān)測層的基本情況見表1。

表1 地下水位和沉降監(jiān)測層基本情況Table 1 Basic information of groundwater monitoring wells and subsidence monitoring layers

3.2 數(shù)據(jù)預(yù)處理

地下水長期持續(xù)開采、地應(yīng)力載荷等長時間尺度因素會對數(shù)據(jù)分析產(chǎn)生干擾，因此首先對地下水位和地層形變數(shù)據(jù)進(jìn)行Mann-Kendall 趨勢檢驗，以判定其變化趨勢。結(jié)果見表2，可以看出，各時間序列的ZMK值均為負(fù)值且均小于-3.3，為下降趨勢，ZMK通過了α=0.001水平的顯著性檢驗，下降趨勢顯著。地下水位方面，D-1 和D-2 井監(jiān)測層位的深部承壓含水層為地下水主要開采層，2008-2018 年地下水位均呈持續(xù)下降趨勢，地下水位年均下降速率分別為1.780 3 m/a和1.184 9 m/a；D-3 中部承壓含水層在2014 年以前水位持續(xù)下降，2014 年以后受地下水禁限采政策影響開采量逐漸減少，該含水層水位出現(xiàn)緩慢回升，2008-2018 年地下水位年均下降速率為0.252 1 m/a；D-4 和D-5 井監(jiān)測層位的淺部含水層地下水開采量較小，主要受大氣降水影響，10 a 平均下降速率分別為0.416 7 m/a和0.491 7 m/a。地層形變方面，F(xiàn)-1 地層自監(jiān)測以來持續(xù)發(fā)生壓縮沉降，由于下更新統(tǒng)土體固結(jié)程度較高，可壓縮性較低，屬于地面沉降次要貢獻(xiàn)層，10 a 平均沉降速率為1.772 8 mm/a；F-2 與F-3 地層是地面沉降的主要貢獻(xiàn)層，兩者沉降量占地面總沉降量的80%以上，其中F-2 地層10 a 平均沉降速率為13.046 9 mm/a，2014 年以后受地下水位回升影響，F(xiàn)-3 地層沉降速率有所減緩，但仍然保持沉降，10 a 平均沉降速率為

表2 Mann-Kendall 趨勢檢驗 ZMK 和 β值Table 2 ZMK and β of the Mann-Kendall trend test

7.158 7 mm/a；F-4 和F-5 地層屬于弱沉降層，地層沉降量較小，個別年份還會有所反彈，10 a 平均沉降速率分別為0.866 4 mm/a 和0.709 0 mm/a。為了消除變化趨勢的影響，根據(jù)R.T.Hanson 等[27]的數(shù)據(jù)前處理方法，將地下水位和地層形變數(shù)據(jù)減去三次多項式回歸方程，取其殘差作為分析地下水位和地層形變關(guān)系的時間序列數(shù)據(jù)。

4 結(jié)果與分析

4.1 連續(xù)小波變換

將全部的地下水位和地層形變量時間序列分別進(jìn)行連續(xù)小波變換(圖2)，以分析其周期性波動特征和顯著周期時段等信息。

地下水位方面(圖2a)，埋深81 m 以下的中-深層承壓水時間序列D1、D2 和D3 存在1 a 左右的主震蕩周期，連續(xù)分布且周期顯著；60 m 以淺含水層地下水位在大部分時域上周期性不顯著，其中淺層承壓水時間序列D4 僅在2008-2012 年存在1 a 左右的震蕩周期，潛水時間序列D5 僅在2014-2015 年存在0.15～0.75 a 的主震蕩周期。2000 年以后，由于淺部含水層水位下降和水質(zhì)惡化，北京平原區(qū)地下水開采深度不斷增加，中-深層承壓水成為主要開采地下水[28]，尤其在農(nóng)作物灌溉期地下水超采現(xiàn)象嚴(yán)重，造成深部含水層水位具有明顯的季節(jié)性波動特征[29]，是深部含水層較淺部含水層水位震蕩周期更為顯著的主要原因。

圖2 各時間序列連續(xù)小波變換Fig.2 Continuous wavelet transforms of the time series

地層形變方面(圖2b)，根據(jù)監(jiān)測數(shù)據(jù)埋深在63～209 m 的地層沉降占比在80% 以上，屬于沉降主要貢獻(xiàn)層，其中地層形變時間序列F2 和F3 在全時域存在1 a 左右的主震蕩周期，且通過95%紅噪聲檢驗的區(qū)域與同層水位具有相似形態(tài)，表明該嚴(yán)重沉降層的地下水位波動可能與地層形變量存在時頻域相關(guān)；埋深超過209 m 的地層形變量相對較小，其時間序列F1 在2010-2017 年存在1 a 左右的主震蕩周期，通過95%紅噪聲檢驗的區(qū)域與同層水位較為相似；埋深小于63 m 的地層形變量最小屬于弱沉降層，其中時間序列F4 在2010-2017 年存在1 a 左右的主震蕩周期，時間序列F5 在2013-2015 年存在0.5～3.0 a 的主震蕩周期，與同層水位小波變換形態(tài)相比均具有更顯著的周期性。

4.2 交叉小波變換

對相同層位的地下水位和地層形變量時間序列進(jìn)行交叉小波變換(圖3)，以分析其共振周期、顯著時段和相位關(guān)系等，研究地面沉降對地下水位變化的時滯響應(yīng)特征。

由圖3a-圖3c 可知，埋深在63 m 以下的地層地下水位和地層形變量通過95%紅噪聲檢驗且共振周期在0.75～1.29 a 的區(qū)域能量密度最高，顯著時段連續(xù)分布整個序列時域，表明兩者共振關(guān)系顯著。滯后性方面，埋深63～209 m 地層形變量相對地下水位的遲滯時間相近，分別為(7.16±7.09) d 和(9.66±6.62) d(表3)；埋深209 m 以下地層形變滯后時間更長，為(16.58±8.91) d。由此可以看出，深部主要沉降層的地層形變與地下水位動態(tài)相干性較強(qiáng)，形變周期相對地下水位周期具明顯的滯后性，但不同層位滯后時間具有差異性。

圖3 地下水位與地層形變量交叉小波變換Fig.3 Cross wavelet transforms of groundwater levels and stratum deformation values

由圖3d 和圖3e 可知，埋深在63 m 以淺的地層地下水位和地層形變量通過95%紅噪聲檢驗的區(qū)域在時頻域上不連續(xù)，在部分時域上具有顯著共振周期。埋深32～63 m 地層形變滯后時間為(21.01±21.06) d(表3)，相較于深部主要沉降層有更長的滯后時間；埋深2～32 m 地層交叉相位的相位誤差值遠(yuǎn)大于平均相位值，即在共振周期顯著的時頻域里形變滯后時間普遍存在負(fù)值，說明淺層地層的形變量與潛水水位變化不具有可信的時滯關(guān)系。

表3 地下水位和地層形變量交叉小波變換周期與相位特征Table 3 Periodicities and phases of cross-wavelet transforms of groundwater levels and stratum deformation values

為進(jìn)一步驗證埋深63 m 以淺地層形變是否與不同層位含水層地下水位存在遲滯關(guān)系，對埋深63 m 以淺地層的地下水位和形變量進(jìn)行錯層交叉小波變換(圖4)。

圖4 淺部錯層地下水位與形變量交叉小波變換Fig.4 Cross wavelet transforms of groundwater levels and stratum deformation values in shallow staggered layers

圖4a 看出，埋深在81～119 m 含水層地下水位與埋深在32～63 m 地層形變量在全時域具有0.74～1.25 a的顯著共振周期，形變時滯為(32.02±9.67) d(表4)，說明中層承壓水的水位波動對上部鄰近地層的形變也具有一定影響，但形變滯后時間較長，引發(fā)的形變量較小；圖4b-圖4c 看出，埋深在60 m 以淺地層的含水層地下水位與相鄰地層形變量通過95%紅噪聲檢驗的區(qū)域形變時滯分別為(3.90±64.87) d 和(40.87±57.38) d(表4)，表明潛水和淺層承壓水與相鄰地層的形變不具有明確的時滯關(guān)系。

表4 淺部錯層地下水位和形變量交叉小波變換周期與相位特征Table 4 Periodicities and phases of cross-wavelet transforms of groundwater levels and stratum deformation values in shallow staggered layers

5 結(jié) 論

a.研究區(qū)作為地下水主采層的中-深層承壓水含水層水位具有1 a 左右的主震蕩周期，在全時域連續(xù)分布且周期顯著，符合地下水位季節(jié)性波動的規(guī)律特征；潛水和淺部承壓水含水層水位僅部分年份通過95%紅噪聲檢驗，在大部分時域無顯著周期。

b.深部嚴(yán)重沉降層的形變量具有1 a 左右的主震蕩周期，在全時域連續(xù)分布且周期顯著，地下水位與形變量通過95%紅噪聲檢驗的共振周期為0.75～1.29 a且連續(xù)分布，驗證了兩者具有可靠相關(guān)性，地層形變滯后于水位變化，地層由淺到深形變滯后時間分別為(16.58±8.91) d、(7.16±7.09) d 和(9.66±6.62) d。

c.淺部弱沉降層中，埋深在2～32 m 地層形變量在COI 區(qū)域內(nèi)存在0.5～3.0 a 的主震蕩周期，與地下水位不具有時滯關(guān)系；埋深在32～63 m 地層形變量具有1 a 左右的主震蕩周期，與第一組含水層水位不具有時滯關(guān)系，而與中層承壓水含水層水位存在顯著共振周期，形變時滯為(32.02±9.67) d。

d.運(yùn)用小波變換方法分析北京平各莊地面沉降監(jiān)測站不同深度地層的形變滯后特征，精細(xì)刻畫地面沉降與地下水的相互關(guān)系，為構(gòu)建地面沉降精細(xì)化模型、提高地面沉降預(yù)測精度以及研究更有效的地面沉降防控措施提供新的技術(shù)思路。但地面沉降的發(fā)展受巖性、厚度、地應(yīng)力和固結(jié)程度等多種因素的控制，地面沉降滯后性研究還需大量數(shù)據(jù)的支撐。后續(xù)將開展其他監(jiān)測站的數(shù)據(jù)分析，擴(kuò)大時間尺度，用更充足的數(shù)據(jù)去論證地面沉降對地下水的響應(yīng)特征；另一方面將加強(qiáng)對研究區(qū)鉆探、物化探及相關(guān)試驗數(shù)據(jù)的搜集，結(jié)合小波變換分析結(jié)果，進(jìn)一步研究不同環(huán)境因素對地面沉降滯后性的影響。