2. 中国地震台网中心, 北京 100045
2. China Earthquake Networks Center, Beijing 100045, China
大型含水层系统中的地下水常被作为水源,是水循环中的重要组成部分(Kim et al,2019)。地下水质量和水位变化监测主要基于钻孔等定点直接测量(张磊等,2021;门北方,2018),其成本相对较高,因而可靠的低成本间接手段越来越受到研究人员和生产部门的关注。目前,基于合成孔径雷达(InSAR)、全球定位系统(GPS)或大地测量等方法观测到的变形模式或重力异常调整,可以在区域尺度监测地下水的赋存和水位变化等(Amelung et al,1999;Argus et al,2014;Kim et al,2019;李嘉等,2020),但是对于科学研究和生产实践,区域性地下水水位变化监测在空间尺度上仍然偏大。
地震波对于地下介质结构的变化比较敏感(Huang et al,2010;Mordret et al,2010;Yeh et al,2013)。降水、河流补给及地下水抽采等对区域地下水调整的影响可以在地震波速结构变化中有较好的反映(Mainsant et al,2012)。固定地震台站连续观测记录可以持续记录非地震的地震动,即背景噪声。连续观测的地震背景噪声互相关函数可用于持续监测地下速度结构的变化(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)。
地震背景噪声测量已经被引入到监测地下水位调整引起的地下速度结构变化中(Kim et al,2019)。目前,国内应用地震背景噪声监测水库水位相关的波速变化情况取得了一定的成果。例如,陈蒙等(2013)运用背景噪声技术对云南省大银甸水库周围波速变化进行研究,发现水位变化与波速变化存在明显的相关性,并认为这种相关性是由水库水体的卸载作用造成的;安艳茹等(2015)研究了紫坪铺水库区域在蓄水与泄水过程中库区介质的变化特征,表明地下介质的相对波速变化与水位变化存在较为明显的相关性,且在时间上有一定延迟,认为其可能与水的渗透作用有关。
本文基于相对圈闭的临沂地区国家测震台网JUX与JUN台站记录的连续波形数据,利用背景噪声互相关技术研究地下介质波速变化,并结合临沂地区地下水水位和降水等观测数据,分析降水和地下水抽取等诸多因素影响下的地下水位变化与地震波速变化的关系。
1 研究区域及资料选取山东省临沂地区位于鲁中南低山丘陵区东南部和鲁东丘陵南部,地势西北高、东南低,自北向南有沂山、蒙山、尼山等3条主要山脉,呈NE-SW向延伸,发育沂河与沭河(图 1)。区域水系相对圈闭(李致家等,2005),地下水位变化主要由于天然降水和工、农业使用。
![]() |
图 1 研究区域概况 蓝色三角为地震台站JUN与JUX;黑色圆圈为台站所在县城;No.14为莒南鲁14号井;黑色方块为莒县气象观测站(编号54936) |
本地区有国家数字测震台网的多个宽频带地震台、中国地震局地下水位动态监测台莒南鲁14号井以及中国气象局的莒县气象观测站等。本文收集并使用了如下3种类型的数据资料:
(1) 地震波形资料为2009年1月1日—2011年1月1日山东省临沂市莒南县的JUN台与日照市莒县的JUX台的连续记录,2个地震台站间距为34.5km。数据资料来自于国家数字测震台网数据备份中心(郑秀芬等,2009),2个台站数据在2009、2010年分别有19天和13天出现缺失或损坏,可用数据的天数总计698天。
(2) 地下水水位连续数据来自莒南鲁14号井,该井为研究区内距离最近的与地震监测相关的水井。
(3) 日降水数据由莒县气象观测站(编号54936)所记录,数据来自国家气象科学数据中心。
2 地震数据处理本文使用NoisePy软件(Jiang et al,2020)处理数据,并计算相对速度变化。
2.1 连续记录数据预处理数据预处理(Bensen et al,2007)主要包括均值和趋势变化去除、波形尖灭、仪器响应去除、带通滤波、时间域归一化及谱白化处理等。具体如下:
(1) 为了使数据标准化,进行除均值来去除波形数据本身具有的非零均值。通常波形数据会存在一个长周期的线性趋势,从而影响数据的分析,故需进行去趋势处理。
(2) 在对数据进行谱域操作(如FFT、滤波等)时,若数据的两端不为零,则会出现谱域假象,因而实际数据经常需要做尖灭处理,使得数据两端在短时间窗内逐渐变成零值。
(3) 为了消除仪器本身的影响,还原真实的地面运动,进行了仪器响应去除。
(4) 本文研究的是浅层介质波速的变化,关注相对高频的信息,故选取滤波参数范围为0.1~2Hz。
(5) 时间域归一化是数据预处理中最重要的步骤,为了减少地震事件和台站附近非平稳噪声源相互关系的影响,采用一位归一化方法(Derode et al,1999)进行处理,采样频率为10Hz。
(6) 谱白化即为频率域归一化,可以拓宽背景噪声在互相关计算中的频带,也可以减小微震信号对计算结果的干扰。
2.2 背景噪声数据的互相关及时间叠加对2个台站经过预处理的数据进行逐日互相关计算,获得相应的日互相关函数。选择合适长度天数的日互相关函数叠加作为当前互相关函数(CCF),其可表征一段时间的地下介质状态;叠加整个研究时间范围内所有的日互相关函数,将其作为参考互相关函数(REF),其可表征地下介质的背景状态。
本文背景噪声互相关处理分为如下几个步骤:
(1) 将每日波形数据截取为48段,每段30min,每段有50%重叠;
(2) 在频率域对2个台站的每段数据进行互相关运算,叠加48段互相关函数作为日互相关函数;
(3) 叠加8天日互相关函数作为CCF(图 2);
![]() |
图 2 JUN与JUX台叠加不同天数的日互相关函数所得当前互相关函数CCF (a)4天叠加;(b)8天叠加;(c)32天叠加;所有互相关波形均进行滤波处理 |
(4) 叠加2年所有日互相关函数作为REF(图 3)。
![]() |
图 3 2009年1月1日—2011年1月1日JUN与JUX台参考互相关函数REF |
由不同叠加天数的互相关函数形态(图 2) 可以看出,随着叠加天数增加,CCF波形趋于稳定;而且对比图 3,CCF波形趋势变化也更趋近于REF。本文选取8天的日互相关叠加作为CCF,用于后续解算。
图 4给出了2009年1月1—8日,共8天长度的互相关叠加函数CCF与所有天数的REF对比,从中可以看出,两者波形具有较好的相似度。2个互相关波形对称性较弱,正时间段的振幅较高,负时间段的振幅较低,应与噪声源方位性分布有关(Chen et al,2010);短时间段CCF的噪声源方位性效应更明显。
![]() |
图 4 8天与全部698天的日互相关函数叠加所得CCF对比 红色虚线表示8天日互相关函数叠加所得CCF;黑色实线表示全部698天日互相关函数叠加所得CCF;虚框内为40~55s放大图;波峰时差为Δt |
由图 4可见,每对CCF与REF波形曲线在不同的时间t上有不同的走时延迟Δt。实际延时提取更多使用谱域方法进行,可以明确相关函数中相关信号的带宽(Poupinet et al,1984),且能够有效提升计算效率。该方法最早被用于研究地震对地下介质影响所引起的相对速度变化。
本文相对速度变化计算使用的是移动窗互谱法(MWCS,Moving Window Cross Spectral)。Clarke等(2011)详细阐述了MWCS方法的原理及步骤,为利用背景噪声研究地下介质状态变化打下坚实基础,本文不再赘述。
2.3.1 计算走时延迟变化MWCS方法第一步为计算CCF与REF的走时延迟Δt。将每个互相关函数划分为不同重叠的窗口N,进行去平均及尖灭处理,并对CCF与REF分别进行傅立叶变换,则互谱X(f)表示为
$ X(f)=F_{\text {ref }}(f) \times F_{\text {cur }}^{*}(f) $ | (1) |
式中,Fref(f)与Fcur(f)分别表示REF与CCF的傅立叶变换,f为频率,*表示复共轭。
走时延迟Δt在互谱中出现于相位谱上,因此将互谱X(f)进一步写为振幅谱|X(f)|与相位谱φ(f)相乘的形式,即
$ X(f)=|X(f)| \mathrm{e}^{i \varphi(f)} $ | (2) |
式(2)中相位谱φ(f)在离散情况下可以写为
$ \begin{gathered} \varphi_{j}=m f_{j} \end{gathered} $ | (3) |
$ m=2 \pi \Delta t_{i} $ | (4) |
其中,Δti(i表示第i个窗口)即为2个波形信号特定时刻t的走时延迟;j=l,…,h。可基于最小二乘加权线性回归法获得的斜率估计m(Clarke et al,2011),有
$m=\frac{\sum\limits_{j=l}^{h} \omega_{j} f_{j} \varphi_{j}}{\sum\limits_{j=l}^{h} \omega_{j} f_{j}^{2}} $ | (5) |
其次,计算相关误差em
$ e_{m}=\sqrt{\sum_{j}\left(\frac{\omega_{j} f_{j}}{\sum_{i} \omega_{i} f_{i}^{2}}\right)^{2} \sigma_{\varphi}^{2}} $ | (6) |
式(5)、(6)中
$\sigma_{\varphi}^{2}=\frac{\sum_{j}\left(\varphi_{j}-m f_{j}\right)^{2}}{N-1} $ | (7) |
其中,ωj为加权系数,有
$ \omega_{j}=\sqrt{\frac{C_{j}^{2}}{1-C_{j}^{2}} \sqrt{\left|X_{j}\right|}} $ | (8) |
其中,Cj为相关函数振幅谱的归一化相干函数离散值,Xj为互谱离散值。
2.3.2 计算相对速度变化在一阶近似下,区域内地震速度扰动Δv/v是均匀的,并且表现为CCF相对于REF的拉伸Δt/t(Clarke et al,2011)。这种拉伸在数值上与速度扰动相反(Poupinet et al,1984),即
$ \Delta t / t=-\Delta v / v $ | (9) |
综合考虑台站间距、波的速度、信号的能量强弱等因素,本文选择加载移动窗长为25s,对应于互相关10~35s部分,此部分包含面波及尾波。
3 结果与讨论 3.1 2009—2010年地下介质相对速度变化运用MWCS方法计算得到2009年1月1日—2011年1月1日JUN与JUX台相对速度变化结果,如图 5所示。经过多次测试,选取8天日互相关叠加作为CCF最为合适,避免了速度变化点过于密集而影响分析。由图 5 (a)可以看出,CCF与REF相关系数较高,约在0.85以上,结果较为可靠,误差置信水平为0.96。由图 5 (b)可看出,2009年4—9月、2010年5—10月相对速度变化降低,2009年1—3月、2009年10月—2010年4月相对速度变化升高;速度变化幅度约为±0.2%;变化整体呈波浪形特征,夏季相对速度变化降低,冬季相对速度变化升高,说明地下介质波速变化与季节变化相关(Meier et al,2010)。值得一提的是,2010年11—12月正值冬季,而相对速度变化趋势与2009年11—12月相比恰恰相反,可能与地震等自然灾害发生或噪声源受到较大影响有关(刘志坤等,2010;Hobiger et al,2014;Taira et al,2015;刘国明,2015)。
![]() |
图 5 相关系数与相对速度变化 (a)CCF与REF相关系数;(b)运用MWCS方法计算的相对速度变化结果;圆点表示相对速度变化值,短棒表示误差 |
地下水位数据来自于莒南鲁14号井,按每分钟水位测定记录了2009年1月1日—2011年1月1日的水位值。分别选择每天0点与12点的数据成图,得到地下水水位变化,如图 6、图 7所示。
![]() |
图 6 2009年1月1日—2011年1月1日莒南鲁14井水位变化(0点水位值) (a)含趋势;(b)去趋势 |
![]() |
图 7 2009年1月1日—2011年1月1日莒南鲁14井水位变化(12点水位值) (a)含趋势;(b)去趋势 |
由图 6(a)可见,水位呈整体上升趋势,但有季节性变化。为了更清晰地看到变化规律,采用最小二乘法对数据进行了去趋势处理。由图 6(b)可见地下水水位呈明显季节性变化,2009年4—10月、2010年6—10月地下水水位较低,随时间呈上升趋势;2009年11月—2010年5月水位值较高,呈下降趋势。前人对不同地区水位研究也得到呈季节性变化的相似结论,如Argus等(2014)利用GPS观测方法推断出加利福尼亚总蓄水量随季节变化。
图 7可见整体水位仍呈上升趋势,但去趋势后水位季节性变化不如0点明显,认为白天水位变化受人为影响较大,夜间水位较为稳定,误差较小。
为了理解地下水位季节性的变化,统计了莒县日降水量变化。由图 7(a)可见,2009年4—10月、2010年4—10月降水量大,其余时间降水量较小;夏季降水量远高于冬季。与图 6(b)对比,认为降水对地下水的影响有一定的滞后性。这里需要指出的是,由于地下水观测点距离地震台站的位置相对较远,降水对地下水位的影响仅具参考意义。
3.3 相对速度变化与地下水水位变化对比将相对速度变化、地下水水位变化与日降水量变化进行对比,如图 8所示,从中可以发现,三者变化均呈明显的季节性特征。相对速度变化与水位变化呈负相关的趋势,但存在一定的时间上的延迟量,意味着地下水的影响不是即时的。例如,2009年10月—2010年4月水位呈下降趋势,同时间段内速度变化则先上升、后下降。
![]() |
图 8 降水变化、地下水水位变化及相对速度变化对比 (a)莒县区域日降水量变化;(b)2009—2010年区域0点水位变化;(c)相对速度变化;短棒表示误差棒,红色实线表示趋势 |
地下介质速度变化不仅受到地下水位调整的影响,地震孕育过程、地球自转和潮汐效应等均会对其有一定的影响。
一般地震活动相关的介质速度变化更多地体现在地震发生前后(刘志坤等,2010)。临沂地区地震活动性不强,且本文获得的地下介质速度变化体现出季节性特征。
地下水潮汐效应对地下水位变化会产生影响,但是目前的研究表明其周期为6天左右(王学静等,2013),这种短周期效应需要更密集的台网观测才能够给出。
另外,地球自转可以引起地下介质速度变化,但地球自转带来的季节性变化主要来源于太阳辐射光在南、北半球表面上的分布不平衡(宋贯一,2011),这种影响很弱,不足以解释地下水位季节性变化的现象。
4 结论本文基于2009—2010年的地震背景噪声数据,利用移动窗互谱法计算得到相对速度变化,结合降水和地下水水位监测数据,研究了山东省临沂地区地下介质速度变化与地下水水位变化的相关性。研究结果表明,地下介质速度变化呈季节性,夏季相对速度变化降低,冬季相对速度变化升高;地下水水位变化亦呈季节性,夏季水位值较低,随后逐渐增加,冬季水位值较高,随后逐渐减少,受莒县降水量季节性变化的影响;介质速度变化与地下水水位变化呈现此消彼长的趋势。介质速度变化与地下水水位的明确相关性表明,可以基于地震背景噪声对区域水位变化进行监控,从而服务于水资源利用及可能的地质灾害预测。
安艳茹、张晓东, 2015, 利用背景噪声研究紫坪铺水库水位加卸载对库区地下介质波速的影响, 中国地震, 31(4): 616-628. DOI:10.3969/j.issn.1001-4683.2015.04.002 |
陈蒙, Hillers G, 王宝善, 等, 2013. 基于地震背景噪声的与水库水位变化相关的地下波速变化的研究. 见: 中国地球物理2013——第十二分会场论文集. 昆明: 中国地球物理学会, 2.
|
李嘉、唐河、饶维龙等, 2020, 南水北调工程对华北平原水储量变化的影响, 中国科学院大学学报, 37(6): 775-783. |
李致家、李志龙、孔祥光等, 2005, 沂沭河水系水文模型与洪水预报研究, 水力发电, 31(7): 25-27. DOI:10.3969/j.issn.0559-9342.2005.07.008 |
刘国明, 2015, 利用背景噪声研究地震波速度变化及其在长白山火山监测中的应用, 地震地质, 37(2): 576-587. |
刘志坤、黄金莉, 2010, 利用背景噪声互相关研究汶川地震震源区地震波速度变化, 地球物理学报, 53(4): 853-863. DOI:10.3969/j.issn.0001-5733.2010.04.010 |
门北方, 2018, 大庆油田钻井对地下水环境的影响, 石油石化节能, 8(5): 45-46. |
宋贯一, 2011, 地球自转速度季节性变化的主要原因解析, 地球物理学进展, 26(2): 450-455. |
王学静、李海龙、柳富田等, 2013, 利用地下水潮汐效应估算含水层参数及预测水位变化, 工程勘察, 41(9): 28-31, 89. |
温扬茂、高松、许才军, 2019, 利用双台站背景噪声分析2017年墨西哥MW7.1地震震源区的地震波速变化, 地球物理学报, 62(8): 3024-3033. |
张磊、任妹娟、何计彬等, 2021, 一孔多层监测井中地下水分层监测系统开发, 桂林理工大学学报, 41(2): 291-295. |
郑秀芬、欧阳飚、张东宁等, 2009, "国家数字测震台网数据备份中心"技术系统建设及其对汶川大地震研究的数据支撑, 地球物理学报, 52(5): 1412-1417. |
Amelung F, Galloway D L, Bell J W, et al, 1999, Sensing the ups and downs of Las Vegas: InSAR reveals structural control of land subsidence and aquifer-system deformation, Geology, 27(6): 483-486. DOI:10.1130/0091-7613(1999)027<0483:STUADO>2.3.CO;2 |
Argus D F, Fu Y N, Landerer F W, 2014, Seasonal variation in total water storage in California inferred from GPS observations of vertical land motion, Geophys Res Lett, 41(6): 1971-1980. DOI:10.1002/2014GL059570 |
Bensen G D, Ritzwoller M H, Barmin M P, et al, 2007, Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys J Int, 169(3): 1239-1260. DOI:10.1111/j.1365-246X.2007.03374.x |
Brenguier F, Campillo M, Takeda T, et al, 2014, Mapping pressurized volcanic fluids from induced crustal seismic velocity drops, Science, 345(6192): 80-82. DOI:10.1126/science.1254073 |
Chen J H, Froment B, Liu Q Y, et al, 2010, Distribution of seismic wave speed changes associated with the 12 May 2008 MW7.9 Wenchuan earthquake, Geophys Res Lett, 37(18): L18302. |
Clarke D, Zaccarelli L, Shapiro N M, et al, 2011, Assessment of resolution and accuracy of the moving window cross spectral technique for monitoring crustal temporal variations using ambient seismic noise, Geophys J Int, 186(2): 867-882. DOI:10.1111/j.1365-246X.2011.05074.x |
Derode A, Tourin A, Fink M, 1999, Ultrasonic pulse compression with one-bit time reversal through multiple scattering, J Appl Phys, 85(9): 6343-6352. DOI:10.1063/1.370136 |
Gerstoft P, Shearer P M, Harmon N, et al, 2008, Global P, PP, and PKP wave microseisms observed from distant storms, Geophys Res Lett, 35(23): L23306. DOI:10.1029/2008GL036111 |
Hadziioannou C, Larose E, Baig A, et al, 2011, Improving temporal resolution in ambient noise monitoring of seismic wave speed, J Geophys Res, 116(B7): B07304. |
Hobiger M, Wegler U, Shiomi K, et al, 2014, Single-station cross-correlation analysis of ambient seismic noise: application to stations in the surroundings of the 2008 Iwate-Miyagi Nairiku earthquake, Geophys J Int, 198(1): 90-109. DOI:10.1093/gji/ggu115 |
Huang H, Yao H J, Van Der Hilst R D, 2010, Radial anisotropy in the crust of SE Tibet and SW China from ambient noise interferometry, Geophys Res Lett, 37(21): L21310. |
Jiang C X, Denolle M A, 2020, NoisePy: A new high-performance Python tool for ambient-noise seismology, Seismol Res Lett, 91(3): 1853-1866. DOI:10.1785/0220190364 |
Kim D, Lekic V, 2019, Groundwater variations from autocorrelation and receiver functions, Geophys Res Lett, 46(23): 13722-13729. DOI:10.1029/2019GL084719 |
Lin F C, Moschetti M P, Ritzwoller M H, 2008, Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps, Geophys J Int, 173(1): 281-298. DOI:10.1111/j.1365-246X.2008.03720.x |
Mainsant G, Larose E, Brönnimann C, et al, 2012, Ambient seismic noise monitoring of a clay landslide: toward failure prediction, J Geophys Res Earth Surf, 117(F1): F01030. |
Meier U, Shapiro N M, Brenguier F, 2010, Detecting seasonal variations in seismic velocities within Los Angeles Basin from correlations of ambient seismic noise, Geophys J Int, 181(2): 985-996. |
Mordret A, Jolly A D, Duputel Z, et al, 2010, Monitoring of phreatic eruptions using interferometry on retrieved cross-correlation function from ambient seismic noise: results from Mt. Ruapehu, New Zealand, J Volcanol Geoth Res, 191(1-2): 46-59. DOI:10.1016/j.jvolgeores.2010.01.010 |
Poupinet G, Ellsworth W L, Frechet J, 1984, Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California, J Geophys Res, 89(B7): 5719-5731. DOI:10.1029/JB089iB07p05719 |
Shapiro N M, Campillo M, 2004, Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, Geophys Res Lett, 31(7): L07614. |
Taira T, Brenguier F, Kong Q K, 2015, Ambient noise-based monitoring of seismic velocity changes associated with the 2014 MW6.0 South Napa earthquake, Geophys Res Lett, 42(17): 6997-7004. DOI:10.1002/2015GL065308 |
Yeh Y L, Wen S, Lee K J, et al, 2013, Shear-wave velocity model of the Chukuo fault zone, Southwest Taiwan, from cross correlation of seismic ambient noise, J Asian Earth Sci, 75: 174-182. DOI:10.1016/j.jseaes.2013.07.023 |
Zheng S H, Sun X L, Song X D, et al, 2008, Surface wave tomography of China from ambient seismic noise correlation, Geochem Geophys Geosyst, 9(5): Q05020. |