張磊 邵永新 馬建英 李恩建 金大利 王靜

1）中國(guó)礦業(yè)大學(xué)（北京），北京市海淀區(qū)學(xué)院路丁11號(hào) 100083

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

3）山東省地震局，山東濟(jì)南 250014

0 引言

伸縮儀是精密測(cè)量地殼巖體兩點(diǎn)間水平距離相對(duì)變化的儀器，適用于觀測(cè)地殼應(yīng)變和固體潮水平分量的連續(xù)變化，為研究地震孕育過(guò)程中水平應(yīng)變的變化規(guī)律提供數(shù)據(jù)，也為地球彈性研究提供重要數(shù)據(jù)（熊仲華，2006）。抽水開(kāi)采是洞體應(yīng)變觀測(cè)中的一種重要干擾源（金克儉等，1986），也是伸縮應(yīng)變觀測(cè)較為常見(jiàn)的干擾因素（李國(guó)斌等，2004；李國(guó)斌，2006；汪翠枝等，2009）。已有研究表明（裴曉峰，1995），薊縣地震臺(tái)定點(diǎn)形變觀測(cè)所受干擾程度最大、發(fā)生最為頻繁的是機(jī)井抽水的干擾。SSY型伸縮儀安裝于薊縣臺(tái)東約4.5km處的小辛莊山洞內(nèi)，距洞口220～260m處，自2001年5月開(kāi)始正式觀測(cè)，基線(xiàn)長(zhǎng)39.7m，NS測(cè)向方位角144.87°，EW測(cè)向方位角59.94°。小辛莊山洞周邊存在多口用于日常飲水和農(nóng)田灌溉的抽水井，水井距觀測(cè)山洞最近約300m，井深一般在100～200m，多分布在山洞的南北兩側(cè)，伸縮應(yīng)變受抽水干擾嚴(yán)重，表現(xiàn)為相對(duì)拉張變化，尤其以NS向更為敏感。

薊縣地震臺(tái)位于首都圈中東部，“燕山沉降帶”的“京東臺(tái)凹”的北部，是全國(guó)定點(diǎn)地殼形變基本臺(tái)。臺(tái)站附近斷裂構(gòu)造較為復(fù)雜，存在一系列斷裂和褶皺，臺(tái)站附近有薊縣山前斷裂、楊莊斷裂和黃崖關(guān)斷裂。小辛莊山洞山體為震旦系霧迷山組灰?guī)r（含燧石），巖層層面走向335°，傾角20°～33°，測(cè)洞南距薊縣斷裂1.5km。儀器室內(nèi)的日溫差小于0.03℃，年溫差小于0.05℃，相對(duì)濕度達(dá)98%。桑梓井位于小辛莊山洞西南約37km處的薊縣桑梓，自2004年開(kāi)始正式觀測(cè)，井深300m左右，含水層屬承壓水，其北部為薊縣山前斷裂（圖1）。

圖1 天津地傾斜臺(tái)站與斷裂分布

為了研究薊縣臺(tái)小辛莊山洞伸縮應(yīng)變受周邊短時(shí)地下水抽取的影響，本文利用水位的變化獲取與地下水抽取有關(guān)的信息，但是小辛莊山洞周?chē)鷽](méi)有水位觀測(cè)井，故文中利用天津地震前兆臺(tái)網(wǎng)中薊縣桑梓井水位數(shù)據(jù)，通過(guò)分析小辛莊伸縮應(yīng)變與桑梓井水位的年動(dòng)態(tài)變化特征，判斷二者對(duì)短時(shí)抽取地下水的響應(yīng)是否存在一致性，進(jìn)而建立水位變化量與應(yīng)變變化量的關(guān)系，研究伸縮應(yīng)變受周邊短時(shí)抽水干擾變化幅度，以期對(duì)在季節(jié)性抽水時(shí)段伸縮應(yīng)變出現(xiàn)短期加速拉張變化的異常判定提供幫助。

1 年動(dòng)態(tài)特征分析

本文整理了薊縣臺(tái)小辛莊山洞伸縮應(yīng)變和薊縣臺(tái)桑梓井水位2004年1月1日～2014年6月30日以來(lái)日均值數(shù)據(jù)（圖2），由圖可知，二者均表現(xiàn)出了較好的年動(dòng)態(tài)變化特征：

（1）小辛莊伸縮應(yīng)變自2004年以來(lái)NS、EW兩方向年變形態(tài)呈現(xiàn)出一定的規(guī)律性變化。NS向每年年初變化基本平穩(wěn)，4月前后受季節(jié)性抽水的影響呈現(xiàn)拉張變化，并在此基礎(chǔ)上波動(dòng)變化，6～8月份雨季到來(lái)受降雨影響呈現(xiàn)壓縮變化，此后逐漸恢復(fù)呈相對(duì)拉張變化；EW向每年年初整體上呈現(xiàn)壓縮變化，4月前后受地下水抽取影響呈現(xiàn)相對(duì)拉張變化，之后受降雨影響呈現(xiàn)壓縮變化，此后逐漸恢復(fù)呈現(xiàn)相對(duì)拉張變化。

（2）桑梓水位自2004年開(kāi)始觀測(cè)以來(lái)具有較好的年周期變化特征。每年3月前后為水位的高值期，之后轉(zhuǎn)降；6月前后降到年度最低值，下降幅度為5m左右，之后轉(zhuǎn)升；到11月前后達(dá)到年度第2次高水位狀態(tài)，上升幅度為3m左右，之后再次下降；到12月初達(dá)到年度第2次低水位期，下降幅度為1m左右，之后又開(kāi)始上升，一直到轉(zhuǎn)年的高水位期。

圖2 薊縣臺(tái)桑梓井水位和小辛莊伸縮應(yīng)變NS向、EW向日均值曲線(xiàn)

圖3表明，薊縣臺(tái)小辛莊伸縮應(yīng)變觀測(cè)對(duì)抽水響應(yīng)敏感，尤其是NS向更為突出。受季節(jié)性抽水影響，薊縣臺(tái)小辛莊伸縮應(yīng)變?cè)?月前后加速拉張變化顯著（圖4），水位在3月前后持續(xù)高值平穩(wěn)變化后轉(zhuǎn)降。對(duì)比2004年以來(lái)歷年年初桑梓井水位與小辛莊伸縮應(yīng)變數(shù)據(jù)（圖4）發(fā)現(xiàn)，伸縮兩分向加速拉張與水位加速下降在時(shí)間上具有很好的同步性。雖然薊縣臺(tái)小辛莊山洞與桑梓井空間上還存在一定的距離，但二者對(duì)季節(jié)性抽水的響應(yīng)表現(xiàn)出很好的同步性變化。楊耀棟等（2011）發(fā)現(xiàn)，天津平原區(qū)深層地下水受開(kāi)采影響表現(xiàn)為徑流-越流-開(kāi)采型動(dòng)態(tài)特征，因受灌溉強(qiáng)開(kāi)采影響，低水位多出現(xiàn)于5～6月份，有時(shí)出現(xiàn)在9～10月份；豐水期8月停采后，水位逐漸回升，有時(shí)出現(xiàn)高水位，大多至翌年1～3月水位最高，高水位期較最大降水期滯后5～7個(gè)月，年水位變幅4～6m。桑梓水位自2004年觀測(cè)以來(lái)的年變化形態(tài)與上述區(qū)域深層地下水動(dòng)態(tài)變化特征一致。另外，陳鵬等（2004）在薊縣地震臺(tái)開(kāi)展抽水實(shí)驗(yàn)研究發(fā)現(xiàn)，小辛莊伸縮應(yīng)變觀測(cè)受抽水干擾的幅度既與抽水井到形變監(jiān)測(cè)儀器距離有關(guān)，也與抽水井水位的最大降深有關(guān)。

2 回歸分析

薊縣臺(tái)小辛莊伸縮應(yīng)變加速拉張變化與桑梓井水位加速下降變化同步性良好，本文利用最小二乘線(xiàn)性回歸分析方法研究二者之間的關(guān)系，建立伸縮應(yīng)變變化量與水位變化量之間的線(xiàn)性回歸模型，并利用t檢驗(yàn)方法對(duì)模型進(jìn)行顯著性檢驗(yàn)，最后將2014年水位實(shí)測(cè)變化量代入回歸模型計(jì)算應(yīng)變量的理論變化區(qū)間，驗(yàn)證2014年薊縣臺(tái)小辛莊伸縮應(yīng)變加速拉張變化與地下水抽取的關(guān)系。

2.1 回歸分析方法（李博納等，2006）

對(duì)于給定數(shù)據(jù)點(diǎn) （xi，yi），i＝1，2…n，滿(mǎn)足如下關(guān)系

圖3 2012年3月12～13日薊縣臺(tái)小辛莊伸縮應(yīng)變NS、EW向受抽水干擾數(shù)據(jù)變化

圖4 薊縣臺(tái)小辛莊伸縮應(yīng)變NS向、EW向與多年桑梓井水位日均值曲線(xiàn)的對(duì)比

其最小二乘回歸直線(xiàn)可以表示為

方差σ2的無(wú)偏估計(jì)

相關(guān)系數(shù)為

利用t檢驗(yàn)方法進(jìn)行回歸模型顯著性檢驗(yàn)。假設(shè)：H0：β1＝0，H1：β1≠ 0，當(dāng)H0成立時(shí)

對(duì)于給定顯著性水平α，記t為統(tǒng)計(jì)量（7）的觀測(cè)值，若

則拒絕H0：β1＝0，即認(rèn)為自變量x對(duì)因變量y有顯著性影響。

給定顯著性水平α，在x＝x0處，標(biāo)準(zhǔn)差為

2.2 回歸分析結(jié)果

2.2.1 回歸模型的建立

由于小辛莊伸縮應(yīng)變兩方向加速拉張與桑梓井水位加速下降變化在時(shí)間上具有很好的同步性，本文統(tǒng)計(jì)了2004～2013年年初桑梓水位最高時(shí)水位和伸縮應(yīng)變NS向、EW向觀測(cè)值以及伸縮應(yīng)變兩方向加速拉張變化結(jié)束時(shí)的應(yīng)變最大值和水位值，計(jì)算水位最高至伸縮兩方向拉張結(jié)束時(shí)段內(nèi)水位和應(yīng)變觀測(cè)的變化量（表1、表2）。

根據(jù)表1、表2可以得到水位變化量與伸縮應(yīng)變變化量的對(duì)應(yīng)關(guān)系（表3）。表3中xNi、yNi分別為研究水位變化量與伸縮NS向應(yīng)變變化關(guān)系的樣本點(diǎn)，xNi為薊縣臺(tái)桑梓井水位變化量，yNi為薊縣臺(tái)小辛莊伸縮NS向應(yīng)變變化量；xEi、yEi為研究水位變化量與伸縮 EW向應(yīng)變變化關(guān)系的樣本點(diǎn)，xEi為薊縣臺(tái)桑梓井水位變化量，yEi為薊縣臺(tái)小辛莊伸縮EW向應(yīng)變變化量。xNi，yNi與xEi、yEi的分布如圖5所示，兩組數(shù)據(jù)的相關(guān)系數(shù)分別為 0.8342、0.5884，可見(jiàn)薊縣臺(tái)桑梓井水位變化量與小辛莊伸縮應(yīng)變變化量相關(guān)性良好。

表1 薊縣臺(tái)桑梓井水位最高值和小辛莊伸縮應(yīng)變NS向拉張最大值時(shí)水位與應(yīng)變的變化

表2 薊縣臺(tái)桑梓井水位最高值和小辛莊伸縮應(yīng)變EW向拉張最大值時(shí)水位與應(yīng)變的變化

表3 桑梓井水位變化量與小辛莊伸縮應(yīng)變變化量的對(duì)應(yīng)關(guān)系

由表3數(shù)據(jù)，利用2.1中的最小二乘線(xiàn)性回歸方法，根據(jù)式（2～4）計(jì)算得到最小二乘回歸直線(xiàn)相關(guān)參數(shù)，如表4所示。

表4 最小二乘線(xiàn)性回歸結(jié)果參數(shù)表

由表4可得薊縣臺(tái)桑梓井水位變化量與小辛莊伸縮NS向、EW向應(yīng)變變化量之間的最小二乘回歸直線(xiàn)（圖5）表示為

相應(yīng)的式（10）回歸直線(xiàn)方差的無(wú)偏估計(jì)分別為

圖5 薊縣臺(tái)小辛莊伸縮應(yīng)變變化量與桑梓井水位變化量線(xiàn)性擬合分析結(jié)果

2.2.2 回歸模型顯著性檢驗(yàn)

本文利用t檢驗(yàn)方法對(duì)回歸模型（10）進(jìn)行顯著性檢驗(yàn)，依據(jù)式（7～9）分別用表示回歸模型式（10）中直線(xiàn)xE的t統(tǒng)計(jì)量的觀測(cè)值。則有：

（1）對(duì)于薊縣臺(tái)桑梓井水位變化量與小辛莊伸縮NS向應(yīng)變變化量的回歸直線(xiàn)：44.6858＋328.7437xN，有n＝10，假設(shè)：H0：β1＝0；H1：β1≠0，對(duì)于給定顯著性水平α＝0.05，t0.025（n-2）＝t0.025（8）＝2.3060，t統(tǒng)計(jì)量的觀測(cè)值，則拒絕H0，認(rèn)為自變量xN對(duì)因變量yN有顯著性影響即水位變化對(duì)伸縮NS向應(yīng)變變化影響顯著，置信水平達(dá)到了95%。

（2）對(duì)于薊縣臺(tái)桑梓井水位變化量與小辛莊伸縮EW向應(yīng)變變化量的回歸直線(xiàn)：有n＝10，給定顯著性2.0582＜t0.025（8），不滿(mǎn)足顯著性水平為α＝0.05的假設(shè)檢驗(yàn)；于是取α＝0.1，有t0.05（n-2）＝1.8595，），則在給定顯著性水平α＝0.1條件下，自變量xE對(duì)因變量yN有顯著性影響即水位變化對(duì)伸縮EW向應(yīng)變變化影響顯著，置信水平達(dá)到了90%。

上述t檢驗(yàn)結(jié)果表明，水位變化對(duì)伸縮應(yīng)變變化影響顯著，并且水位變化量與應(yīng)變變化量的回歸直線(xiàn)方程N(yùn)S向置信水平高于EW向，這與實(shí)際觀測(cè)中伸縮應(yīng)變NS向?qū)Τ樗憫?yīng)更為敏感這一實(shí)測(cè)現(xiàn)象也是一致的。

2.2.3 應(yīng)變量變化區(qū)間分析

薊縣臺(tái)小辛莊伸縮應(yīng)變觀測(cè)受抽水干擾影響嚴(yán)重，一般在每年4月前后開(kāi)始加速拉張變化（圖6），但2014年伸縮應(yīng)變兩方向自3月中旬開(kāi)始加速拉張變化，時(shí)間較往年同期有所提前，且變化較為突出（圖7）。震例研究表明（汪翠枝等，2009），薊縣伸縮應(yīng)變震前一般表現(xiàn)為短臨前兆異常，NS向在文安地震前、EW向在盧龍地震前均出現(xiàn)短期加速拉張變化。針對(duì)這種異常變化，首先對(duì)比小辛莊伸縮應(yīng)變和桑梓井水位日均值曲線(xiàn)，發(fā)現(xiàn)水位加速轉(zhuǎn)降與伸縮兩方向加速拉張變化開(kāi)始的時(shí)間基本一致。另外，調(diào)查發(fā)現(xiàn)，此期間儀器工作正常，周邊也并未發(fā)現(xiàn)新的干擾源，故初步判定此異常變化可能與季節(jié)性抽水有關(guān)。

圖6 薊縣臺(tái)小辛莊伸縮應(yīng)變NS、EW向日均值多年對(duì)比曲線(xiàn)

本文依據(jù)表1、表2，統(tǒng)計(jì)了2014年水位最高時(shí)水位和伸縮NS向、EW向觀測(cè)值、伸縮應(yīng)變兩方向加速拉張變化結(jié)束時(shí)的應(yīng)變最大值和水位值，分別計(jì)算出2014年該時(shí)段內(nèi)薊縣臺(tái)桑梓井水位和小辛莊伸縮 NS向、EW向的實(shí)測(cè)應(yīng)變變化量（表5），分別用（xN0，yN0）和（xE0，yE0）表示，相應(yīng)的標(biāo)準(zhǔn)差分別用δN和δE表示。如表5所示，依據(jù)式（10）和式（9）有：

（1）當(dāng)薊縣臺(tái)桑梓井水位變化量xN0＝2.8841時(shí)，代入式得到2014年伸縮NS向應(yīng)變變化量的理論估計(jì)值標(biāo)準(zhǔn)差δN＝408.8522，2014年薊縣臺(tái)小辛莊伸縮應(yīng)變NS向的實(shí)測(cè)變化量yN0＝674.4375，落入?yún)^(qū)間即（583.9633，1401.6677）內(nèi)，置信水平為 95%；

圖7 2014年薊縣臺(tái)小辛莊伸縮應(yīng)變NS向、EW向與桑梓井水位日均值對(duì)比曲線(xiàn)

表5 2014年桑梓井水位最高值和小辛莊伸縮應(yīng)變拉張最大值時(shí)水位與應(yīng)變變化

（2）當(dāng)薊縣臺(tái)桑梓井水位變化量xE0＝2.1929時(shí)，代入式得到2014年伸縮EW向應(yīng)變變化量的理論估計(jì)值、標(biāo)準(zhǔn)差δE＝472.3435，2014年薊縣臺(tái)小辛莊伸縮EW向的實(shí)測(cè)應(yīng)變變化量yE0＝481.05，落入?yún)^(qū)間即（-93.7461，850.9409）內(nèi)，置信水平為 90%。

對(duì)于給定任意xN0和xE0，yN0和yE0的變化區(qū)間理論上下限曲線(xiàn)如圖5中虛線(xiàn)所示。2014年薊縣臺(tái)桑梓井水位的實(shí)測(cè)變化量和小辛莊伸縮NS向、EW向的實(shí)測(cè)應(yīng)變變化量（2.8841，674.4375）、（2.1929，481.05）在理論控制線(xiàn)內(nèi)。因此，伸縮應(yīng)變?cè)?014年3月中旬以來(lái)出現(xiàn)的加速拉張變化與季節(jié)性抽水有關(guān)，作為地震前兆異常的信度不高。

2.3 地下水抽取對(duì)小辛莊伸縮應(yīng)變觀測(cè)影響的機(jī)制分析

地下水抽取過(guò)程實(shí)質(zhì)上為以井孔為中心的含水層的疏干過(guò)程，由于含水層失水產(chǎn)生地表沉降，形成以井孔為中心的局部收縮區(qū)域，受其影響井孔中心收縮區(qū)以外受到拉張作用，這就解釋了小辛莊山洞附近短時(shí)抽水對(duì)伸縮應(yīng)變的拉張影響。Terzaghi于1925年提出有效應(yīng)力原理（王大純等，2002）可以解釋抽水引起地下水位動(dòng)態(tài)變化進(jìn)而引發(fā)地面沉降的問(wèn)題。有效應(yīng)力即為實(shí)際作用于砂層骨架上的應(yīng)力，等于總應(yīng)力減去孔隙水壓力。抽水造成地下水水位下降，孔隙水壓力降低，有效應(yīng)力增加，砂層顆粒發(fā)生位移，顆粒接觸面積增加，排列更趨緊密，孔隙度降低，砂層受到壓密，從而產(chǎn)生地面沉降，含水層的孔隙度、給水度、滲透系數(shù)等參數(shù)均將變小。砂礫類(lèi)巖土基本呈彈性形變，待砂層孔隙水壓恢復(fù)后，砂層能大體恢復(fù)原狀。

已有水文資料表明薊縣臺(tái)小辛莊山洞附近地下水類(lèi)型屬于潛水，利用地下水潛水完整井穩(wěn)定流Dupuit公式可以得到抽水井影響半徑，該區(qū)域旱季和雨季井水位降深對(duì)形變干擾的影響范圍在800～1000m（陳鵬等，2004；李恩建等，2006）。楊耀棟等（2011）的研究也表明，天津地區(qū)抽水開(kāi)采是影響深層地下水位動(dòng)態(tài)的主要因素，地下水開(kāi)采與全淡水區(qū)的淺層地下水位動(dòng)態(tài)具有良好的相關(guān)性。文中回歸分析結(jié)果也表明，桑梓水位的短時(shí)下降變化與伸縮應(yīng)變的短時(shí)加速拉張變化相關(guān)性良好，二者之間的回歸模型置信水平達(dá)到了90%，同時(shí)也說(shuō)明了利用桑梓水位作為區(qū)域地下水動(dòng)態(tài)背景信息研究伸縮應(yīng)變受抽水影響的短時(shí)變化是具有實(shí)際意義的。

3 結(jié)論與討論

薊縣小辛莊伸縮應(yīng)變觀測(cè)受抽水干擾影響顯著，每年年初受季節(jié)性抽水影響NS、EW兩方向會(huì)表現(xiàn)出不同程度的短期加速拉張變化。桑梓水位多年來(lái)年動(dòng)態(tài)特征良好且與區(qū)域深層地下水動(dòng)態(tài)特征一致，在季節(jié)性抽水時(shí)段，水位加速下降與伸縮應(yīng)變加速拉張變化時(shí)間上具有很好的同步性。本文利用天津地震前兆臺(tái)網(wǎng)的觀測(cè)數(shù)據(jù)開(kāi)展了地下水抽取與伸縮應(yīng)變短時(shí)變化關(guān)系的研究，得到如下結(jié)果：

（1）文中統(tǒng)計(jì)并計(jì)算了2004～2013年水位最高值至伸縮應(yīng)變拉張結(jié)束時(shí)水位變化量和應(yīng)變的變化量，利用最小二乘線(xiàn)性回歸方法得到了水位變化量與應(yīng)變變化量的回歸直線(xiàn)（式（10）），二者呈正相關(guān)，即伸縮應(yīng)變受抽水影響的程度隨著水位變化量的增大而增大，這與陳鵬等（2004）抽水實(shí)驗(yàn)的研究結(jié)果也是一致的。

（2）水位變化量與伸縮NS向和EW向應(yīng)變變化量的線(xiàn)性回歸模型分別滿(mǎn)足顯著性水平為α＝0.05和α＝0.1的t檢驗(yàn)，即回歸模型的置信水平分別為95%和90%。回歸模型置信水平較高，水位變化量與NS向應(yīng)變變化量的回歸模型置信水平更高，這與實(shí)際觀測(cè)過(guò)程中伸縮應(yīng)變NS向?qū)Τ樗憫?yīng)更為敏感的客觀事實(shí)一致。

（3）伸縮應(yīng)變兩方向在2014年3月中旬開(kāi)始加速拉張變化，考量年變形態(tài)發(fā)現(xiàn)，時(shí)間較往年有所提前，且變化較為突出。而調(diào)查中并未發(fā)現(xiàn)新的干擾源，且儀器工作正常。利用文中的回歸分析方法計(jì)算得到2014年伸縮兩方向應(yīng)變變化量的理論變化區(qū)間，并且2014年應(yīng)變的實(shí)測(cè)變化量在變化區(qū)間內(nèi)，證明了該變化為季節(jié)性抽水所致，作為地震前兆異常的信度較低。

桑梓水位反映了區(qū)域地下水的動(dòng)態(tài)背景，在季節(jié)性抽水時(shí)段桑梓水位變化量與小辛莊伸縮應(yīng)變變化量的回歸模型置信水平達(dá)到了90%以上，回歸模型顯示伸縮應(yīng)變變化受水位變化影響顯著。由于小辛莊山洞附近缺乏相應(yīng)的水位觀測(cè)，因此可以利用桑梓水位作為輔助資料研究伸縮應(yīng)變受季節(jié)性抽水短時(shí)干擾的變化幅度。隨著觀測(cè)資料的積累，樣本量逐漸增加，回歸模型可以不斷進(jìn)行修正并進(jìn)行實(shí)例驗(yàn)證，這有利于輔助判定在季節(jié)性抽水時(shí)段伸縮應(yīng)變出現(xiàn)短期加速拉張變化的異常信度。另外，在受抽水干擾影響顯著的形變觀測(cè)臺(tái)站，可以開(kāi)展水位觀測(cè)項(xiàng)目，以輔助研究抽水對(duì)地形變觀測(cè)的短時(shí)影響。