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

1)中國科學院大學，地球與行星科學學院，北京 100049

2)中國地震臺網(wǎng)中心，北京 100045

0 引言

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

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

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

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

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

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

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

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

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

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

(3)日降水數(shù)據(jù)由莒縣氣象觀測站(編號54936)所記錄，數(shù)據(jù)來自國家氣象科學數(shù)據(jù)中心。

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

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

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

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

(1)為了使數(shù)據(jù)標準化，進行除均值來去除波形數(shù)據(jù)本身具有的非零均值。通常波形數(shù)據(jù)會存在一個長周期的線性趨勢，從而影響數(shù)據(jù)的分析，故需進行去趨勢處理。

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

(3)為了消除儀器本身的影響，還原真實的地面運動，進行了儀器響應去除。

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

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

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

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

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

本文背景噪聲互相關處理分為如下幾個步驟：

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

(2)在頻率域?qū)?個臺站的每段數(shù)據(jù)進行互相關運算，疊加48段互相關函數(shù)作為日互相關函數(shù)；

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

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

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

圖2 JUN與JUX臺疊加不同天數(shù)的日互相關函數(shù)所得當前互相關函數(shù)CCF

圖3 2009年1月1日—2011年1月1日JUN與JUX臺參考互相關函數(shù)REF

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

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

2.3 到時差提取與相對速度變化計算

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

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

2.3.1 計算走時延遲變化

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

(1)

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

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

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

(2)

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

φj=mfj

(3)

m=2πΔti

(4)

其中， Δti(i表示第i個窗口)即為2個波形信號特定時刻t的走時延遲；j=l，…，h。可基于最小二乘加權線性回歸法獲得的斜率估計m(Clarke et al，2011)，有

(5)

其次，計算相關誤差em

(6)

(7)

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

(8)

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

2.3.2 計算相對速度變化

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

Δt/t=-Δv/v

(9)

綜合考慮臺站間距、波的速度、信號的能量強弱等因素，本文選擇加載移動窗長為25s，對應于互相關10～35s部分，此部分包含面波及尾波。

3 結果與討論

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

圖5 相關系數(shù)與相對速度變化

3.2 臨近地下水位變化

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

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

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

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

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

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

3.3 相對速度變化與地下水水位變化對比

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

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

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

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

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

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

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

4 結論

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