樊俊屹, 劉春國(guó), 陶志剛

(中國(guó)地震臺(tái)網(wǎng)中心,北京100045)

0 引言

通過(guò)井水位觀測(cè)獲取含水巖體的應(yīng)力、應(yīng)變狀態(tài)變化,是研究地震孕育與發(fā)展過(guò)程的重要手段之一。井水位的動(dòng)態(tài)變化是水巖系統(tǒng)對(duì)外部環(huán)境激勵(lì)的綜合響應(yīng),主要影響因素有構(gòu)造應(yīng)力、氣壓、固體潮、降雨及抽水等[1]。其中構(gòu)造應(yīng)力引起的水位變化是我們觀測(cè)的目標(biāo),氣壓和固體潮長(zhǎng)期作用于井-含水層系統(tǒng),具有自身固有的變化規(guī)律,通常作為水位背景動(dòng)態(tài)的影響因素來(lái)看待,且井水位的響應(yīng)系數(shù)(氣壓系數(shù)和潮汐系數(shù))常用于衡量其應(yīng)力、應(yīng)變的響應(yīng)能力。而降雨與抽水則是水位觀測(cè)的主要干擾因素,與降雨相比,抽水干擾幅度大、無(wú)規(guī)律,其更難被排除。根據(jù)2020年全國(guó)地震地下水位觀測(cè)網(wǎng)干擾因素統(tǒng)計(jì),全國(guó)有多口井受到抽水干擾,其中約有10%的井水位觀測(cè)受抽水干擾較為嚴(yán)重,出現(xiàn)大幅的階變與突跳,增加了從水位觀測(cè)數(shù)據(jù)中提取地震異常信息的難度。

通過(guò)對(duì)仙游井的觀測(cè)系統(tǒng)進(jìn)行排查,并完成其抽水實(shí)驗(yàn)分析,結(jié)果顯示仙游井是受打井、抽水干擾較為嚴(yán)重的觀測(cè)井[2]。這種井水位動(dòng)態(tài)抽水干擾在時(shí)域上的影響是顯著的[3-6],但是在頻域上是如何影響,及其影響下的頻帶差異性以及與應(yīng)力應(yīng)變響應(yīng)能力變化的研究相對(duì)較少。

小波變換具有多分辨率分析的特點(diǎn),在時(shí)頻兩域都具有表征信號(hào)局部特征的能力[7-9],因此本文擬采用連續(xù)小波分析、離散小波分析與調(diào)和分析方法來(lái)分析福建仙游井在打井抽水前后各頻帶的影響程度和水位與氣壓、固體潮響應(yīng)特性的變化,尋找受干擾較少的頻段,以期為受干擾的水位資料分析處理與異常提取提供一種新的思路。

1 仙游井水位觀測(cè)概況

仙游井位于福建莆田市仙游縣楓亭鎮(zhèn)寶坑村,屬非自流井,井深153.37 m。根據(jù)鉆孔資料,0～10.51 m為河相沖積物,10.51 m以下為花崗巖?；◢弾r局部裂隙發(fā)育,賦水性較好,138.51～138.67 m為主要觀測(cè)含水帶。該井自2013年8月起開(kāi)始觀測(cè)水位,觀測(cè)儀器為珠海泰德公司生產(chǎn)的TDL-15水位儀,水位儀的分辨力為1 mm,觀測(cè)精度優(yōu)于2 cm,量程為0～10 m。圖1為其觀測(cè)井柱狀圖。

圖1 觀測(cè)井柱狀圖[2]Fig.1 Bar graph of the observation well[2]

2017年7月5日以前水位動(dòng)態(tài)年變規(guī)律較清晰,年動(dòng)態(tài)總體受降雨控制,每年9—10月出現(xiàn)波峰、5月出現(xiàn)波谷,年變幅度小于1 m。多年動(dòng)態(tài)受區(qū)域性淺層地下水開(kāi)采影響呈現(xiàn)緩慢下降趨勢(shì),水位固體潮效應(yīng)和氣壓效應(yīng)明顯,受干擾前全年平均潮汐系數(shù)為2.3 mm/10-8,觀測(cè)精度達(dá)到0.24%。

2017年7月5日14—17時(shí)水位出現(xiàn)大幅下降,幅度達(dá)到2.5 m;7月6日8—12時(shí)水位出現(xiàn)大幅下降-回升變化,幅度約4 m;7月13—23日水位上升超量程,之后調(diào)整傳感器繼續(xù)正常觀測(cè);10月31日達(dá)到峰值后緩慢下降,期間伴有多次階降。2018年5月以后數(shù)據(jù)突跳、階變頻次及幅度大幅變化減少,動(dòng)態(tài)相對(duì)趨于穩(wěn)定,穩(wěn)定后水位與2017年同期相比上升約6.5 m(如圖2所示)。根據(jù)張清秀等[2]的研究結(jié)果,水位出現(xiàn)大幅上升與2017年7月4—6日在仙游井東南面108 m處新打100 m機(jī)井有關(guān)。新打機(jī)井可能成為上部淺層含水層的涌入通道,使其經(jīng)常出現(xiàn)下降變化。通過(guò)新打機(jī)井的抽水試驗(yàn)來(lái)確定與機(jī)井的不定期、不定量抽水有關(guān)。

圖2 2017—2018年仙游井水位小時(shí)值變化Fig.2 Hourly value variation of water level in Xianyou well (2017—2018)

2 時(shí)頻特征分析

收集仙游井2017年1月1日—2018年12月31日水位整點(diǎn)值數(shù)據(jù)。由于該段數(shù)據(jù)中存在缺數(shù)問(wèn)題,因此在進(jìn)一步開(kāi)展水位時(shí)頻特征分析前,首先的對(duì)水位數(shù)據(jù)進(jìn)行錯(cuò)誤數(shù)據(jù)去除、缺數(shù)插值等預(yù)處理[10],以保證時(shí)頻處理的合理性與完備性。

2.1 總體時(shí)頻特征

小波變換包括連續(xù)小波變換和離散小波變換兩種,兩者的區(qū)別在于尺度參數(shù)是否離散化。小波連續(xù)變換通常使用基數(shù)小于2的指數(shù)尺度,離散小波變換則通常使用基數(shù)為2的指數(shù)尺度。連續(xù)小波變換具有頻帶更為連續(xù)的特點(diǎn),特別適合全局時(shí)頻特征的總體把握。

對(duì)于任意能量有限信號(hào)f(t),其連續(xù)小波變換的定義為:

式中:ψ(t)是母小波或者基本小波;將小波母函數(shù)ψ(t)進(jìn)行伸縮和平移,設(shè)其伸縮因子(尺度因子)為a和平移因子為t,連續(xù)小波變換中a和t均取連續(xù)變化的值,將a和t離散化則為離散小波變換。

要得到水位觀測(cè)的振幅和相位兩方面的信息,就要選擇復(fù)值小波,因?yàn)閺?fù)值小波具有虛部,可以對(duì)相位進(jìn)行很好的表達(dá)[11]。Morlet小波不但具有非正交性而且還是由Gaussian調(diào)節(jié)的指數(shù)復(fù)值小波,因此采用Morlet小波對(duì)2017—2018年全時(shí)段整點(diǎn)值序列進(jìn)行連續(xù)小波變換,獲得全時(shí)段連續(xù)小波功率譜圖[圖3(a)]。從圖3(a)可以看出2017年6月8日—2018年8月17日全頻帶能量均較強(qiáng),特別在低頻帶存在10～30天顯著的長(zhǎng)周期信號(hào)。在干擾發(fā)生前,功率譜上呈現(xiàn)顯著的串珠狀排列的固體潮日波、半日波[圖3(b)]。干擾最強(qiáng)的周期成分能量是正常時(shí)段最強(qiáng)的半日波能量的20倍左右,因此固體潮在全時(shí)段功率譜圖上不明顯。干擾趨緩段[圖3(c)],不再存在能量大于0.1的周期成分,半日波有規(guī)律排列,能量較為顯著;日波起始不清晰,2018年10月7日以后趨于清晰,與正常時(shí)段相比,低頻周期信號(hào)豐富,且出現(xiàn)沒(méi)有規(guī)律,2.15～5.68天、6.53～9.89天、9.89～16.07天的周期成分能量與半日波基本相當(dāng)或略大。

2.2 干擾后效細(xì)節(jié)分析

針對(duì)上述干擾趨于穩(wěn)定后與正常時(shí)段相比,低頻周期信號(hào)豐富,且出現(xiàn)沒(méi)有規(guī)律的現(xiàn)象。本文使用離散小波方法分析抽水給地下水動(dòng)態(tài)帶來(lái)的后效,分別選取正常時(shí)段中的一個(gè)月(2017年1月)和干擾趨緩段中的一個(gè)月(2018年12月)的觀測(cè)數(shù)據(jù)來(lái)分析。

對(duì)這兩個(gè)月的水位整點(diǎn)值數(shù)據(jù)進(jìn)行小波分解,小波函數(shù)為db4、分解為8層[7]。正常期和趨緩期地下水位小波分解后細(xì)節(jié)部分對(duì)比曲線如圖4所示。由圖可見(jiàn):在中頻帶4～32 h[圖4(b)、(c)、(d)]趨緩期和正常期的振幅幅度一致,峰值時(shí)間不同,存在明顯日波與半日波,但是兩個(gè)時(shí)期相位略有改變。在高頻帶[圖4(a)]存在多個(gè)強(qiáng)信號(hào),出現(xiàn)沒(méi)有規(guī)律;在低頻帶[圖4(e)、(f)、(g)、(h)]仍然存在多個(gè)強(qiáng)信號(hào),幅度差異較大,趨穩(wěn)期振幅是正常期振幅的2～4倍。

圖4 正常期和趨緩期地下水位小波分解后主要波動(dòng)頻段Fig.4 Main fluctuation frequency bands of ground water level after wavelet decomposition in normal and slowdown periods

3 氣壓與固體潮響應(yīng)特性

3.1 水位與氣壓、固體潮相關(guān)

抽水干擾趨穩(wěn)后,抽水干擾仍對(duì)水位的低頻帶存在影響,但仙游井水位對(duì)固體潮、氣壓的響應(yīng)特性是否會(huì)發(fā)生變化?

選擇正常期(2017年1月)、干擾期(2017年10月)、恢復(fù)期(2018年12月)為代表,同時(shí)獲取該井同期觀測(cè)的氣壓整點(diǎn)值數(shù)據(jù),利用Mapsis2000軟件計(jì)算同期重力理論固體潮整點(diǎn)值數(shù)據(jù)。采用db4小波函數(shù),將水位、氣壓、理論固體潮分解為8層[12],分別計(jì)算每一層細(xì)節(jié)部分水位與氣壓、水位與固體潮的相關(guān)系數(shù),繪制正常期、干擾期、趨緩期水位與氣壓、水位與固體潮不同頻段相關(guān)系數(shù)變化如圖5所示。

圖5 正常期、干擾期及趨緩期的水位與固體潮、氣壓的頻域相關(guān)性Fig.5 Frequency domain correlation of water level with earth tide and pressure in normal period, disturbance period, and slowdown period

正常期2～32 h頻帶,水位與固體潮高度相關(guān),相關(guān)系數(shù)0.8～1.0;32 h以上周期成分相關(guān)性急劇下降,64～256 h基本不相關(guān);氣壓正好相反,2～32 h基本不相干,而64～256 h則顯著相關(guān)。

干擾期水位與固體潮相關(guān)系數(shù)變化曲線上出現(xiàn)兩個(gè)小峰,同時(shí)它與氣壓的相關(guān)系數(shù)變化呈現(xiàn)與正常期相反的變化趨勢(shì)。水位與固體潮、氣壓在整個(gè)頻段的相關(guān)系數(shù)分別低于0.2和0.3,顯示二者關(guān)聯(lián)性較差。

趨緩期水位固體潮的相關(guān)系數(shù)變化形態(tài)又恢復(fù)到正常期的變化形態(tài)。2～4 h水位與固體潮相關(guān)不顯著,8～32 h水位與固體潮相關(guān)較為顯著,但與正常期相比偏低,64～256 h基本不相關(guān);與氣壓的相關(guān)系數(shù)仍延續(xù)了干擾期的變化規(guī)律,基本不相關(guān)。

綜上所述,仙游井在打井及抽水干擾后,水位動(dòng)態(tài)趨于穩(wěn)定,水位對(duì)固體潮在8～32 h頻段的效應(yīng)逐漸恢復(fù),而氣壓效應(yīng)則基本消失。這可能與恢復(fù)期仍存在大量的低頻成分,而水位對(duì)氣壓的響應(yīng)主要在低頻段有關(guān)。

3.2 M2潮汐參數(shù)變化

以上水位與固體潮的相關(guān)性分析表明,干擾前后水位對(duì)8～32 h的固體潮周期成分響應(yīng)能量變化總體上較為穩(wěn)定。選擇固體潮分波中最為顯著的諧波M2水位響應(yīng)潮汐系數(shù)和相位來(lái)研究其在整個(gè)時(shí)段的變化。

利用滑動(dòng)Venidicov調(diào)和分析方法,滑動(dòng)步長(zhǎng)1天,窗長(zhǎng)30天,獲得M2波潮汐因子、相位滯后序列。M2波觀測(cè)精度為M2波潮汐因子計(jì)算中誤差與M2波潮汐因子的比值,可以反映水位M2波潮汐觀測(cè)的干擾程度。去除干擾強(qiáng)烈段,繪制2017年1月1日—2018年12月31日M2潮汐因子及相位滯后隨時(shí)間變化對(duì)比曲線(圖6)。從圖中可以看出,7—8月觀測(cè)精度大于0.5,最大達(dá)到3.42,水位M2波潮汐信號(hào)淹沒(méi)在強(qiáng)干擾信號(hào)中,計(jì)算的潮汐因子和相位滯后沒(méi)有參考價(jià)值,2018年下半年緩和期M2波潮觀測(cè)精度較好,但仍不及未打井前的1/10。

圖6 井水位M2波觀測(cè)振幅與相位隨時(shí)間變化Fig.6 Variation of observation amplitude and phase of M2 wave of well water level with time

從圖6可以看出,水位M2潮汐相位波動(dòng)比振幅波動(dòng)顯著,反映了觀測(cè)含水層花崗巖裂隙水非均質(zhì)特性。打井抽水前M2潮汐振幅變化很小,打井后一直到2018年5月井水位M2潮汐振幅及相位均處于劇烈變化階段,總體上大幅上升、下降變化基本同步,部分時(shí)段出現(xiàn)不同步現(xiàn)象,初步分析認(rèn)為打井及抽水引發(fā)了垂向流[13]。干擾趨緩期振幅及相位波動(dòng)接近打井前水平,相位略有抬升、振幅有所下降,反映了打井后淺部地下水混入,原觀測(cè)含水系統(tǒng)的封閉性變差,潮汐應(yīng)變響應(yīng)能力有所下降。

4 結(jié)語(yǔ)

采用多種時(shí)頻分析的方法,對(duì)仙游井水位打井抽水干擾前后的時(shí)頻特征變化和氣壓與固體潮效應(yīng)特性變化的同步時(shí)頻分析,結(jié)果顯示:

(1) 打井/抽水干擾嚴(yán)重時(shí)段,井水位動(dòng)態(tài)在時(shí)域上出現(xiàn)大幅階變和突跳變化,在頻域上亦影響顯著,抽水干擾信號(hào)完全掩蓋地下水位的固有周期信號(hào)及對(duì)氣壓、固體潮的響應(yīng)效應(yīng)。

(2) 打井/抽水干擾趨于緩和,中高頻段(8～23 h)基本修復(fù),而低頻段(2～8 h)仍然存在多個(gè)強(qiáng)信號(hào)(幾倍于之前),顯示打井抽水干擾后,地下水含水系統(tǒng)的修復(fù)一個(gè)緩慢的過(guò)程。

(3) 井水位M2潮汐振幅及相位變化分析表明,打井/抽水引發(fā)井-含水系統(tǒng)的以徑向?yàn)橹鞔瓜驗(yàn)檩o的排水,干擾趨緩后,井水位M2潮汐振幅及相位波動(dòng)趨緩,井水位M2波觀測(cè)精度與打井前相比有所下降,潮汐應(yīng)變響應(yīng)能力略有下降。

綜上所述,打井抽水干擾對(duì)井水位的高、中、低頻段影響是不相同的,其影響具有一定的后效性,在水位資料異常識(shí)別與分析中應(yīng)選取干擾較少的頻段來(lái)進(jìn)行分析,對(duì)于仙游井水位來(lái)說(shuō)中頻段是受影響最小、最有益于開(kāi)展監(jiān)測(cè)預(yù)報(bào)工作的頻段。