利用華北地區(qū)承壓井水位資料反演含水層體應(yīng)變

楊柳 馬建英 曹井泉 邵永新 劉文兵2，

1)天津市地震局，天津市河西區(qū)友誼路19號 300201

2)中國地震局地球物理研究所，北京100081

0 引言

有研究認(rèn)為，承壓井水位的變化反映了含水層孔隙壓力的變化，而含水層孔隙壓力的變化與含水層所受的壓力狀態(tài)有密切的關(guān)系。因此，觀測并記錄鉆井內(nèi)水位的變化，是測量地下深部彈性應(yīng)力變化的方法之一(黃輔瓊等，2004)。為了尋求承壓井水位與含水層應(yīng)力應(yīng)變的關(guān)系，國內(nèi)外學(xué)者開展了大量工作。Narasimhan等(1984)給出了體應(yīng)變與固體潮引起的水位變化之間的定量關(guān)系式。黃輔瓊等(2004)利用華北地區(qū)40多口深井水位動態(tài)觀測資料，定性地分析了大華北地區(qū)的構(gòu)造應(yīng)力場狀態(tài)。張昭棟等(1999)利用含水層參數(shù)、固體潮效應(yīng)和氣壓效應(yīng)等3種方法反演含水層應(yīng)力應(yīng)變，并進(jìn)行了分析比較。孫小龍等(2011)利用華北地區(qū)63口井的水位變化和井含水層水文地質(zhì)參數(shù)資料，反演出華北地區(qū)構(gòu)造應(yīng)力場變化圖像。

由于利用含水層的水文地質(zhì)參數(shù)和水位氣壓效應(yīng)來反演含水層應(yīng)力應(yīng)變，是建立在水平層狀承壓含水層這種平面假設(shè)的基礎(chǔ)上的，并且在選取水文地質(zhì)參數(shù)時一般取其可能值或平均值，容易使計算結(jié)果出現(xiàn)誤差(張昭棟等，1999)。因此，在借鑒前人研究成果的基礎(chǔ)上，本文篩選了華北地區(qū)反應(yīng)固體潮較好的承壓井水位資料，利用井水位的固體潮效應(yīng)來反演華北地區(qū)現(xiàn)有承壓井孔的含水層體應(yīng)變。該方法不但結(jié)合了含水層的實(shí)際情況而且反映了體積力產(chǎn)生的應(yīng)力應(yīng)變，計算結(jié)果更為可靠。

1 研究方法

研究表明，兩個因素可引起承壓井水頭的變化:承壓井含水層內(nèi)水量的變化和含水層所受的應(yīng)力應(yīng)變的變化。對于封閉性較差的含水層主要是第一種原因，例如降雨或同層抽水引起含水層水量發(fā)生變化;對于封閉性良好的含水層主要是第二種原因。對于完全封閉的含水層，若只有一口井，其水頭的變化只取決于含水層所受的應(yīng)力應(yīng)變的變化。實(shí)際上由于含水層的封閉性不可能是完全理想的，所以水井含水層壓力水頭的變化基本上都是兩種原因的綜合(張昭棟等，1999)。

此處，本文僅考慮第二種理想的封閉性良好的含水層。根據(jù)孔隙彈性介質(zhì)理論，作用于含水層某一平面上的荷載分別由固體顆粒與顆粒間的孔隙流體共同承擔(dān)。由此，可以得到井水位變化與體應(yīng)變之間的關(guān)系式(Rhoads et al，1979)

其中，ρ為含水層內(nèi)水的密度;g為地球表面的重力加速度;n為含水層的孔隙度;Em和Ew分別為巖石固體顆粒和孔隙流體的體積模量;dh為含水層井水位的變化量;ΔΘ為體應(yīng)變量;δ為潮汐因子，表示單位體應(yīng)變引起的井水位變化。

由式(1)可以看出，對于封閉性較好的深層承壓井來說，通過計算井孔水位的潮汐因子和水位變化量，就可以求得該井的含水層體應(yīng)變的變化量。

2 資料的選取

“九五”、“十五”期間數(shù)字觀測技術(shù)在全國前兆觀測網(wǎng)絡(luò)中得到了廣泛的應(yīng)用，華北地區(qū)的井水位觀測已經(jīng)由模擬觀測逐步改造為以數(shù)字觀測為主的觀測網(wǎng)絡(luò)。數(shù)字化改造使得前兆資料的觀測、傳輸高度自動化，觀測數(shù)據(jù)的采樣密度和靈敏度也大幅提高(田山等，2009)。目前，井水位數(shù)字化觀測已全部實(shí)現(xiàn)了分鐘計值采樣，部分測項的采樣密度達(dá)到以秒計值，這為提取水位固體潮的潮汐因子提供了數(shù)據(jù)保障。

本文收集了2003年以來華北地區(qū)(34°～42.5°N，111°～123°E)反應(yīng)固體潮的承壓井的基礎(chǔ)資料和水位、氣壓的整點(diǎn)值數(shù)據(jù)，并篩選出數(shù)據(jù)質(zhì)量較好的33口井孔的水位資料。觀測井主要分布在河北、北京、天津、遼寧、內(nèi)蒙、山東、河南和山西等8個地區(qū)(圖1)。

3 潮汐因子的計算

計算潮汐因子，首先要提取研究井孔水位的固體潮信息。水位波動的影響因素主要有固體潮、地震波、氣壓、降雨、海潮等，對于封閉的含水層，這些影響因素是通過對含水層的作用導(dǎo)致其壓縮或膨脹，從而引起水頭的變化(史浙明等，2012)。由于水位固體潮的變化周期約為12～25h，而地震波以及降雨荷載與固體潮的周期相差較大，且海潮負(fù)荷僅對距離海洋較近區(qū)域的井水位產(chǎn)生明顯的影響。因此，利用地下水水位的固體潮反演含水層體應(yīng)變時，我們只排除氣壓的影響，并對受海潮影響較大區(qū)域的井孔排除了海潮負(fù)荷，從而得到33口承壓水位井的潮汐因子。

圖1 華北地區(qū)33口數(shù)字化水位井分布

3.1 數(shù)據(jù)的預(yù)處理

首先對觀測資料進(jìn)行預(yù)處理，消除因觀測條件改變或其他干擾造成的數(shù)據(jù)突跳、階變等情況，然后應(yīng)用3次樣條插值法將觀測儀器故障等原因造成的缺數(shù)補(bǔ)齊。圖2(b)為河北無極井預(yù)處理后的整點(diǎn)值水位數(shù)據(jù)曲線。

3.2 去除氣壓效應(yīng)

在各種影響因素中，氣壓對地下水位觀測的影響較為關(guān)鍵(萬永芳等，2009;劉學(xué)領(lǐng)等，2010)。目前，在去除地下水位觀測氣壓效應(yīng)的分析研究中，較為常用的方法有線性回歸、多元回歸等?；貧w方法的優(yōu)點(diǎn)是方法成熟、易于計算，但也存在一些不足。比如在計算過程中，主要是通過相關(guān)分析求取1個氣壓系數(shù)，進(jìn)而去除氣壓的影響，并沒有考慮氣壓隨時間的變化對水位的影響;此外，井水位的氣壓效應(yīng)普遍存在時間滯后現(xiàn)象(趙丹等，2013)，而上述回歸法并沒有考慮這種影響。應(yīng)用卷積回歸法來去除氣壓效應(yīng)就可解決上述問題，該方法是通過調(diào)整最大響應(yīng)時間以生成配對的水位和氣壓觀測值的最佳階躍響應(yīng)函數(shù)(氣壓響應(yīng)函數(shù))，再由氣壓響應(yīng)函數(shù)計算出修正后的水位變化。BETCO程序(Toll et al，2007)應(yīng)用卷積回歸法實(shí)現(xiàn)了估計氣壓和井水位變化之間的時間滯后效應(yīng)，并可以自動去除氣壓對觀測水位的影響(圖2)。

3.3 去除長周期趨勢和濾波

由于所要提取的固體潮周期為12～25h，因此為了得到更清晰的這種短周期的水位變化曲線，首先要去除水位的長周期趨勢變化。本文采用按月分段的一般多項式擬合方法來去除水位的趨勢項，圖3(a)為去除趨勢后的水位變化曲線。

圖2 2013年1月1～31日無極井氣壓和水位

圖3 2013年1月1日～31日無極井水位

去除水位長周期趨勢變化后，水位曲線中仍然存在較短周期的波動變化，因此還需要進(jìn)行濾波處理。對水位整點(diǎn)值做48h的滑動均值并濾波，無極井水位濾波后的曲線(圖3(c))與理論曲線(圖3(d))相比潮汐的變化形態(tài)已經(jīng)比較一致。

3.4 計算潮汐因子

井孔水位的固體潮可以分離成363個分波。調(diào)和分析結(jié)果表明，在這些分波中最主要的分波為半日波M2波，在水位固體潮振幅中，該波所占的比例達(dá)30% ～40%;其次是K1、S2、O1和N2等波(張國民等，2001)。因此，本文采用Venedikov調(diào)和分析方法，利用2008年1月～2012年12月的水位數(shù)據(jù)計算了33口水位井的M2波潮汐因子。

3.5 剔除海潮負(fù)荷

由于海潮負(fù)荷的體應(yīng)變響應(yīng)與距離有密切關(guān)系，觀測點(diǎn)距海岸線30km以內(nèi)，體應(yīng)變響應(yīng)隨著距離的減小而迅速增大(曹井泉等，2010);在30～90km范圍內(nèi)，隨著距離的增大，體應(yīng)變響應(yīng)逐漸減小;超過100km時，體應(yīng)變響應(yīng)明顯減弱。因此，本文選取了距離海岸線100km以內(nèi)的水位井做剔除海潮負(fù)荷的計算，包括寶坻、張道口、靜海、玉田、沈家臺、魯02、魯07等井。

本文應(yīng)用Schwiderski全球海潮模型，采用積分格林函數(shù)方法計算上述7口承壓井水位M2潮汐波的應(yīng)變負(fù)荷潮，在計算水位M2潮汐波的振幅和相位的基礎(chǔ)上，進(jìn)行海潮負(fù)荷改正。首先，以38.70°N、118.40°E點(diǎn)為海潮負(fù)荷點(diǎn)，根據(jù)李艷蕓等(2006)利用Coherens模型給出的渤海海潮圖，計算出7口水位觀測井的M2波海潮負(fù)荷矢量振幅和相位滯后。然后結(jié)合之前的水位調(diào)和分析結(jié)果，計算出觀測殘差矢量和剩余殘差矢量，最終求出經(jīng)海潮負(fù)荷改正后的潮汐因子和相位差，渤海海潮圖及具體的計算步驟和公式見李艷蕓等(2006)及曹井泉等(2010)文獻(xiàn)。

在此選取了2012年1月1日～3月31日共3個月的數(shù)據(jù)進(jìn)行海潮負(fù)荷改正的計算，表1顯示了海潮負(fù)荷改正前、后所得到的水位M2波潮汐因子和相位差的計算結(jié)果。由表1可見，距離海岸線50km以內(nèi)的魯02井和張道口井，海潮負(fù)荷改正前、后的M2潮汐因子計算結(jié)果相差較大;其他井孔的計算結(jié)果隨著距離的增大，差距在逐漸變小。而且經(jīng)過海潮改正，這幾口井的相位差均表現(xiàn)為負(fù)值且數(shù)值離散度較小，說明計算取得了較好的海潮負(fù)荷改正效果。

表1 水位M2潮汐波海潮負(fù)荷改正結(jié)果

4 水位趨勢變化的提取

分析33口井的原始觀測曲線圖，水位的變化可以分為趨勢性上升、趨勢性下降和起伏波動等3種類型，有的井水位還會在趨勢性變化上疊加年周期變化(圖4)。在利用井水位反演含水層體應(yīng)變時，水位下降表明含水層壓性減弱，張性增強(qiáng);水位上升表明含水層張性減弱，壓性增強(qiáng)。因此，承壓井水位的變化決定了井所在區(qū)域應(yīng)力場的變化。經(jīng)過統(tǒng)計，本文研究的水位井中呈現(xiàn)下降變化的有24口，占總數(shù)的73%。如果這種變化完全是由于含水層應(yīng)力狀態(tài)引起的，那么說明整個研究區(qū)域處于長期拉張狀態(tài)，這與研究區(qū)域的應(yīng)力應(yīng)變背景并不相符。進(jìn)一步分析這24口井的呈下降趨勢的水位資料，其中16口井的水位呈線性下降特征，且幅度較大，這些井主要集中在北京、河北、天津、山西和魯豫交界等地區(qū)。本文從“中國水文信息網(wǎng)”①水利部水文局，2003～2010，我國北方平原區(qū)地下水通報(中國水文信息網(wǎng))搜集了這幾個區(qū)域的年降雨量數(shù)據(jù)(圖5)，由圖5可見，年降雨量曲線并沒有持續(xù)下降，反而是呈現(xiàn)一種起伏性上升，因此井孔水位的持續(xù)下降并不是由降雨量變化引起的。

圖4 不同趨勢類型水位變化曲線

圖5 2006～2012年不同區(qū)域年降雨量

華北平原地處我國北方半干旱地區(qū)，隨著經(jīng)濟(jì)建設(shè)的不斷發(fā)展，深層承壓水被大量超采，引發(fā)了地下水位降落漏斗、地面沉降和咸水界面下移、地面塌陷等一系列地質(zhì)環(huán)境問題(楊麗芝等，2013)。楊明波等(2009)曾對北京地區(qū)地下水位動態(tài)特征的成因作了詳細(xì)分析，認(rèn)為影響地下水位的主要因素是地下水開采和降雨，地下水位的趨勢性變化則與區(qū)域地下水超量開采有直接關(guān)系，而具有明顯長期下降趨勢的北京昌平井正是處于嚴(yán)重超采的區(qū)域。因此，進(jìn)一步搜集了華北地區(qū)不同區(qū)域的年度地下水單位面積蓄變量①，其中北京、河北、山西等地的累計年蓄變量呈明顯的趨勢性下降(圖6)，與該區(qū)域內(nèi)水位井的趨勢性下降形態(tài)相一致;天津地區(qū)的沉降漏斗區(qū)基本覆蓋南部平原區(qū)，整個區(qū)域地下水均處于超采狀態(tài)(董克剛，2010);而魯豫交界地區(qū)深井水位的下降主要受周邊地?zé)衢_采的影響(孫小龍等，2013)。

圖6 2003～2010年北京、河北、山西地下水單位面積年累計蓄變量

通過上述分析可以判斷，華北地區(qū)多年接近線性下降的井水位絕大部分都是由于外界干擾導(dǎo)致的，并不是區(qū)域應(yīng)力應(yīng)變的反映。由于具有這種變化的井孔數(shù)量較多，變化幅度也較大，因此有必要對這些數(shù)據(jù)進(jìn)行處理，否則會影響我們對含水層應(yīng)力應(yīng)變情況的分析，故本文采用K-L法最佳直線擬合去除了這16口井水位的線性變化趨勢(圖7)。

圖7 寶坻井水位及處理后的趨勢性變化

上述處理過后，由于水位數(shù)據(jù)中除了含有長趨勢的變化，還包含著年變動態(tài)、月變和日變以及其他頻率的干擾波動。因此若要提取較長時間內(nèi)比如年度或半年尺度上的水位變化，還需要從水位的原始曲線中提取出長趨勢變化動態(tài)。本文采用小波趨勢分析方法提取33口井水位的趨勢變化，根據(jù)前人研究成果(萬永芳等，2009)并通過計算對比，選取了db4小波基的5階或6階分解結(jié)果(圖7(d))。

5 反演華北地區(qū)井含水層體應(yīng)變

反演華北地區(qū)井孔含水層體應(yīng)變，首先利用第4節(jié)中小波方法提取的井水位趨勢變化數(shù)據(jù)，計算出各井每半年的水位平均值，將前后兩個水位半年平均值逐步相減，得到每半年各井水位的變化值dh;然后將已經(jīng)計算好的各井潮汐因子值與水位變化量代入式(1)進(jìn)行計算，求得各井所反映的含水層體應(yīng)變量ΔΘ，最后運(yùn)用徑向基插值法，得到華北地區(qū)半年尺度的應(yīng)力應(yīng)變場的變化圖像(圖8)。

圖8是反演得到的華北地區(qū)2009～2012年應(yīng)力應(yīng)變場變化圖像。由圖8可見，變化比較明顯的區(qū)域主要集中在山西中南部和首都圈地區(qū)。首都圈地區(qū)在2009年下半年和2012年都出現(xiàn)了應(yīng)力應(yīng)變集中的情況，這種集中主要表現(xiàn)為一種應(yīng)力應(yīng)變增強(qiáng)與減弱相交的四象限分布特征，與區(qū)域應(yīng)力場背景相符。2010年首都圈地區(qū)發(fā)生了2次4級地震，分別是3月6日河北灤縣MS4.2地震和4月9日河北豐南MS4.1地震，由圖8可以看到，2009年下半年首都圈地區(qū)出現(xiàn)了一定程度的應(yīng)變集中，到2010年上半年該集中有所減弱，這與2次地震的發(fā)震時間相符。2012年的2幅圖在首都圈地區(qū)都清晰地出現(xiàn)了應(yīng)變集中，幅度比2009年的更大，而2012年首都圈地區(qū)相繼發(fā)生了5月28日河北唐山MS4.8地震、6月18日天津?qū)氎鍹S4.0地震和8月26日天津?qū)氎鍹S3.5地震。其中最為顯著的唐山MS4.8地震的震源機(jī)制解為右旋走滑類型(張躍剛等，2013)，該結(jié)果與圖8中顯示的首都圈地區(qū)計算的體應(yīng)變拉張壓縮分布形態(tài)比較相似。

圖8顯示，山西中南部地區(qū)2009和2010年呈現(xiàn)一種壓性減弱(拉性增強(qiáng))的狀態(tài)，2012年下半年變?yōu)閴盒栽鰪?qiáng)(拉性減弱)狀態(tài)。2009～2010年山西地區(qū)發(fā)生的最大地震是1月24日山西河津MS4.8地震，該地震的震源機(jī)制解類型是以正斷為主兼走滑(宋美琴等，2012)，而2009～2010年其他幾次中等地震也均是以正斷或走滑為主的地震活動。2011年以后，山西地區(qū)地震活動較為平靜，沒有顯著的地震事件發(fā)生。

6 結(jié)論與討論

本文選取華北地區(qū)33口水位井的整點(diǎn)值數(shù)據(jù)，經(jīng)過去除氣壓效應(yīng)、長周期趨勢和濾波處理提取了水位固體潮信息，計算出井孔水位的潮汐因子，并采用積分格林函數(shù)方法剔除了距海岸線較近水位井的潮汐因子的海潮負(fù)荷;同時，通過對水位趨勢性變化的分析，粗略地去除了部分受外界影響較大的水位井的線性趨勢;最后以半年為時間尺度對選取的華北地區(qū)井水位進(jìn)行了反演計算，得到了以半年為尺度的應(yīng)力應(yīng)變場的動態(tài)變化圖像。

由計算結(jié)果可見，井孔分布較為密集的首都圈地區(qū)，在顯著地震前呈現(xiàn)出一種拉張、壓縮相交的近似四象限分布的狀態(tài);而對于井孔分布較少的山西帶則基本呈現(xiàn)出單純的拉張或壓縮的應(yīng)力集中狀態(tài)，這種情況與當(dāng)時的區(qū)域應(yīng)力狀態(tài)有一定關(guān)系，但也與本文選取的井孔分布不均勻有關(guān)。張昭棟等(2001)曾利用井水位資料反演了大同-陽高6.1級地震前后大區(qū)域應(yīng)力場的動態(tài)變化圖像(38°～42°N，114°～121°E)，在地震前后該區(qū)域也呈現(xiàn)一種“正應(yīng)力”與“負(fù)應(yīng)力”交織的狀態(tài)，與2012年唐山4.8級地震前后的空間圖像較為相似(圖8)。曹井泉等(2004)利用應(yīng)力變化響應(yīng)系數(shù)反演了1976年唐山7.8級地震前后3年京津唐地區(qū)的應(yīng)力場變化情況，在該次地震前震中區(qū)和周邊地區(qū)出現(xiàn)了應(yīng)力應(yīng)變集中的有序變化。結(jié)合上述震例可以看出，本文對華北地區(qū)井水位含水層體應(yīng)變的反演在一定程度上反映出了該區(qū)域深層的應(yīng)力應(yīng)變狀態(tài)。

在本文的研究時段內(nèi)，華北地區(qū)的許多井孔由于受到區(qū)域地下水超采等情況的影響而出現(xiàn)了趨勢性下降，本文采用K-L直線擬合的方法去除了這一影響。由于這種方法比較簡單，并不能完全模擬環(huán)境干擾對井水位的影響，因此還需要進(jìn)步一研究更貼近實(shí)際的干擾去除方法。