樊俊屹, 劉春國, 陶志剛
(中國地震臺網中心,北京100045)
通過井水位觀測獲取含水巖體的應力、應變狀態變化,是研究地震孕育與發展過程的重要手段之一。井水位的動態變化是水巖系統對外部環境激勵的綜合響應,主要影響因素有構造應力、氣壓、固體潮、降雨及抽水等[1]。其中構造應力引起的水位變化是我們觀測的目標,氣壓和固體潮長期作用于井-含水層系統,具有自身固有的變化規律,通常作為水位背景動態的影響因素來看待,且井水位的響應系數(氣壓系數和潮汐系數)常用于衡量其應力、應變的響應能力。而降雨與抽水則是水位觀測的主要干擾因素,與降雨相比,抽水干擾幅度大、無規律,其更難被排除。根據2020年全國地震地下水位觀測網干擾因素統計,全國有多口井受到抽水干擾,其中約有10%的井水位觀測受抽水干擾較為嚴重,出現大幅的階變與突跳,增加了從水位觀測數據中提取地震異常信息的難度。
通過對仙游井的觀測系統進行排查,并完成其抽水實驗分析,結果顯示仙游井是受打井、抽水干擾較為嚴重的觀測井[2]。這種井水位動態抽水干擾在時域上的影響是顯著的[3-6],但是在頻域上是如何影響,及其影響下的頻帶差異性以及與應力應變響應能力變化的研究相對較少。
小波變換具有多分辨率分析的特點,在時頻兩域都具有表征信號局部特征的能力[7-9],因此本文擬采用連續小波分析、離散小波分析與調和分析方法來分析福建仙游井在打井抽水前后各頻帶的影響程度和水位與氣壓、固體潮響應特性的變化,尋找受干擾較少的頻段,以期為受干擾的水位資料分析處理與異常提取提供一種新的思路。
仙游井位于福建莆田市仙游縣楓亭鎮寶坑村,屬非自流井,井深153.37 m。根據鉆孔資料,0~10.51 m為河相沖積物,10.51 m以下為花崗巖。花崗巖局部裂隙發育,賦水性較好,138.51~138.67 m為主要觀測含水帶。該井自2013年8月起開始觀測水位,觀測儀器為珠海泰德公司生產的TDL-15水位儀,水位儀的分辨力為1 mm,觀測精度優于2 cm,量程為0~10 m。圖1為其觀測井柱狀圖。

圖1 觀測井柱狀圖[2]Fig.1 Bar graph of the observation well[2]
2017年7月5日以前水位動態年變規律較清晰,年動態總體受降雨控制,每年9—10月出現波峰、5月出現波谷,年變幅度小于1 m。多年動態受區域性淺層地下水開采影響呈現緩慢下降趨勢,水位固體潮效應和氣壓效應明顯,受干擾前全年平均潮汐系數為2.3 mm/10-8,觀測精度達到0.24%。
2017年7月5日14—17時水位出現大幅下降,幅度達到2.5 m;7月6日8—12時水位出現大幅下降-回升變化,幅度約4 m;7月13—23日水位上升超量程,之后調整傳感器繼續正常觀測;10月31日達到峰值后緩慢下降,期間伴有多次階降。2018年5月以后數據突跳、階變頻次及幅度大幅變化減少,動態相對趨于穩定,穩定后水位與2017年同期相比上升約6.5 m(如圖2所示)。根據張清秀等[2]的研究結果,水位出現大幅上升與2017年7月4—6日在仙游井東南面108 m處新打100 m機井有關。新打機井可能成為上部淺層含水層的涌入通道,使其經常出現下降變化。通過新打機井的抽水試驗來確定與機井的不定期、不定量抽水有關。

圖2 2017—2018年仙游井水位小時值變化Fig.2 Hourly value variation of water level in Xianyou well (2017—2018)
收集仙游井2017年1月1日—2018年12月31日水位整點值數據。由于該段數據中存在缺數問題,因此在進一步開展水位時頻特征分析前,首先的對水位數據進行錯誤數據去除、缺數插值等預處理[10],以保證時頻處理的合理性與完備性。
小波變換包括連續小波變換和離散小波變換兩種,兩者的區別在于尺度參數是否離散化。小波連續變換通常使用基數小于2的指數尺度,離散小波變換則通常使用基數為2的指數尺度。連續小波變換具有頻帶更為連續的特點,特別適合全局時頻特征的總體把握。
對于任意能量有限信號f(t),其連續小波變換的定義為:
式中:ψ(t)是母小波或者基本小波;將小波母函數ψ(t)進行伸縮和平移,設其伸縮因子(尺度因子)為a和平移因子為t,連續小波變換中a和t均取連續變化的值,將a和t離散化則為離散小波變換。
要得到水位觀測的振幅和相位兩方面的信息,就要選擇復值小波,因為復值小波具有虛部,可以對相位進行很好的表達[11]。Morlet小波不但具有非正交性而且還是由Gaussian調節的指數復值小波,因此采用Morlet小波對2017—2018年全時段整點值序列進行連續小波變換,獲得全時段連續小波功率譜圖[圖3(a)]。從圖3(a)可以看出2017年6月8日—2018年8月17日全頻帶能量均較強,特別在低頻帶存在10~30天顯著的長周期信號。在干擾發生前,功率譜上呈現顯著的串珠狀排列的固體潮日波、半日波[圖3(b)]。干擾最強的周期成分能量是正常時段最強的半日波能量的20倍左右,因此固體潮在全時段功率譜圖上不明顯。干擾趨緩段[圖3(c)],不再存在能量大于0.1的周期成分,半日波有規律排列,能量較為顯著;日波起始不清晰,2018年10月7日以后趨于清晰,與正常時段相比,低頻周期信號豐富,且出現沒有規律,2.15~5.68天、6.53~9.89天、9.89~16.07天的周期成分能量與半日波基本相當或略大。
針對上述干擾趨于穩定后與正常時段相比,低頻周期信號豐富,且出現沒有規律的現象。本文使用離散小波方法分析抽水給地下水動態帶來的后效,分別選取正常時段中的一個月(2017年1月)和干擾趨緩段中的一個月(2018年12月)的觀測數據來分析。
對這兩個月的水位整點值數據進行小波分解,小波函數為db4、分解為8層[7]。正常期和趨緩期地下水位小波分解后細節部分對比曲線如圖4所示。由圖可見:在中頻帶4~32 h[圖4(b)、(c)、(d)]趨緩期和正常期的振幅幅度一致,峰值時間不同,存在明顯日波與半日波,但是兩個時期相位略有改變。在高頻帶[圖4(a)]存在多個強信號,出現沒有規律;在低頻帶[圖4(e)、(f)、(g)、(h)]仍然存在多個強信號,幅度差異較大,趨穩期振幅是正常期振幅的2~4倍。

圖4 正常期和趨緩期地下水位小波分解后主要波動頻段Fig.4 Main fluctuation frequency bands of ground water level after wavelet decomposition in normal and slowdown periods
抽水干擾趨穩后,抽水干擾仍對水位的低頻帶存在影響,但仙游井水位對固體潮、氣壓的響應特性是否會發生變化?
選擇正常期(2017年1月)、干擾期(2017年10月)、恢復期(2018年12月)為代表,同時獲取該井同期觀測的氣壓整點值數據,利用Mapsis2000軟件計算同期重力理論固體潮整點值數據。采用db4小波函數,將水位、氣壓、理論固體潮分解為8層[12],分別計算每一層細節部分水位與氣壓、水位與固體潮的相關系數,繪制正常期、干擾期、趨緩期水位與氣壓、水位與固體潮不同頻段相關系數變化如圖5所示。

圖5 正常期、干擾期及趨緩期的水位與固體潮、氣壓的頻域相關性Fig.5 Frequency domain correlation of water level with earth tide and pressure in normal period, disturbance period, and slowdown period
正常期2~32 h頻帶,水位與固體潮高度相關,相關系數0.8~1.0;32 h以上周期成分相關性急劇下降,64~256 h基本不相關;氣壓正好相反,2~32 h基本不相干,而64~256 h則顯著相關。
干擾期水位與固體潮相關系數變化曲線上出現兩個小峰,同時它與氣壓的相關系數變化呈現與正常期相反的變化趨勢。水位與固體潮、氣壓在整個頻段的相關系數分別低于0.2和0.3,顯示二者關聯性較差。
趨緩期水位固體潮的相關系數變化形態又恢復到正常期的變化形態。2~4 h水位與固體潮相關不顯著,8~32 h水位與固體潮相關較為顯著,但與正常期相比偏低,64~256 h基本不相關;與氣壓的相關系數仍延續了干擾期的變化規律,基本不相關。
綜上所述,仙游井在打井及抽水干擾后,水位動態趨于穩定,水位對固體潮在8~32 h頻段的效應逐漸恢復,而氣壓效應則基本消失。這可能與恢復期仍存在大量的低頻成分,而水位對氣壓的響應主要在低頻段有關。
以上水位與固體潮的相關性分析表明,干擾前后水位對8~32 h的固體潮周期成分響應能量變化總體上較為穩定。選擇固體潮分波中最為顯著的諧波M2水位響應潮汐系數和相位來研究其在整個時段的變化。
利用滑動Venidicov調和分析方法,滑動步長1天,窗長30天,獲得M2波潮汐因子、相位滯后序列。M2波觀測精度為M2波潮汐因子計算中誤差與M2波潮汐因子的比值,可以反映水位M2波潮汐觀測的干擾程度。去除干擾強烈段,繪制2017年1月1日—2018年12月31日M2潮汐因子及相位滯后隨時間變化對比曲線(圖6)。從圖中可以看出,7—8月觀測精度大于0.5,最大達到3.42,水位M2波潮汐信號淹沒在強干擾信號中,計算的潮汐因子和相位滯后沒有參考價值,2018年下半年緩和期M2波潮觀測精度較好,但仍不及未打井前的1/10。

圖6 井水位M2波觀測振幅與相位隨時間變化Fig.6 Variation of observation amplitude and phase of M2 wave of well water level with time
從圖6可以看出,水位M2潮汐相位波動比振幅波動顯著,反映了觀測含水層花崗巖裂隙水非均質特性。打井抽水前M2潮汐振幅變化很小,打井后一直到2018年5月井水位M2潮汐振幅及相位均處于劇烈變化階段,總體上大幅上升、下降變化基本同步,部分時段出現不同步現象,初步分析認為打井及抽水引發了垂向流[13]。干擾趨緩期振幅及相位波動接近打井前水平,相位略有抬升、振幅有所下降,反映了打井后淺部地下水混入,原觀測含水系統的封閉性變差,潮汐應變響應能力有所下降。
采用多種時頻分析的方法,對仙游井水位打井抽水干擾前后的時頻特征變化和氣壓與固體潮效應特性變化的同步時頻分析,結果顯示:
(1) 打井/抽水干擾嚴重時段,井水位動態在時域上出現大幅階變和突跳變化,在頻域上亦影響顯著,抽水干擾信號完全掩蓋地下水位的固有周期信號及對氣壓、固體潮的響應效應。
(2) 打井/抽水干擾趨于緩和,中高頻段(8~23 h)基本修復,而低頻段(2~8 h)仍然存在多個強信號(幾倍于之前),顯示打井抽水干擾后,地下水含水系統的修復一個緩慢的過程。
(3) 井水位M2潮汐振幅及相位變化分析表明,打井/抽水引發井-含水系統的以徑向為主垂向為輔的排水,干擾趨緩后,井水位M2潮汐振幅及相位波動趨緩,井水位M2波觀測精度與打井前相比有所下降,潮汐應變響應能力略有下降。
綜上所述,打井抽水干擾對井水位的高、中、低頻段影響是不相同的,其影響具有一定的后效性,在水位資料異常識別與分析中應選取干擾較少的頻段來進行分析,對于仙游井水位來說中頻段是受影響最小、最有益于開展監測預報工作的頻段。