顧悅 裴爍瑾 梁姍姍 樊春燕 周元澤

1)中國(guó)科學(xué)院大學(xué)，地球與行星科學(xué)學(xué)院，北京 100049

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

0 引言

大型含水層系統(tǒng)中的地下水常被作為水源，是水循環(huán)中的重要組成部分(Kim et al，2019)。地下水質(zhì)量和水位變化監(jiān)測(cè)主要基于鉆孔等定點(diǎn)直接測(cè)量(張磊等，2021； 門北方，2018)，其成本相對(duì)較高，因而可靠的低成本間接手段越來(lái)越受到研究人員和生產(chǎn)部門的關(guān)注。目前，基于合成孔徑雷達(dá)(InSAR)、全球定位系統(tǒng)(GPS)或大地測(cè)量等方法觀測(cè)到的變形模式或重力異常調(diào)整，可以在區(qū)域尺度監(jiān)測(cè)地下水的賦存和水位變化等(Amelung et al，1999；Argus et al，2014；Kim et al，2019； 李嘉等，2020)，但是對(duì)于科學(xué)研究和生產(chǎn)實(shí)踐，區(qū)域性地下水水位變化監(jiān)測(cè)在空間尺度上仍然偏大。

地震波對(duì)于地下介質(zhì)結(jié)構(gòu)的變化比較敏感(Huang et al，2010；Mordret et al，2010；Yeh et al，2013)。降水、河流補(bǔ)給及地下水抽采等對(duì)區(qū)域地下水調(diào)整的影響可以在地震波速結(jié)構(gòu)變化中有較好的反映(Mainsant et al，2012)。固定地震臺(tái)站連續(xù)觀測(cè)記錄可以持續(xù)記錄非地震的地震動(dòng)，即背景噪聲。連續(xù)觀測(cè)的地震背景噪聲互相關(guān)函數(shù)可用于持續(xù)監(jiān)測(cè)地下速度結(jié)構(gòu)的變化(Shapiro et al，2004；Gerstoft et al，2008；Lin et al，2008；Zheng et al，2008；Hadziioannou et al，2011；Brenguier et al，2014； 溫?fù)P茂等，2019)。

地震背景噪聲測(cè)量已經(jīng)被引入到監(jiān)測(cè)地下水位調(diào)整引起的地下速度結(jié)構(gòu)變化中(Kim et al，2019)。目前，國(guó)內(nèi)應(yīng)用地震背景噪聲監(jiān)測(cè)水庫(kù)水位相關(guān)的波速變化情況取得了一定的成果。例如，陳蒙等(2013)運(yùn)用背景噪聲技術(shù)對(duì)云南省大銀甸水庫(kù)周圍波速變化進(jìn)行研究，發(fā)現(xiàn)水位變化與波速變化存在明顯的相關(guān)性，并認(rèn)為這種相關(guān)性是由水庫(kù)水體的卸載作用造成的； 安艷茹等(2015)研究了紫坪鋪水庫(kù)區(qū)域在蓄水與泄水過(guò)程中庫(kù)區(qū)介質(zhì)的變化特征，表明地下介質(zhì)的相對(duì)波速變化與水位變化存在較為明顯的相關(guān)性，且在時(shí)間上有一定延遲，認(rèn)為其可能與水的滲透作用有關(guān)。

本文基于相對(duì)圈閉的臨沂地區(qū)國(guó)家測(cè)震臺(tái)網(wǎng)JUX與JUN臺(tái)站記錄的連續(xù)波形數(shù)據(jù)，利用背景噪聲互相關(guān)技術(shù)研究地下介質(zhì)波速變化，并結(jié)合臨沂地區(qū)地下水水位和降水等觀測(cè)數(shù)據(jù)，分析降水和地下水抽取等諸多因素影響下的地下水位變化與地震波速變化的關(guān)系。

1 研究區(qū)域及資料選取

山東省臨沂地區(qū)位于魯中南低山丘陵區(qū)東南部和魯東丘陵南部，地勢(shì)西北高、東南低，自北向南有沂山、蒙山、尼山等3條主要山脈，呈NE-SW向延伸，發(fā)育沂河與沭河(圖1)。區(qū)域水系相對(duì)圈閉(李致家等，2005)，地下水位變化主要由于天然降水和工、農(nóng)業(yè)使用。

圖1 研究區(qū)域概況

本地區(qū)有國(guó)家數(shù)字測(cè)震臺(tái)網(wǎng)的多個(gè)寬頻帶地震臺(tái)、中國(guó)地震局地下水位動(dòng)態(tài)監(jiān)測(cè)臺(tái)莒南魯14號(hào)井以及中國(guó)氣象局的莒縣氣象觀測(cè)站等。本文收集并使用了如下3種類型的數(shù)據(jù)資料：

(1)地震波形資料為2009年1月1日—2011年1月1日山東省臨沂市莒南縣的JUN臺(tái)與日照市莒縣的JUX臺(tái)的連續(xù)記錄，2個(gè)地震臺(tái)站間距為34.5km。數(shù)據(jù)資料來(lái)自于國(guó)家數(shù)字測(cè)震臺(tái)網(wǎng)數(shù)據(jù)備份中心(鄭秀芬等，2009)，2個(gè)臺(tái)站數(shù)據(jù)在2009、2010年分別有19天和13天出現(xiàn)缺失或損壞，可用數(shù)據(jù)的天數(shù)總計(jì)698天。

(2)地下水水位連續(xù)數(shù)據(jù)來(lái)自莒南魯14號(hào)井，該井為研究區(qū)內(nèi)距離最近的與地震監(jiān)測(cè)相關(guān)的水井。

(3)日降水?dāng)?shù)據(jù)由莒縣氣象觀測(cè)站(編號(hào)54936)所記錄，數(shù)據(jù)來(lái)自國(guó)家氣象科學(xué)數(shù)據(jù)中心。

2 地震數(shù)據(jù)處理

本文使用NoisePy軟件(Jiang et al，2020)處理數(shù)據(jù)，并計(jì)算相對(duì)速度變化。

2.1 連續(xù)記錄數(shù)據(jù)預(yù)處理

數(shù)據(jù)預(yù)處理(Bensen et al，2007)主要包括均值和趨勢(shì)變化去除、波形尖滅、儀器響應(yīng)去除、帶通濾波、時(shí)間域歸一化及譜白化處理等。具體如下：

(1)為了使數(shù)據(jù)標(biāo)準(zhǔn)化，進(jìn)行除均值來(lái)去除波形數(shù)據(jù)本身具有的非零均值。通常波形數(shù)據(jù)會(huì)存在一個(gè)長(zhǎng)周期的線性趨勢(shì)，從而影響數(shù)據(jù)的分析，故需進(jìn)行去趨勢(shì)處理。

(2)在對(duì)數(shù)據(jù)進(jìn)行譜域操作(如FFT、濾波等)時(shí)，若數(shù)據(jù)的兩端不為零，則會(huì)出現(xiàn)譜域假象，因而實(shí)際數(shù)據(jù)經(jīng)常需要做尖滅處理，使得數(shù)據(jù)兩端在短時(shí)間窗內(nèi)逐漸變成零值。

(3)為了消除儀器本身的影響，還原真實(shí)的地面運(yùn)動(dòng)，進(jìn)行了儀器響應(yīng)去除。

(4)本文研究的是淺層介質(zhì)波速的變化，關(guān)注相對(duì)高頻的信息，故選取濾波參數(shù)范圍為0.1～2Hz。

(5)時(shí)間域歸一化是數(shù)據(jù)預(yù)處理中最重要的步驟，為了減少地震事件和臺(tái)站附近非平穩(wěn)噪聲源相互關(guān)系的影響，采用一位歸一化方法(Derode et al，1999)進(jìn)行處理，采樣頻率為10Hz。

(6)譜白化即為頻率域歸一化，可以拓寬背景噪聲在互相關(guān)計(jì)算中的頻帶，也可以減小微震信號(hào)對(duì)計(jì)算結(jié)果的干擾。

2.2 背景噪聲數(shù)據(jù)的互相關(guān)及時(shí)間疊加

對(duì)2個(gè)臺(tái)站經(jīng)過(guò)預(yù)處理的數(shù)據(jù)進(jìn)行逐日互相關(guān)計(jì)算，獲得相應(yīng)的日互相關(guān)函數(shù)。選擇合適長(zhǎng)度天數(shù)的日互相關(guān)函數(shù)疊加作為當(dāng)前互相關(guān)函數(shù)(CCF)，其可表征一段時(shí)間的地下介質(zhì)狀態(tài)； 疊加整個(gè)研究時(shí)間范圍內(nèi)所有的日互相關(guān)函數(shù)，將其作為參考互相關(guān)函數(shù)(REF)，其可表征地下介質(zhì)的背景狀態(tài)。

本文背景噪聲互相關(guān)處理分為如下幾個(gè)步驟：

(1)將每日波形數(shù)據(jù)截取為48段，每段30min，每段有50%重疊；

(2)在頻率域?qū)?個(gè)臺(tái)站的每段數(shù)據(jù)進(jìn)行互相關(guān)運(yùn)算，疊加48段互相關(guān)函數(shù)作為日互相關(guān)函數(shù)；

(3)疊加8天日互相關(guān)函數(shù)作為CCF(圖2)；

(4)疊加2年所有日互相關(guān)函數(shù)作為REF(圖3)。

由不同疊加天數(shù)的互相關(guān)函數(shù)形態(tài)(圖2)可以看出，隨著疊加天數(shù)增加，CCF波形趨于穩(wěn)定； 而且對(duì)比圖3，CCF波形趨勢(shì)變化也更趨近于REF。本文選取8天的日互相關(guān)疊加作為CCF，用于后續(xù)解算。

圖2 JUN與JUX臺(tái)疊加不同天數(shù)的日互相關(guān)函數(shù)所得當(dāng)前互相關(guān)函數(shù)CCF

圖3 2009年1月1日—2011年1月1日J(rèn)UN與JUX臺(tái)參考互相關(guān)函數(shù)REF

圖4給出了2009年1月1—8日，共8天長(zhǎng)度的互相關(guān)疊加函數(shù)CCF與所有天數(shù)的REF對(duì)比，從中可以看出，兩者波形具有較好的相似度。2個(gè)互相關(guān)波形對(duì)稱性較弱，正時(shí)間段的振幅較高，負(fù)時(shí)間段的振幅較低，應(yīng)與噪聲源方位性分布有關(guān)(Chen et al，2010)； 短時(shí)間段CCF的噪聲源方位性效應(yīng)更明顯。

圖4 8天與全部698天的日互相關(guān)函數(shù)疊加所得CCF對(duì)比

2.3 到時(shí)差提取與相對(duì)速度變化計(jì)算

由圖4 可見(jiàn)，每對(duì)CCF與REF波形曲線在不同的時(shí)間t上有不同的走時(shí)延遲Δt。實(shí)際延時(shí)提取更多使用譜域方法進(jìn)行，可以明確相關(guān)函數(shù)中相關(guān)信號(hào)的帶寬(Poupinet et al，1984)，且能夠有效提升計(jì)算效率。該方法最早被用于研究地震對(duì)地下介質(zhì)影響所引起的相對(duì)速度變化。

本文相對(duì)速度變化計(jì)算使用的是移動(dòng)窗互譜法(MWCS，Moving Window Cross Spectral)。Clarke等(2011)詳細(xì)闡述了MWCS方法的原理及步驟，為利用背景噪聲研究地下介質(zhì)狀態(tài)變化打下堅(jiān)實(shí)基礎(chǔ)，本文不再贅述。

2.3.1 計(jì)算走時(shí)延遲變化

MWCS方法第一步為計(jì)算CCF與REF的走時(shí)延遲Δt。將每個(gè)互相關(guān)函數(shù)劃分為不同重疊的窗口N，進(jìn)行去平均及尖滅處理，并對(duì)CCF與REF分別進(jìn)行傅立葉變換，則互譜X(f)表示為

(1)

式中，F(xiàn)ref(f)與Fcur(f)分別表示REF與CCF的傅立葉變換，f為頻率，*表示復(fù)共軛。

走時(shí)延遲Δt在互譜中出現(xiàn)于相位譜上，因此將互譜X(f)進(jìn)一步寫為振幅譜|X(f)|與相位譜φ(f)相乘的形式，即

X(f)=|X(f)|eiφ(f)

(2)

式(2)中相位譜φ(f)在離散情況下可以寫為

φj=mfj

(3)

m=2πΔti

(4)

其中， Δti(i表示第i個(gè)窗口)即為2個(gè)波形信號(hào)特定時(shí)刻t的走時(shí)延遲；j=l，…，h?？苫谧钚《思訖?quán)線性回歸法獲得的斜率估計(jì)m(Clarke et al，2011)，有

(5)

其次，計(jì)算相關(guān)誤差em

(6)

(7)

其中，ωj為加權(quán)系數(shù)，有

(8)

其中，Cj為相關(guān)函數(shù)振幅譜的歸一化相干函數(shù)離散值，Xj為互譜離散值。

2.3.2 計(jì)算相對(duì)速度變化

在一階近似下，區(qū)域內(nèi)地震速度擾動(dòng)Δv/v是均勻的，并且表現(xiàn)為CCF相對(duì)于REF的拉伸Δt/t(Clarke et al，2011)。這種拉伸在數(shù)值上與速度擾動(dòng)相反(Poupinet et al，1984)，即

Δt/t=-Δv/v

(9)

綜合考慮臺(tái)站間距、波的速度、信號(hào)的能量強(qiáng)弱等因素，本文選擇加載移動(dòng)窗長(zhǎng)為25s，對(duì)應(yīng)于互相關(guān)10～35s部分，此部分包含面波及尾波。

3 結(jié)果與討論

3.1 2009—2010年地下介質(zhì)相對(duì)速度變化

圖5 相關(guān)系數(shù)與相對(duì)速度變化

3.2 臨近地下水位變化

地下水位數(shù)據(jù)來(lái)自于莒南魯14號(hào)井，按每分鐘水位測(cè)定記錄了2009年1月1日—2011年1月1日的水位值。分別選擇每天0點(diǎn)與12點(diǎn)的數(shù)據(jù)成圖，得到地下水水位變化，如圖6、圖7 所示。

圖6 2009年1月1日—2011年1月1日莒南魯14井水位變化(0點(diǎn)水位值)

圖7 2009年1月1日—2011年1月1日莒南魯14井水位變化(12點(diǎn)水位值)

由圖6(a)可見(jiàn)，水位呈整體上升趨勢(shì)，但有季節(jié)性變化。為了更清晰地看到變化規(guī)律，采用最小二乘法對(duì)數(shù)據(jù)進(jìn)行了去趨勢(shì)處理。由圖6(b)可見(jiàn)地下水水位呈明顯季節(jié)性變化，2009年4—10月、2010年6—10月地下水水位較低，隨時(shí)間呈上升趨勢(shì)； 2009年11月—2010年 5月水位值較高，呈下降趨勢(shì)。前人對(duì)不同地區(qū)水位研究也得到呈季節(jié)性變化的相似結(jié)論，如 Argus等(2014)利用GPS觀測(cè)方法推斷出加利福尼亞總蓄水量隨季節(jié)變化。

圖7可見(jiàn)整體水位仍呈上升趨勢(shì)，但去趨勢(shì)后水位季節(jié)性變化不如0點(diǎn)明顯，認(rèn)為白天水位變化受人為影響較大，夜間水位較為穩(wěn)定，誤差較小。

為了理解地下水位季節(jié)性的變化，統(tǒng)計(jì)了莒縣日降水量變化。由圖7(a)可見(jiàn)，2009年4—10月、2010年4—10月降水量大，其余時(shí)間降水量較小； 夏季降水量遠(yuǎn)高于冬季。與圖6(b)對(duì)比，認(rèn)為降水對(duì)地下水的影響有一定的滯后性。這里需要指出的是，由于地下水觀測(cè)點(diǎn)距離地震臺(tái)站的位置相對(duì)較遠(yuǎn)，降水對(duì)地下水位的影響僅具參考意義。

3.3 相對(duì)速度變化與地下水水位變化對(duì)比

將相對(duì)速度變化、地下水水位變化與日降水量變化進(jìn)行對(duì)比，如圖8 所示，從中可以發(fā)現(xiàn)，三者變化均呈明顯的季節(jié)性特征。相對(duì)速度變化與水位變化呈負(fù)相關(guān)的趨勢(shì)，但存在一定的時(shí)間上的延遲量，意味著地下水的影響不是即時(shí)的。例如，2009年10月—2010年4月水位呈下降趨勢(shì)，同時(shí)間段內(nèi)速度變化則先上升、后下降。

圖8 降水變化、地下水水位變化及相對(duì)速度變化對(duì)比

3.4 其他因素對(duì)地下介質(zhì)速度變化的影響

地下介質(zhì)速度變化不僅受到地下水位調(diào)整的影響，地震孕育過(guò)程、地球自轉(zhuǎn)和潮汐效應(yīng)等均會(huì)對(duì)其有一定的影響。

一般地震活動(dòng)相關(guān)的介質(zhì)速度變化更多地體現(xiàn)在地震發(fā)生前后(劉志坤等，2010)。臨沂地區(qū)地震活動(dòng)性不強(qiáng)，且本文獲得的地下介質(zhì)速度變化體現(xiàn)出季節(jié)性特征。

地下水潮汐效應(yīng)對(duì)地下水位變化會(huì)產(chǎn)生影響，但是目前的研究表明其周期為6天左右(王學(xué)靜等，2013)，這種短周期效應(yīng)需要更密集的臺(tái)網(wǎng)觀測(cè)才能夠給出。

另外，地球自轉(zhuǎn)可以引起地下介質(zhì)速度變化，但地球自轉(zhuǎn)帶來(lái)的季節(jié)性變化主要來(lái)源于太陽(yáng)輻射光在南、北半球表面上的分布不平衡(宋貫一，2011)，這種影響很弱，不足以解釋地下水位季節(jié)性變化的現(xiàn)象。

4 結(jié)論

本文基于2009—2010年的地震背景噪聲數(shù)據(jù)，利用移動(dòng)窗互譜法計(jì)算得到相對(duì)速度變化，結(jié)合降水和地下水水位監(jiān)測(cè)數(shù)據(jù)，研究了山東省臨沂地區(qū)地下介質(zhì)速度變化與地下水水位變化的相關(guān)性。研究結(jié)果表明，地下介質(zhì)速度變化呈季節(jié)性，夏季相對(duì)速度變化降低，冬季相對(duì)速度變化升高； 地下水水位變化亦呈季節(jié)性，夏季水位值較低，隨后逐漸增加，冬季水位值較高，隨后逐漸減少，受莒縣降水量季節(jié)性變化的影響； 介質(zhì)速度變化與地下水水位變化呈現(xiàn)此消彼長(zhǎng)的趨勢(shì)。介質(zhì)速度變化與地下水水位的明確相關(guān)性表明，可以基于地震背景噪聲對(duì)區(qū)域水位變化進(jìn)行監(jiān)控，從而服務(wù)于水資源利用及可能的地質(zhì)災(zāi)害預(yù)測(cè)。