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

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

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

0 引言

有研究認(rèn)為，承壓井水位的變化反映了含水層孔隙壓力的變化，而含水層孔隙壓力的變化與含水層所受的壓力狀態(tài)有密切的關(guān)系。因此，觀測(cè)并記錄鉆井內(nèi)水位的變化，是測(cè)量地下深部彈性應(yīng)力變化的方法之一(黃輔瓊等，2004)。為了尋求承壓井水位與含水層應(yīng)力應(yīng)變的關(guān)系，國(guó)內(nèi)外學(xué)者開(kāi)展了大量工作。Narasimhan等(1984)給出了體應(yīng)變與固體潮引起的水位變化之間的定量關(guān)系式。黃輔瓊等(2004)利用華北地區(qū)40多口深井水位動(dòng)態(tài)觀測(cè)資料，定性地分析了大華北地區(qū)的構(gòu)造應(yīng)力場(chǎ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)力場(chǎng)變化圖像。

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

1 研究方法

研究表明，兩個(gè)因素可引起承壓井水頭的變化:承壓井含水層內(nèi)水量的變化和含水層所受的應(yīng)力應(yīng)變的變化。對(duì)于封閉性較差的含水層主要是第一種原因，例如降雨或同層抽水引起含水層水量發(fā)生變化;對(duì)于封閉性良好的含水層主要是第二種原因。對(duì)于完全封閉的含水層，若只有一口井，其水頭的變化只取決于含水層所受的應(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)可以看出，對(duì)于封閉性較好的深層承壓井來(lái)說(shuō)，通過(guò)計(jì)算井孔水位的潮汐因子和水位變化量，就可以求得該井的含水層體應(yīng)變的變化量。

2 資料的選取

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

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

3 潮汐因子的計(jì)算

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

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

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

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

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

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

3.3 去除長(zhǎng)周期趨勢(shì)和濾波

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

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

圖3 2013年1月1日～31日無(wú)極井水位

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

3.4 計(jì)算潮汐因子

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

3.5 剔除海潮負(fù)荷

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

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

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

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

4 水位趨勢(shì)變化的提取

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

圖4 不同趨勢(shì)類型水位變化曲線

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

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

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

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

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

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

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

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

圖8是反演得到的華北地區(qū)2009～2012年應(yīng)力應(yīng)變場(chǎng)變化圖像。由圖8可見(jiàn)，變化比較明顯的區(qū)域主要集中在山西中南部和首都圈地區(qū)。首都圈地區(qū)在2009年下半年和2012年都出現(xiàn)了應(yīng)力應(yīng)變集中的情況，這種集中主要表現(xiàn)為一種應(yīng)力應(yīng)變?cè)鰪?qiáng)與減弱相交的四象限分布特征，與區(qū)域應(yīng)力場(chǎng)背景相符。2010年首都圈地區(qū)發(fā)生了2次4級(jí)地震，分別是3月6日河北灤縣MS4.2地震和4月9日河北豐南MS4.1地震，由圖8可以看到，2009年下半年首都圈地區(qū)出現(xiàn)了一定程度的應(yīng)變集中，到2010年上半年該集中有所減弱，這與2次地震的發(fā)震時(shí)間相符。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ū)計(jì)算的體應(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年其他幾次中等地震也均是以正斷或走滑為主的地震活動(dòng)。2011年以后，山西地區(qū)地震活動(dòng)較為平靜，沒(méi)有顯著的地震事件發(fā)生。

6 結(jié)論與討論

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

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

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