近年来, 川西地区已成为青藏高原动力学和大陆内部块体边界动力学研究的热点, 与之相应的地球物理探测工作也十分活跃^[1~7].其中, 根据深地震测深剖面资料, 王椿镛等^{[4, 5]}研究了川西地区的地壳上地幔结构, 并用接收函数方法得到了川西高原壳内低速层的证据.Yao等^[8]利用中美合作的流动地震观测台网所记录的环境噪声及面波数据反演了川滇地区三维地壳上地幔速度结构.他们的结果显示龙门山断裂带是四川盆地和松潘-甘孜地块之间重要的地壳构造边界, 两侧的地壳结构有明显差异.但是, 限于台站覆盖的区域和空间分辨率, 已有的研究结果对于全面了解川西高原地壳上地幔速度结构的横向变化仍存在较大的局限性.鉴于川西地区处于研究青藏高原物质向东运移的关键部位, 其高分辨率的地壳上地幔三维速度结构探测研究对于青藏高原动力学和大陆内部块体边界动力学来说无疑具有至关重要的意义.

另一方面, 川西高原也是中国大陆地震活动最为频繁的一个地区.川西地区高分辨率的地壳上地幔结构对于研究这一地区地震成因背景和地震形成的机理具有不可替代的作用.特别是, 2008年5月12日在龙门山断裂带发生了汶川M_s8.0地震之后, 这进一步激发了人们对这一地区地震深部构造背景研究的兴趣.例如, 朱介寿^[9]依据多种地球物理资料分析了汶川地震及其邻域的岩石圈速度结构, 探讨了汶川地震的动力学背景.其他学者也陆续发表了一系列有关龙门山断裂带及其邻域的地壳上地幔三维速度结构的最新研究结果^[10~13].但是, 比较不同作者给出的结果后, 不难发现它们之间存在较为明显的差异和矛盾.

自2006年10月起, 在国家重大基础研究项目(973)的支持下, 中国地震局地质研究所地震动力学国家重点实验室在川西地区(26°N~32°N, 100°E~105°E)布设了由297台宽频带数字地震仪组成的流动观测台阵(以下简称川西台阵), 获得了汶川M_s8.0地震前后完整的地震活动的观测记录.这为利用天然地震数据研究川西地区三维高分辨率地壳上地幔速度结构图像积累了新的宝贵资料^[14].利用川西台阵观测的远震P波波形数据和非线性接收函数方法, 刘启元等^[15]已给出了穿越汶川地震震中区(沿31°N)的二维地壳上地幔S波速度结构及地壳的平均泊松比.利用远震P波走时层析成像方法, 郭飚等^[16]给出了川西地区(29°N~32°N, 100°E~105°E) 400km深度范围内地壳上地幔的三维P波速度结构.

近年来, 环境噪声的层析成像技术获得了快速发展, 并成为广泛关注的热点^{[17, 18]}.这为利用密集台阵进行地下速度结构的高分辨率地震成像提供了新的技术途径^{[19, 20]}.本文的目的在于, 利用川西台阵记录的环境噪声数据和相应的相速度层析成像技术, 研究川西台阵下方的三维地壳结构.由于各种地球物理反演方法往往存在程度不同的非唯一性, 对同一区域采用不同数据源和多种地球物理反演技术研究其地下速度结构无疑成为解决这一问题最有效的途径.本文得到的结果为川西高原和四川盆地地壳结构的研究提供了新的独立观测证据, 并为进一步的深入研究提供了重要约束.本文所用的北部156个台站主要覆盖龙门山断裂带及其邻区, 这与南部的安宁河断裂带地区的构造有较大差别.本文的重点仅限于讨论北部的成像结果.

图 1给出了在29°N~32°N和100°E~105°E区域范围内川西台阵的台站分布情况.在此研究区内, 共计有156个台站.其台站间距约为5~40km.刘启元等^[14]已详细介绍了川西台阵的台站观测系统, 不再赘述.如图 1所示, 这些台站覆盖了川滇菱形块体、松潘-甘孜地块及属于华南地块的四川盆地, 而鲜水河和龙门山断裂带则构成了它们之间的地质构造边界.图 1同时给出了陈九辉等^[21]利用双差定位技术给出的汶川地震主震及部分余震的重新定位结果.本文的工作主要利用2007年1~12月川西台阵记录的环境噪声数据.

图 1

Fig. 1

图 1 川西台阵的台站分布 红色三角表示地震台站，红色圆表示主震，空心圆表示余震，黑色实线表示断层；F1:龙门山断裂带，F2:鲜水河断裂带，F3:安宁河断裂带，F4:理塘断裂带；1:汶川一茂县断裂，f2:映秀一北川断裂，f3:灌县一江油断裂 Fig. 1 Station map of the western Sichuan seismic array Red triangles represent the stations, red circle represents the main shock, white circles represent aftershocks, black solid line represents the fault; F1:Longmen Shan fault, F2 : Xianshuihe fault, F3 : Anninghe fault, F4 : Litang fault; f1 : Wenchuan-Maoxian fault, f2 : Yingxiu-Beichuan fault, f3 : Guanxian-Jiangyou fault

若干文献已就环境噪声层析成像的理论和方法做过详细讨论^[22~28].这些研究结果是本文工作的理论基础.本文的数据处理主要依据Yao等^[19]发展的由环境噪声数据测量瑞利波相速度频散曲线及相速度层析成像的方法.它主要包括以下步骤: (1)数据预处理; (2)环境噪声垂直分量波形数据的互相关和叠加; (3)瑞利波频散曲线的测量和误差估计; (4)二维相速度层析成像.其原理和参数选择可简述如下.

首先, 以十天作为一个数据段, 对每个数据段分别进行重采样(5 Hz)和去倾处理.经过带宽0.5~60s的带通滤波后, 数据需要进一步通过滑动绝对平均进行时间域归一化处理, 其中时窗长度可以取带通滤波最大周期的一半.需要强调的是, 上述预处理过程对于去除地震信号及台站附近非平稳噪声的影响是十分必要的.

理论上^[25~28], 对于给定的观测台站A和B, 环境噪声的互相关函数与格林函数之间的关系可由

(1)

给出.这里, C_AB(t)为台站A和台站B环境噪声的互相关函数; G_AB(t)为台站B的精确格林函数(以A为源); G_BA(-t)则是台站A的逆向格林函数(以B为源); Ĝ_AB(t)和Ĝ_BA(-t)表示相应的估计值(经验格林函数), t表示时间.

由于格林函数的空间互易性, 利用

(2)

可由环境噪声的互相关函数测量台站B的经验格林函数(以台站A为源).

研究表明^[25~28]:如果噪声源在空间上分布是均匀的, 那么噪声数据正向和逆向互相关的结果应该是对称的, 这样得到的经验格林函数和精确格林函数的相位信息是相同的.但实际上, 噪声源的空间分布通常并不均匀, 以至于不能保证由实际观测数据得到的正向和逆向的互相关序列的对称性.为此, Lin等^{[29, 30]}建议用噪声数据正向及逆向互相关函数的反序平均结果计算相应的经验格林函数.他们证明这不但能够有效地降低噪声源分布不均匀造成的影响, 而且可同时提高结果的信噪比.本文工作中先将互相关函数求时间导数, 再将得到的正向及逆向求导结果做反序平均, 最终可获得等价的经验格林函数.图 2给出了利用上述方法由川西台阵环境噪声数据得到的部分台站为期一年的经过带通滤波的反序平均前的瑞利波经验格林函数.图 2表明, 它们不但显示了很好的对称性, 而且具有较高的信噪比.

图 2

Fig. 2

图 2 川西台阵部分台站的瑞利波经验格林函数纵坐标为台站对的间距 Fig. 2 Rayleigh wave EGFs from ambient noise recorded at some station pairs in the western Sichuan seismic array The vertical axis shows the inter-station distance.

需要说明的是: (1)噪声源方位分布的不均会使经验格林函数产生相移, 但一般较少, 这样所引起的相速度测量误差一般在1%以内^[31]; (2)式(2)中的导数计算会使经验格林函数与相应的互相关函数之间产生π/2的相位差, 虽然这对群速度的测量不会有影响, 但对于相速度必须给予具体考虑^[19]; (3)由环境噪声提取的经验格林函数存在振幅谱的畸变, 这对于依赖振幅信息的群速度频散曲线测量是不利的^{[19, 23]}, 但经验格林函数的相位信息与真实格林函数的相位信息是非常相似的, 并能够通过图像转换技术获得准确可靠的相速度频散曲线^[19].

姚华建等^{[32, 33]}曾提出过一种基于图像分析的相速度频散曲线快速提取方法.与传统的等值线描点方法相比, 该方法不但可以快速追踪整条频散曲线, 而且能够提高相速度频散曲线的测量精度.本文的相速度频散曲线测量采用了同样方法.

基于时间域经验格林函数的相速度频散曲线测量^[19]需要满足相应的远场近似条件

(3)

这里, c_AB为平均相速度, T为测量周期, λ为波长, Δ为台站对的间距.实际上, 平均相速度是T的函数, 它可以由

(4)

给出, 其中t为经验格林函数在周期T附近窄带滤波后波峰的传播时间.

另外, 根据(3)式, 我们可以确定不同台站对所能得到的最大周期范围, 并确定提取相速度频散曲线的时间窗.当周期大于T时, 由于噪声源的季节性变化将导致频散曲线测量结果稳定性变差.由于川西台阵台站间的最大间距约为550km, 我们提取频散曲线的最大周期限于40s以内.

由于川西台阵的台站对间距变化较大(从几公里到550多公里不等), 为了有效地提取短周期相速度频散曲线, 依据台站对的间距, 我们将其经验格林函数大致分为两类: (1)对于Δ < 100km的情况, 面波群速度窗设定为1.0~3.75km/s, 并取相应的频散分析的最大周期为20s; (2)对于Δ≥100km的情况, 面波群速度窗设定为1.75~4.5km/s, 并取相应的频散分析的最大周期为50s.这种做法不仅可以提高效率, 而且有助于准确提取有效信号.

如图 1所示, 川西台阵覆盖区域的地形复杂, 且地形高差变化很大.尤其是在龙门山断裂带的两侧, 在水平距离仅为50~60km的范围内, 海拔高程由约400m陡增到约4000m.对于这样巨大的地形变化, 在计算台站间距时必须考虑相应的高程修正.本文的台站间距校正由

(5)

给出.其中, Δ₀为由台站经纬度计算得到的距离, δh为台站对的高程差, Δ_c为校正后的台站间距.

经过筛选, 在2~40s范围内, 对156个台站, 我们共得到11358条瑞利波相速度频散曲线.这里所用的频散曲线是通过人工挑选得到的, 没有定量的标准.图 3给出了各个周期的射线路径数目.图 3表明, 我们在各个周期内都获得了足够多的频散曲线.图 4则给出了不同周期(2s、8s、12s、18s、25s和35s)瑞利波相速度测量射线路径分布.图 4表明, 不同周期瑞利波射线路径对本文研究的区域均有密集采样, 而这对进行地壳不同深度高分辨率相速度层析成像非常重要.

图 3

Fig. 3

图 3 不同周期瑞利波射线路径统计 Fig. 3 Number of the ray paths of Rayleigh waves at different periods

图 4

Fig. 4

图 4 不同周期瑞利波相速度测量的射线路径PathN为射线路径数. Fig. 4 Ray-path coverage for the Rayleigh phase velocity measurements at different periods PathN represents the number of the ray paths.

作为一个例子, 图 5给出了由台站对KCD09-KDF07的环境噪声数据得到的瑞利波相速度频散曲线.其校正后的台站间距为298km, 瑞利波相速度频散曲线相应的最大周期为28s.由图 5b可见, 经过噪声数据正向及逆向波形的反序平均明显地改善了相应的经验格林函数的对称性.

图 5

Fig. 5

图 5 瑞利波相速度频散曲线测量示意图 (a)台站分布示意图，红色三角（KCD09)和蓝色三角(KDF07)表示本示意图所用的台站对；（b)红色波形和蓝色波形分别表示经验格林函数的正逆部分，其中虚线波形和实线波形分别表示加时间窗前后的波形，黑色波形表示正逆波形反序叠加后得到的经验格林函数的对称形式；(c)相速度频散曲线测量（黑色实线）图，红色星号标示出所提取的频散点. Fig. 5 Example of inter-station Rayleigh wave phase velocity dispersion curve (a) Station map with the station pair of KCD09 (red triangle) and KDF07 (blue triangle) showed in this example; (b) Red and blue waveforms are positive and negative components of EGF, respectively; Meanwhile the dashed and solid lines are EGFs before and after windowing, respectively; The black waveform is the summation of positive and negative components of EGF; (c) Phase velocity dispersion measurements (black solid line) from the EGF with red asterisks as the extracted dispersion points.

利用由观测台阵噪声数据得到的瑞利波经验格林函数及相速度频散曲线, 我们可以进行台阵覆盖区域下方地壳结构的相速度层析成像.本文的瑞利波相速度层析成像采用Tarantola等^{[34, 35]}发展的反演方法.仿照Yao等^[19], 相应的目标函数由

(6)

给出.其中, t和s分别为预测的相走时和相慢度矢量, s_p为先验相慢度矢量, t_obs为观测相位走时向量, C_D和C_M分别为数据和模型协方差矩阵.模型空间的协方差可由

(7)

给出.其中, σ_s表示先验慢度的不确定性, σ_c表示先验相速度的不确定性, L为模型的相关长度.

令σ_c等于各个周期相速度值的两倍标准偏差, c₀取不同周期相速度的平均结果, 模型的相关长度l=c₀T/2, T为相速度的周期.对于短周期(2~10s)数据, 如果L过小, 有可能导致虚假异常.为此, 我们设定L=min (30km, c₀T/2), 以避免出现奇异点.

在本文的相速度层析成像中, 我们采用了Montagner^[36]发展的区域化方法, 并按照0.25°× 0.25°对研究区域(29°N~32°N, 100°E~105°E)进行网格划分, 得到周期2~35s的相速度扰动图像.这个方法不但在射线路径覆盖密集的区域可以给出高精度的相速度图像, 而且可用于研究介质的方位各向异性.但是, 本文将仅限于介质各向同性的情况.

图 6给出了利用前文描述的相速度层析成像技术得到的相速度分布图像.研究表明^{[37, 38]}, 基阶瑞利波相速度对大约1/3波长深度附近介质的S波速度结构最为敏感, 且周期越长, 其垂向分辨率越差.对于图 6, 我们可有如下的观察:

图 6

Fig. 6

图 6 不同周期瑞利波相速度分布 红色三角表示地震台站，空心圆表示主震，白色实线表示断层.图的顶部给出了相应的周期和平均相速度. Fig. 6 Rayleigh-wave phase velocity maps at different periods Red triangles represent the stations, white circle represents the main shock, white solid lines represent the faults. The period and average phase velocity are given at the top of each figure, respectively.

(1) 图 6a和图 6b分别给出了T=2s和8s的瑞利波相速度分布图像.大体上, 它们分别反映了地壳2km和8km深度范围内的平均速度结构的横向变化.对比图 1可知, 图 6a的瑞利波相速度横向变化与地表构造特征是吻合的.不难看出, 研究区近地表及地壳上部速度结构呈现出明显受断层控制的区块特征.特别是, 以龙门山断裂带为界, 松潘-甘孜地块与四川盆地的相速度分别具有明显的高速和低速异常.这种差异应主要起因于四川盆地前陆存在较厚的(约10km左右)沉积层^[39].在上述深度范围内, 各区块内的浅表及地壳上部速度结构均呈现不同程度的横向不均匀性.随着周期的增大, 这种横向不均匀性总体上明显变弱.与此同时, 在鲜水河断裂带和理塘断裂带, 相速度表现出不均匀的低速异常, 而龙门山断裂带则表现出不均匀的高速异常, 并在龙门山断裂带、鲜水河断裂带和安宁河断裂带交界处呈现出较复杂的速度结构.

(2) 图 6c和图 6d分别给出了T=12s和18s的瑞利波相速度分布图像.大体上, 它们分别反映了上地壳12km和19km深度范围内速度结构的横向分块不均匀的特征, 在龙门山断裂带, 则显示较为一致的高速异常.在T=12s的图像上, 四川盆地内的低速异常主要分布于盆地前陆靠近龙门山断裂带的区域, 并在31°N以北的区域显示了较大的异常.在T=18s的图像上, 松潘-甘孜地块, 特别是川滇地块呈现出明显的低速异常, 表明这两个地块内壳内低速体的顶部边界应该位于大约19km的深度.这与刘启元等^[15]接收函数反演的结果基本一致.

(3) 图 6e和图 6f分别给出了T=25s和35s的瑞利波相速度分布图像.大体上, 它们分别反映中下地壳(或上地幔顶部, 如在四川盆地)28km和42km深度范围内速度结构的横向变化.其最为突出的特征是, 四川盆地的中下地壳呈现明显的高速异常, 而松潘-甘孜地块, 特别是川滇地块, 呈现明显的低速异常.沿龙门山断裂带, 以汶川地震的震中为界, 相速度分布呈现出明显的南北分段现象.上述特征与郭飚等^[16]给出的P波走时层析成像结果基本吻合.

综合不同周期相速度分布, 不难看出, 当周期大于18s时, 四川盆地表现为高速异常, 而松潘-甘孜地块和川滇地块表现为低速异常, 特别是川滇地块具有比松潘-甘孜地块更低的相速度, 表明四川盆地的中地壳和下地壳强度应明显高于川滇地块和松潘-甘孜地块相应深度的地壳强度.这与刘启元等^[15]根据接收函数反演所做的推断是一致的.

需要指出是, 周期35s的相速度图像显示, 在松潘-甘孜地块和川滇地块内, 存在跨鲜水河断裂的大面积低速异常.这也可能是, 随着周期的增大和路径数量的减少, 相速度层析成像的分辨率明显变差造成的一种假象.

为了检验本文相速度层析成像的空间分辨能力, 利用文献[17]的检测板测试技术, 图 7给出了检测板测试结果.考虑到研究区复杂的地形地貌及地壳结构的横向变化, 对于不同周期, 我们给出不同的初始速度和速度扰动.具体取值如表 1所示.

图 7

Fig. 7

图 7 检测板测试 (a)，(b)初始模型；(c)~(h)不同周期的检测板测试结果. Fig. 7 Checkerboard tests (a), (b) Inputmodels; (c)~(h) Inversion results at different periods.

表 1(Table 1)

表 1 检测板测试的初始模型设定 Table 1 Input models for checkerboard tests

表 1 检测板测试的初始模型设定 Table 1 Input models for checkerboard tests

由图 7可见, 射线路径的分布直接影响反演恢复的结果.对比图 6可以看出, 只要保证射线路径分布合理, 我们目前所用的反演方法能够得到比较可靠而且稳定的结果.短周期2~8s的良好恢复情况说明, 即使遇到浅层速度变化较大的地区(如本文所研究的川西地区), 只要保证射线路径分布合理, 我们目前所用的反演方法能够得到比较可靠而且稳定的结果.随着周期的增大, 射线路径分布呈现出中心密集、边缘稀疏的形态.与此同时, 只有检测板中间区域可以较好地恢复.图 7h表明, 对于周期35s的情况, 除了射线相对密集的有限区域尚有一定的恢复能力外, 其他部分的恢复能力均比较差.尽管如此, 一般来说, 利用本文采用的相速度成像方法在30km深度范围内仍可以给出较为可靠的地壳速度成像的结果.我们认为, 原则上, 密集地震观测台阵记录的环境噪声数据对于研究上地壳及中上地壳的速度结构具有重要的应用价值.

综上所述, 利用近年来新发展起来的环境噪声面波层析成像方法和2007年1~12月川西台阵(29°N~32°N, 100°E~105°E)区域记录环境噪声数据, 我们给出了川西地区周期2~35s瑞利波相速度分布.我们的结果表明, 观测台阵覆盖的川滇地块、松潘-甘孜地块和四川盆地的地壳速度结构具有显著差异和各自不同的特点.对于不同周期的相速度分布, 它们具体表现为:

(1) 短周期(2~8s)相速度分布与地表构造特征相吻合, 作为川滇地块、松潘-甘孜地块和四川盆地之间的边界断裂, 龙门山断裂带和鲜水河断裂带对上述三个地块上地壳的速度结构具有明显的控制作用, 四川盆地前陆低速特征表明相应区域存在较厚的(约10km)沉积盖层;

(2) 中周期(12~18s)相速度分布表明, 川滇地块和松潘-甘孜地块中上地壳速度结构存在明显的不均匀横向变化, 并形成了尺度不同、高低速相间的分块结构, 而四川盆地中地壳整体上表现出相对高速;

(3) 长周期(25~35s)相速度分布表明, 松潘-甘孜地块, 特别是川滇地块中下地壳表现为广泛的明显低速异常, 意味着它们的中下地壳相对软弱, 而四川盆地的中下地壳呈现整体性的相对高速, 意味着四川盆地具有相对坚硬的中下地壳, 并且以汶川地震的震中为界, 龙门山断裂带的地壳结构显示了北段为高速异常、南段为低速异常的分段特征;

(4) 周期大于12s相速度分布表明, 随着周期的加大, 龙门山断裂带与鲜水河断裂带及安宁河断裂带交汇处东侧高速异常区逐渐增大, 表明川滇地块运动方向的改变应与四川盆地坚硬地壳的阻挡作用密切有关.

本文的结果表明, 密集台阵观测的环境噪声数据对研究区高分辨率地壳速度结构的层析成像具有重要的应用价值和前景.与传统面波方法相比, 由于能够从环境噪声数据中提取短周期面波频散数据, 这种方法的优势在于可以显著提高上地壳速度结构成像的空间分辨率.但是, 对于密集台阵观测系统来说, 由于台站间距的限制, 从环境噪声数据中提取可靠的周期较长的面波频散数据目前尚有一定的困难.除非与其他地震观测数据结合, 这种方法目前尚难以独立用于超过一般地壳深度范围的地震层析成像.

本文是中国地震局地质研究所与美国麻省理工学院合作研究项目的部分成果.本文第1和第3作者感谢访美工作期间所有美方人员给予的帮助.