利用背景噪声数据提取地震台站间面波的可靠性分析——以中国大陆中东部地区的宽频带台站为例＊

郑 现 赵翠萍 周连庆 郑斯华

（中国北京100036中国地震局地震预测研究所）

引言

对随机噪声场中两接收点记录的噪声进行互相关计算，可以提取两个接收点间结构体的脉冲响应（格林函数）.这种思想最早在超声波的实验室测试中得到验证（Weaver，Lobkis，2001；Lobkis，Weaver，2001）.之后很快被应用于地震学领域，通过台站对记录到的远震尾波（Campillo，Paul，2003）和背景噪声（Shapiro，Campillo，2004）的互相关计算，证明了散射波场中可以提取地震面波.其中，利用背景噪声提取面波由于不依赖于地震，资料获取不受时间和空间的限制，在层析成像应用方面较有优势.

Shapiro等（2005）成功地将利用背景噪声提取的格林函数展开层析成像的研究，并得到了美国加州地区与已知地质构造吻合较好的7.5s和15s的瑞雷波群速度分布图像.随之在世界各地区开展了一系列利用背景噪声进行层析成像的研究.研究尺度从小区域到整个大陆（Kang，Shin，2006；Lin et al，2007；Yang et al，2007；Bensen et al，2008；Lin et al，2008；Zheng et al，2008；Sun et al，2010），研究内容包括瑞雷波和勒夫波的群速度和相速度.背景噪声层析成像的周期范围通常在8—40s范围内，但Kang和Shin（2006）得到了2—4s的瑞雷波群速度图像；Nishida等（2009）利用台站记录到的“地球哼鸣”（earth’s hum）信号进行互相关分析，也得到全球范围周期在300s以上的层析成像结果.在数据处理方面，Bensen等（2007）提出了一套比较完整灵活、适应性好的处理方法，此后的噪声层析成像研究多参照了这一套方法.另外，根据背景噪声不依赖震源的优势也可以对地下介质进行连续监测.其方法的主要原理是通过背景噪声提取经验格林函数，对波形的尾波进行干涉测量得到地下介质的相对速度变化.目前，该方法已经应用于地震（Xu，Song，2009；Chen et al，2010；刘志坤，黄金莉，2010；赵盼盼，2011）、火山（Brenguier et al，2008a）和断层（Wegler，Sens-Schonfelder，2007；Brenguier et al，2008b）等方面的监测.

背景噪声层析成像法克服了地震面波层析成像法的一些局限性，首先它不受地震数量不足、分布不均的制约，只要有数量足够、分布均匀的宽频带地震台，在一些地震活动性较弱的地区也可以展开高精度的成像研究.此外，由于不依赖震源信息，避免了定位不精确带来的误差.该方法在短周期成像结果的分辨率比地震面波层析成像法有着绝对的优势，能够更好地反映地壳浅层的结构特征.

由于面波层析成像方法的优越性，近些年在中国大陆也开展了一系列区域性（Yao et al，2006；Li et al，2009；Fang et al，2010；Zheng et al，2010；李昱等，2010）和大陆范围（Zheng et al，2008；Sun et al，2010）的层析成像工作，取得了一定的成果.本研究利用中国中东部地区100个分布均匀的宽频带地震台记录到的21个月的连续波形数据，进行互相关和叠加计算，提取瑞雷波经验格林函数（empirical Green function，简写为EGF），利用时频分析法（Levshin，Ritzwoller，2001）提取台站对间的频散曲线，并采用2种检测方法对经验格林函数和频散曲线的可靠性进行评估.此外也对叠加时长和时间稳定性进行了分析，对使用不同长度时间序列互相关进行叠加后EGF大于10的路径数目进行了统计和讨论，以期为今后的研究提供借鉴.

1 数据与处理

本研究从中国大陆中东部地区国家台网和区域台网（图1）中，挑选出100个分布基本均匀的宽频带台站（平均台间距200km），利用这些台站记录的背景噪声数据，进行面波可靠性研究.所使用的台站包括国家地震台网的61个台站和区域地震台网的39个台站（郑秀芬等，2009）.地震计的主体类型是CTS-1、KS2000、BBVS-120和CMG-3ESPB，绝大部分仪器的频带宽度的低频端为120s，只有少数仪器的可达360s或仅为60s.本文使用这些台站21个月（2009年1月—2010年9月）的垂直分量连续波形数据，经过1Hz重采样后进行研究（郑现等，2012）.

图1 中国大陆中东部地区主要构造单元及台站分布图灰色粗线勾勒出主要块体，蓝色细线表示主要断层，白色三角为本研究所使用的国家台网台站和区域台网台站.OB：鄂尔多斯盆地；SCB：四川盆地；JB：江汉盆地；NJB：苏北盆地Fig.1 Major geological structures and seismic stations in central and eastern China Gray thick lines mark major blocks，blue thin lines represent major faults and white triangles represent stations of China national and regional networks.OB：Ordos Basin；SCB：Sichuan Basin；JB：Jianghan Basin；NJB：north Jiangsu Basin

数据处理的方法主要分为3个步骤（Bensen et al，2007）：① 单台数据处理；② 互相关和叠加计算；③ 频散曲线提取.

1.1 单台数据处理

单台数据处理的目的主要是消除地震信号、仪器故障等因素对背景噪声的影响.首先把数据分成长度为一天的数据单元，对数据进行去均值、去零漂，并截取成同样长度的记录.然后进行5—150s带通滤波和时域归一化.本研究根据Bensen等（2007）对归一化方法所做的讨论，使用滑动绝对值平均法对所使用的数据进行处理.该方法在固定长度的时间窗中，滑动计算波形的绝对值的平均值，权重的计算公式为

并用原始时间序列乘以权重的倒数，得到新的时间序列d′n＝dn／wn.本研究的归一化时间窗（2L＋1）设定为40s.最后，为了抑制某些信号的干扰，拓宽背景噪声的频带宽度，对数据进行频谱白化处理.

1.2 互相关和叠加计算

对所有台站对的数据进行互相关和叠加计算可得到台站对间的经验格林函数（EGF）.经验格林函数是沿时间轴双向传播的，由于台站两侧噪声源方位角分布不均匀，波形正负时间坐标轴的分量振幅是不对称的.为了避免这种不均匀性的影响并增加EGF的信噪比（SNR）（Bensen et al，2007；Lin et al，2008），我们对沿正负时间轴的波形取均值，从而得到经验格林函数的“对称”分量，在接下来的计算中只使用该对称分量.图2给出了以内蒙台网NMJIN台站为中心台，与其它台站组成的台站对间的经验格林函数沿时间轴双向和取平均后的结果.在整个研究区，路径达2 000km以上的范围内，瑞雷波的信号都十分清晰.

图2 以NMJIN台为中心的互相关波形（经过10—40s滤波）（a）双向传播的互相关波形；（b）互相关波形的对称分量Fig.2 Wave cross-correlations between station NMJIN and other stations（band-pass filtering 10—40s）（a）Positive and negative lags of the cross-correlations；（b）Symmetric component of the cross-correlations

1.3 群速度频散曲线提取

本研究使用时频分析法（Levshin，Ritzwoller，2001），采用反频散滤波或相位匹配滤波去除潜在干扰，获得了100个台站两两台站之间的瑞雷波群速度频散曲线.图3a中我们给出了3组穿过不同构造块体的大弧路径，即穿过四川盆地和华南造山系的SCSPA－－GDSHD，穿过秦岭造山带和华南造山系的GSTSS－－FZCM2，以及穿过华北－渤海湾盆地和华南造山系的GDHUJ－－HEZJK.图3b，d表明，这3条路径由经验格林函数提取的频散曲线与由参考的3D剪切波速度模型（Shapiro，Ritzwoller，2002）提取的频散曲线在形态上基本一致.其中，沿GDSHD－－SCSPA路径和GDHUJ－－HEZJK的频散曲线与理论曲线不完全吻合，恰恰反映了这条具体路径介质速度与区域平均值的差异.

图3 瑞雷波经验格林函数及频散曲线示意图（图3d引自郑现等，2012）（a）与（b），（c），（d）对应的射线路径示意图；（b），（c），（d）为沿（a）图射线路径由时频分析法提取的瑞雷波群速度频散结果（蓝色曲线）和由全球剪切波速度模型（Shapiro，Ritzwoller，2002）计算的理论值（黑色曲线）Fig.3 Example of Rayleigh wave EGFs and dispersion measurements（Fig.3bis quoted from Zheng et al，2012）（a）Location map for（b），（c），（d）showing the 3selected ray paths；（b），（c），（d）Rayleigh wave dispersion measurement（blue curve）retrieved by frequency-time analysis and the prediction caculate（black curve）from the 3-D global shear velocity model（Shapiro，Ritzwoller，2002）along ray paths in（a）

从理论上讲，100个台站会构成4 950［n（n－1）／2］个台站对（Bensen et al，2007）.但是，为了确保结果的可靠性，我们基于信噪比测试和台间距对所有频散曲线进行了筛选.对各周期只保留SNR大于10和台间距大于3倍波长的记录（Bensen et al，2007）.

2 质量评估

背景噪声频散曲线的质量是利用其进一步开展速度结构反演或速度变化监测的关键.因此，在进行层析成像前，需要做一些检测来评估满足上述截断准则的频散曲线的质量.本文采用的检测方法主要有：沿相同路径的格林函数与地震信号对比、空间一致性检测.其中，经验格林函数与地震记录对比检测最为有效，也是对利用噪声提取面波理论的一个验证.

2.1 经验格林函数与地震信号对比

对于经验格林函数与地震记录对比检测，最理想的情况是地震的位置正好与台站位置重合（Bensen et al，2007）.依据这一准则，从中东部地区2009年1月—2010年9月发生的4.5级以上地震中选取了重庆荣昌台附近的荣昌地震（发震时刻2010年9月10日，105.43°E、29.36°N，震源深度6km，MS＝4.5），该地震距荣昌台约2km.我们挑选出相对于这次地震与荣昌台不同方位角和震中距的其它5个台站（图4a），首先分别对这5个台站记录到的荣昌MS4.5

图4 经验格林函数与地震信号对比图（a）荣昌地震（红色圆）和台站CQROC与其它5个台站间的射线路径；（b）沿图（a）射线路径的经验格林函数（绿色线）与5个台站记录到的地震波形（红色线）的对比结果（以20s为中心频率滤波），其中左侧数字为震中距，右侧字母为台站代号；（c）与（b）图类似，经过了10—25s的带通滤波Fig.4 Comparison between EGFs and earthquake seismic records（a）Location map showing the Rongchang earthquake（red circle）and the paths between station CQROC and five other stations；（b）Comparison between EGFs（green lines）and earthquake seismic records（red lines）along ray paths in（a）（after narrowband filtering centered at 20speriod）；（c）Similar to（b），after 10—25sband-pass filtering

地震波形记录进行去仪器响应校正；然后分别对提取的荣昌台与这5个台站之间的EGF和地震记录进行以20s为中心频率滤波（图4b）和10—25s的带通滤波（图4c），并进行归一化处理；最后由于不同路径长度差别很大，为了方便观察，将信号由时间域转换至群速度域.从图4b，c可以看出，各台站记录的地震信号与EGF基本一致，在峰值部分吻合得比较好.这证明利用背景噪声提取的经验格林函数是台站对之间路径上的面波，是真实可靠的.

2.2 空间一致性检测

空间一致性测试是通过比较空间相似路径的形态，来检查频散曲线的稳定性和可靠性.本文挑选出了两个位置邻近的台站SXXAX和SXLIF（间隔108km），它们和一较远的台站JXSHR构成两组相似路径.作为对比，这里另给出了台站JXSHR－－SCSPA路径（图5a），该路径与上述两组路径差别较大.图5b表明，JXSHR－－SXXAX、JXSHR－－SXLIF两组路径间的频散曲线在整个观测周期（8—50s）形态基本相似，而JXSHR-SCSPA的频散曲线则与其差别较大且形态则更为复杂，反映出此路径上介质的强烈非均匀性.这从一个侧面证明了本研究提取的频散曲线的可靠性.

图5 频散曲线的空间连续性实验（a）实验选取的3组射线路径：JXSHR－－SCSPA（长度约1400km）、JXSHR－－SXLIF、JXSHR－－SXXAX（长度约1000km），其中SXLIF和SXXAX的间距为108km；（b）沿（a）图射线路径的频散曲线对比结果Fig.5 Example of spatial consistency of dispersion curves（a）The 3selected paths：JXSHR－－SCSPA （about 1400km long），JXSHR－－SXLIF，JXSHR－－SXXAX（about 1000km long），where the distance between SXLIF and SXXAX is about 108km；（b）Comparison of the dispersion curves along the paths in（a）

2.3 叠加时长及时间稳定性分析

下面通过统计互相关叠加的数据长度与满足一定信噪比的频散曲线数目的关系，来分析叠加时长的选取.

图6 各周期所使用的频散曲线数目随时间序列长度的分布Fig.6 Dispersion measurements’number versus time-series length at different periods

我们对研究区内的所有台站对分别进行了时间长度为3，4，5，6，9，12，15，18和21个月的叠加计算，并基于信噪比和台间距对这些频散曲线进行筛选.各周期只保留SNR大于10且台间距大于3倍波长的记录，统计结果如图6所示.对于本研究所用的中国大陆中东部地区的100个台站分布和仪器频带宽度，相同叠加时长，周期为15，20，10，25，8，30，35，40和45s范围内的频散曲线数目依次大幅减小，如叠加12个月时，15s周期的频散曲线数目可达4 000条，而40s周期的频散曲线数目仅有不到1 500条.随着叠加时间长度的增加，沿各路径经验格林函数的信噪比增高，各周期（8—45s）筛选后用于反演的频散曲线数目增加.叠加长度从3个月增加到12个月时，频散曲线数目上升比较快；叠加长度从12个月增加到21个月时，频散曲线数目上升趋势相对较缓.根据这一现象，我们认为，当利用噪声对中国大陆中东部地区进行面波成像时，叠加时间最好取12个月或以上.

在此基础上，我们选择了SCCD2－－ANKG3和SXTIY－－JXJIJ两组路径（图7），来分析叠加时长为12个月的数据的稳定性.以1个月为滑动窗，对21个月的数据共提取到10条频散曲线.此外，我们也给出了21个月时长数据提取的频散曲线.由图7可见，SXTIY－－JXJIJ路径在周期8—50s范围内整体一致性较好；SCCD2－－ANKG3路径的频散曲线在30s以下周期一致性均较好，提取的EGF和频散曲线基本稳定，受噪声源的影响较小；在30s以上周期，不同时间段的叠加结果比较分散，与21个月的频散曲线也存在不同程度的偏差，说明在该周期范围内频散曲线一定程度上受到了噪声源分布不均的影响.因此，在利用噪声数据滑动提取EGF，分析其时间变化特征以达到监视介质性质变化目的的研究中，由于计算EGF所使用的数据通常比较短，所以应使用30s以下周期的EGF波形进行相关分析.30s以下较短周期提取的EGF对应的大致是地壳深度30km以内的面波，基本上是大部分大陆地区孕震层的下边界.

3 讨论与结论

本文以中国大陆中东部100个分布均匀的宽频带地震台为例，由台站记录到的21个月的垂直分量的连续波形数据，进行互相关和叠加计算，并由时频分析法提取了各台站对间的经验格林函数和频散曲线.通过质量评估，我们检测了利用背景噪声在中东部地区提取经验格林函数和瑞雷面波频散曲线的可靠性和稳定性，对叠加时长也进行了讨论，主要结论为：

图7 频散曲线的时间稳定性检测（a）实验所选取的两组路径；（b）对每条路径分别使用10组长度为12个月的数据（红色线）、21个月的数据（黑色线）互相关叠加计算，得到频散曲线Fig.7 Example of temporal consistency of dispersion curves（a）The 2selected paths；（b）Dispersion measurements calculated by using 12months’data（10lines，all in red）and 21months data（black）separately along the paths in（a）

1）在质量评估测试中，沿不同长度、不同方位角路径的经验格林函数和地震记录波形在峰值附近基本一致，表明由背景噪声记录进行互相关叠加后可以得到地下介质结构的格林函数.而沿空间上近似路径提取的频散曲线、不同时间段数据提取的频散曲线形态大体一致，也说明本研究得到的格林函数和频散曲线是可靠的.

2）当叠加的时间长度由3个月增加到21个月时，经验格林函数的信噪比增大，各周期可以使用的频散曲线数目也随之增多.根据频散曲线数目随叠加时长的变化趋势，我们认为，叠加时长至少为12个月才能得到数目足够和可靠的频散曲线，来进行进一步的面波速度结构成像研究.

3）本文结果表明，叠加时间长度为12个月时，可以确保在30s以下周期沿各路径的频散曲线时间上稳定.对于30s以上周期，提取频散曲线的叠加时间要足够长，应在12个月以上.作者认为，联合天然地震的长周期面波进行成像，弥补噪声在长周期的不可靠性，以取得各个周期可靠的面波和剪切波结构，是目前的研究方向.

中国地震局地球物理研究所“国家数字测震台网数据备份中心”为本文提供了地震波形数据.文中所有图件的绘制均使用了Wessel和Smith 1998年的GMT软件.在此表示诚挚的谢意.

李昱，姚华健，刘启元，陈九辉，van der Hilst R D，李顺成，黄慧，郭飚，王峻，齐少华.2010.川西地区台阵环境噪声瑞利波相速度层析成像［J］.地球物理学报，53（4）：842－－852.

刘志坤，黄金莉.2010.利用背景噪声互相关研究汶川地震震源区地震波速度变化［J］.地球物理学报，53（4）：853－－863.

赵盼盼.2011.汶川地震区地壳速度相对变化的地震环境噪声自相关研究［D］.北京：中国地震局地质研究所：1－－50.

郑秀芬，欧阳飚，张东宁，姚志祥，梁建宏，郑洁.2009.“国家数字测震台网数据备份中心”技术系统建设及其对汶川大地震研究的数据支撑［J］.地球物理学报，52（5）：1412－－1417.

郑现，赵翠萍，周连庆，郑斯华.2012.中国大陆中东部地区基于背景噪声的瑞利波层析成像［J］.地球物理学报，55（6）：1919－－1928.

Bensen G B，Ritzwoller M H，Barmin M P，Levshin A L，Lin F，Moschetti M P，Shapiro N M，Yang Y.2007.Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements［J］.Geophys J Int，169（3）：1239－－1260.

Bensen G B，Ritzwoller M H，Shapiro N M.2008.Broadband ambient noise surface wave tomography across the United States［J］.J Geophys Res，113（B5306）：1－－21.

Brenguier F N，Shapiro N M，Campillo M，Ferrazzini V，Duputel Z，Coutant O，Nercessian A.2008a.Towards forescasting volcanic eruptions using seismic noise［J］.Nat Geosci，1（2）：126－－130.

Brenguier F N，Campillo M，Hadziioannou C，Shapiro N M，Nadeau R M，Larose E.2008b.Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations［J］.Science，321（5895）：1478.Campillo M，Paul A.2003.Long-range correlations in the diffuse seismic coda［J］.Science，299（5606）：547－－549.

Chen J H，Froment B，Liu Q Y，Campillo M.2010.Distribution of seismic wave speed changes associated with the 12 May 2008 MW7.9Wenchuan earthquake［J］.Geophys Res Lett，37（18）：L18302，doi：10.1029／2010GL044582.

Fang L H，Wu J P，Ding Z F，Panza G F.2010.High resolution Rayleigh wave group velocity tomography in North China from ambient seismic noise［J］.Geophys J Int，181（2）：1171－－1182.

Kang T S，Shin J S.2006.Surface-wave tomography from ambient seismic noise of accelerograph networks in southern Korea［J］.Geophys Res Lett，33（L17303）：1－－5.

Levshin A L，Ritzwoller M H.2001.Automated detection，extraction，and measurement of regional surface waves［J］.Pure Appl Geophys，158（8）：1531－－1545.

Li H，Su W，Wang C，Huang Z.2009.Ambient noise Rayleigh wave tomography in western Sichuan and eastern Tibet［J］.Earth Planet Sci Lett，282（1－－4）：201－－211.

Lin F，Ritzwoller M H，Townend J，Bannister S，Savage M K.2007.Ambient noise Rayleigh wave tomography of New Zealand［J］.Geophys J Int，170（2）：649－－666.

Lin F，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［J］.Geophys J Int，173（1）：281－－298.

Lobkis O I，Weave R L.2001.On the emergence of the Green’s function in the correlations of a diffuse field［J］.J Acoust Soc Amer，110（6）：3011－－3017.

Nishida K，Montagner J P，Hitoshi K.2009.Global surface wave tomography using seismic hum［J］.Science，326（5949）：112.

Shapiro N M，Ritzwoller M H.2002.Monte-Carlo inversion for a global shear-velocity model of the crust and upper mantle［J］.Geophys J Int，151（1）：88－－105.

Shapiro N M，Campillo M.2004.Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise［J］.Geophys Res Lett，31（7）：1－－5.

Shapiro N M，Campillo M，Stehly L，Ritzwoller M H.2005.High-resolution surface-wave tomography from ambient seismic noise［J］.Science，305（5681）：233－－236.

Sun X L，Song X D，Zheng S H，Yang Y J，Ritzwoller M H.2010.Three dimensional shear wave velocity structure of the crust and upper mantle beneath China from ambient noise surface wave tomography［J］.Earthq Sci，23（5）：449－－463.

Weaver R L，Lobkis O I.2001.Ultrasonics without a source Thermal fluctuation correlations at MHz frequencies［J］.Phys Rev Lett，87（13）：134301.

Wegler U，Sens-Schonfelder C.2007.Fault zone monitoring with passive image interferometry［J］.Geophys J Int，168（3）：1029－－1033.

Xu Z J，Song X D.2009.Temporal changes of surface wave velocity associated with major Sumatra earthquakes from ambient noise correlation［J］.Proc Natl Acad Sci USA，106（34）：14207－－14212，doi：10.1073／pnas.0901164106.

Yang Y，Ritzwoller M H，Levshin A L，Shapiro N M.2007.Ambient noise Rayleigh wave tomography across Europe［J］.Geophys J Int，168（1）：259－－274.

Yao H J，van der Hilst R D，de Hoop M V.2006.Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-I.Phase velocity maps［J］.Geophys J Int，166（2）：732－－744.

Zheng S，Sun X，Song X，Yang Y，Ritzwoller M H.2008.Surface wave tomography of China from ambient seismic［J］.Geochem Geophys Geosyst，9（5）：1525－－2027.

Zheng Y，Yang Y，Ritzwoller M H，Zheng X，Xiong X，Li Z.2010.Crustal structure of the northeastern Tibetan plateau，the Ordos block and the Sichuan basin from ambient noise tomography［J］.Earthq Sci，23（5）：465－－476.