图1 SsPmp的射线路径示意图
Fig. 1 Raypath diagram of SsPmp wave
摘要 提出虚拟地震台阵双差测深法。通过多道互相关, 测量密集观测台阵中相邻台站间SsPmp波与Ss波到时差的差值, 然后求解最小二乘意义下每个台站SsPmp波与Ss波的到时差。该方法利用台阵优势, 可以有效地降低到时差测量误差, 有利于准确地估计地壳厚度。将该方法应用在模拟数据和“中国地震科学台阵探测计划”台阵数据上, 测量精度显著提高。
关键词 虚拟地震测深法(VDSS); 台阵; 多道互相关; 双差; Ss到时异常
莫霍面是地球的一级界面, 其深度和起伏对推测区域地质结构和演化过程有重要意义。莫霍面深度的探测方法主要有人工震源反射法[1‒4]和远震接收函数法[5‒7]。其中, 接收函数法因成本低、环境污染小而被广泛采用。H-κ叠加[8]是接收函数的常用方法, 给定平均P波速度后, 将Ps波与后面的多次震相PpPs, PpSs和PsPs等的理论到时对应的振幅叠加, 据此估计地壳厚度和P/S平均波速比。然而, 在地壳结构变化剧烈的区域以及有巨厚松散沉积物的区域, 多次转换波的震相不清晰, 造成P波接收函数类方法估计地壳参数的不确定度较大[9‒10]。
虚拟地震测深方法(virtual deep seismic soun-ding, VDSS)是近年发展起来的测量地壳厚度的有效方法。该方法利用同一台站远震SsPmp波与Ss波的到时差(图1)估计地壳厚度。当震中距在30°～55°之间时, SsPmp波是全反射震相, 同时没有Ss波三叉震相的干扰, 振幅大, 信噪比高, 受壳内界面和局部异常体的影响较小, 使用少量地震的资料即可有效地约束地壳结构。VDSS方法已经成功地应用于青藏高原及其周缘地区[11‒16]、华北克拉通[17‒19]以及美国西部[20]、美国东南部[21]和澳大利亚中 部[22‒23]等地区的地壳结构探测中。
图1 SsPmp的射线路径示意图
Fig. 1 Raypath diagram of SsPmp wave
在单层均匀介质的假设情况下, 单台记录中SsPmp波与Ss波的到时差与地壳厚度的转换关系如下式所示:
(1)
其中, Tvdss为同一台站SsPmp波与Ss波的到时差, H为地壳厚度, vp为地壳平均P波速度, p为射线 参数。
VDSS方法通常采用波形拟合法[11,17,22]测量Tvdss。首先根据先验知识, 给定地壳和上地幔的平均密度、P波速度和S波速度, 然后改变地壳厚度, 模拟得到一系列介质响应, 再分别与震相分离估计的Ss子波[24]进行卷积, 得到模拟波形。选取与实际记录互相关系数最大的模拟波形, 将该模型的地壳厚度代入式(1), 计算得到Tvdss的理论值, 作为该台站的测量结果。该方法具有一定的局限性, 当给定不同的先验速度模型时, 波形拟合法测量的最佳拟合厚度存在较大的差异。史克旭等[15]的测试结果显示, 当地壳速度从6.0km/s上升至6.6km/s时, 最佳拟合厚度会增加11.2km。同时, 在地质结构复杂的区域, 波形拟合法的最佳拟合波形不能很好地拟合实际观测值, 得到的Tvdss误差较大。
密集台阵的布设, 不仅增加了数据, 也促进了新方法的产生和应用。Yu等[20]利用多道互相关方法计算每个台站的Ss波到时异常(即实际到时与理论到时之差), 并对到时差与地壳厚度转换加以校正, 消除Ss波实际到时与理论到时差异导致的地壳厚度估计偏差。他们的校正公式如下:
(2)
其中, ΔtSs为Ss波到时异常校正项, 为图1中台站接收点与地表转换点之间Ss波到时异常的差值。魏晓拙等[14]推导出测线垂直于地质块体边界走向及震源位于测线延长线的特殊情况下, Ss波到时异常的表达形式。Liu等[19]按照实际到时, 将全部台站的Ss波对齐后取均值作为子波, 降低了子波估计误差导致的到时误差。
本文提出虚拟地震台阵双差测深法(array-based double difference virtual deep seismic sounding,ABDDVDSS), 基于台阵数据, 首先利用多道互相关计算相邻台站间的Tvdss之差ΔTvdss,ij。ΔTvdss,ij是到时差之差, 是一个双差; 然后求解双差方程组, 计算最小二乘意义下每个台站的Tvdss; 最后计算地壳厚度。本文将此方法应用在模拟数据和“中国地震科学台阵探测计划”台阵数据上, 以期显著地降低Tvdss的误差。
SsPmp震相在莫霍界面发生全反射时, 相对于入射Ss波存在相移, SsPmp的波形与Ss不同。在单层均匀介质的假设情况下, 这种相移随着地壳P波速度、上地幔P波速度和射线参数的增加而降低, 与地壳厚度无关[18,25]。
当两个台站距离较近时, 地壳P波速度、上地幔P波速度和射线参数的变化较小, 两个台站的SsPmp波相移之差小, 波形相似度较高, 适合采用互相关方法计算两个台站间SsPmp波的到时差。互相关方法计算相邻台站间同一震相的到时差, 不需要先验速度模型和Ss子波波形等信息。
同一台站对之间的Ss到时差与SsPmp到时差不同, 因此需要选取合适的时窗, 分别截取只包含Ss震相和SsPmp震相的数据进行互相关, 测量同一震相在台站对之间的到时差, 降低测量误差。将全部台站对间Ss震相的到时差方程联立, 可以求解每个台站在最小二乘意义下的Ss到时异常, 并可按照式(2), 对地壳厚度的计算结果加以校正[20]。
VDSS方法通常测量每个台站的Tvdss, ABDD VDSS方法则测量台站对之间Tvdss的差值。按照实际Ss到时, 将台站i与j的记录对齐, 截取两个台站实际Ss到时之后相同时间窗内的记录进行互相关, 得到的到时差就是台站i与j间SsPmp与Ss到时差之差ΔTvdss,ij。ABDDVDSS方法就是将多个ΔTvdss,ij对应的方程进行联立, 求取最小二乘解, 得到每个台站到时差的参考解
。
因为求解的是到时差方程组, 所以给每个台站的参考解加上相同的偏移量dt后, 仍然是原方程组的解。本文取全部台站参考解总和为0s时的为ABDDVDSS法的相对到时解。当计算莫霍面的起伏形态时, 我们更关注相邻台站间的地壳厚度变化, 可以只使用相对到时解, 而不考虑偏移量 dt。当关注台站下方的莫霍面绝对深度时, 则需要求取偏移量 dt。选取同一地震事件下波形拟合法得到的每个台站的 Tvdss,i 和 ABDDVDSS 法得到的参考解, 通过式(3)求取偏移量, 使得两者残差的平方和最小, 则有
(3)
式(3)中, N为台站数量, 取
即为ABDDVDSS法的绝对到时解。本文中没有特别说明时, ABDDVDSS法测量得到的到时差均指绝对到时解。
根据表1所示的模型参数, 使用反射率法[26]计算介质响应, 选取4s的汉宁窗函数作为入射Ss子波, 模拟波形如图 2 所示。给模拟波形的垂向分量V和径向分量R分别加噪声, 作为5个台站的观测数据(图3)。加入的噪声是标准差为最大振幅0.3倍的白噪声在2～20 s频段带通滤波后的结果。分别采用波形拟合法和ABDDVDSS法测量Tvdss。本文中波形拟合法的正演波形通过反射率法计算得到。
在关于波形拟合方法的测试中, 给定的先验速度模型与表1一致, 通过变化地壳厚度来计算正演波形。截取正演波形在−10～20s内的数据, 与观测数据互相关, 记录互相关系数最大时的模型地壳厚度, 计算对应的 Tvdss 理论值和测量误差。每个含噪声数据的最佳拟合波形如图 3 所示。
表1 波形模拟参数
Table 1 Parameters for wave form simulation
参数 数值 地壳厚度(H)40 km 地壳平均P波速度(vp)6.3 km/s 地壳平均波速比(κ)1.73 地壳平均密度(ρ)2.78 g/cm3 射线参数(p)0.13 s/km 上地幔平均P波速度(vum)8.1 km/s 上地幔平均波速比(κum)1.73 上地幔平均密度(ρum)3.33 g/cm3
黑色点线为 Ss 子波, 是 4s 的汉宁窗函数; 黑色实线和虚线分别为垂向方向的介质响应和模拟波形, 灰色实线和虚线分别为径向的介质响应和模拟波形; 左、右黑色竖直点画线分别为 Ss 和 SsPmp 到时
图2 基于表 1 中参数模拟得到的波形
Fig. 2 Synthetic waveforms using parameters given in Table 1
在ABDDVDSS方法中, 选取Ss到时之后的4～14s时窗, 两两台站的波形进行互相关测量, 得到ΔTvdss,ij (表2), 然后联立对应的方程组, 求取每个台站最小二乘意义下的相对到时解
。再根据式(3)求取每个台站的绝对到时解Tvdss (表3)。
给正演波形的垂向分量V和径向分量R分别加50次如上所述的噪声, 作为50个台站的观测数据。用波形拟合和ABDDVDSS两种方法分别测量两个分量的Tvdss, 计算每个分量的Tvdss误差(图4 (a))和同一台站两个分量的Tvdss之差Tvdss,VR (图 4 (b))。
在不同震中距的记录中, 因为SsPmp相对于 Ss的相移不同, 所以互相关方法得到的到时差存在误差。ABDDVDSS法中, 需要使用台间距较小的两个台站的记录进行互相关来求取到时差, 将上述误差控制在合理的水平。
在30°～55°的震中距范围内, 假设地震深度为630km, 采用TauP程序[27‒28]计算理论射线参数。然后, 根据表1所示的壳幔模型参数, 使用反射率法计算介质响应。选取4s的汉宁窗函数作为入射Ss子波, 得到不同震中距下的模拟记录。将模拟记录中不同震中距下同一分量的数据进行两两互相关, 得到台站对间的ΔTvdss,ij。计算两个台站的真实Tvdss, 相减后与ΔTvdss,ij再次相减, 得到ΔTvdss,ij的测量误差, 结果如图5所示。
(a) 垂向分量; (b) 径向分量。灰色曲线为含噪数据, 黑色曲线为最佳拟合模型对应的模拟数据; 曲线左侧数值为Tvdss的误差
图3 含噪数据的波形拟合结果
Fig. 3 Wave fitting results of noise data
表2 不同台站对间的ΔTvdss,ij
Table 2 ΔTvdss,ij between different station pairs
台站ΔTvdss,ij/s 台站 1台站 2台站 3台站 4台站 5 台站 1‒0.100.020.120.00 台站 20.03‒0.000.00−0.13 台站 30.000.00‒0.00−0.07 台站 4−0.02−0.08−0.06‒−0.15 台站 50.00−0.04−0.020.01‒
说明: 每一格为台站号小的台站的Tvdss减去台站号大的台站的Tvdss, 右上方为垂向分量的结果, 左下方为径向分量的结果。
本文选取“中国地震科学台阵探测计划”二期台阵数据中与文献[15]相同的地震事件的记录, 分别使用波形拟合法和 ABDDVDSS 法测量到时差, 并估计研究区域的地壳厚度。震源位置为 123.15°E, 6.17°N, 震源深度为 632.4km, 矩震级为 6.6, 发震时刻为 2014 年 12 月 2 日 5 点 11 分 33 秒[29]。共有599 个台站记录到该事件, 震中距分布在 29°～43°之间。研究区域的地质构造分区如图 6 所示。
第一步, 估计每个台站的 Ss 波到时异常 tSs。将每个台站记录的垂向分量和径向分量进行偏振分析, 得到伪 S 分量和伪 P 分量。选取每个台站 Ss 波理论到时前后各 10s 内的伪 S 分量数据, 进行两两互相关, 记录 Ss 波的互相关到时差和互相关系数。选取互相关系数高于 0.8 的台站, 对其 Ss 波到时差进行方程组联立求解, 得到每个台站的 Ss 波到时异常值。
表3 波形拟合法与ABDDVDSS法的到时差误差
Table 3 Errors of Tvdss of wave fitting method and ABDDVDSS method
波形分量(测量方法)Tvdss的测量误差/s 台站1台站2台站3台站4台站5 V (ABDD)0.072−0.0240.002−0.0320.092 R (ABDD)−0.034−0.066−0.052−0.002−0.026 V (Fitting)0.036−0.1820.109−0.0730.219 R (Fitting)−0.0730.000−0.073−0.0360.000
说明: ABDD表示ABDDVDSS法的结果, Fitting表示波形拟合法的结果。
第二步, 对 Ss 波到时异常进行质量控制。选取互相关系数高于 0.5 和 0.7 的台站, 对 Ss 波到时差进行方程组联立求解, 得到另两个 Ss 波到时异常值的估计值。当同一台站的 3 个 Ss 波到时异常中任意两个之差值大于 1 s 时, 则对结果的约束不足, 便排除该台站的数据。对于剩余台站, 计算台站附近 1°范围内其他台站 Ss 波到时异常值的均值和标准差, 当该台站的 Ss 波到时异常在均值的两倍标准差以外时, 则可能是台站附近的异常结构或是台站的绝对时间不准, 即排除该台站的数据。进行质量控制后, 剩余 554 个台站数据的 Ss 波异常值分布如图 7所示。按照实际 Ss 波到时, 将剩余台站的伪 S 分量对齐, 取均值作为入射的 Ss 子波。
(a)单分量测量误差分布, 黑色正三角形和倒三角形分别表示 ABDDVDSS 法得到的垂向分量和径向分量 Tvdss 的误差, 灰色符号为波形拟合法相应的结果; (b) Tvdss,VR分布, 黑色正方形为 ABDDVDSS 法的结果, 灰色正方形表示波形拟合法的结果
图4 两种方法的测量误差和 Tvdss,VR 的分布
Fig. 4 Distribution of measurement error and Tvdss,VR for both methods
第三步, 用波形拟合法测量 Tvdss。给定与表 1一致的速度模型, 通过变化地壳厚度来计算介质的响应, 然后与估计的 Ss 子波进行卷积, 得到正演波形。截取正演波形在−10～20s 范围内的垂向分量和径向分量波形, 与观测数据对应的分量互相关, 记录互相关系数最大时的模型地壳厚度, 计算对应的到时差。取两个分量到时差的 Tvdss 均值作为测量结果, 同时计算两个分量的到时差之差 Tvdss,VR。图 8展示同一台站 Tvdss,VR 与两个分量的拟合结果中相对较小互相关系数的分布。
第四步, 对波形拟合得到的 Tvdss 进行质量控制。排除垂向分量或径向分量互相关系数小于 0.7和 Tvdss,VR大于 1s 的数据, 同时排除到时差在 1°以内的其他台站的到时差均值加减两倍标准差以外的数据。经过质量筛选后, 保留 264 个台站的测量结果。波形拟合法得到的 Tvdss 地理分布如图 9(a)所示, Tvdss,VR 地理分布如图 9(c)所示, Tvdss,VR 的概率分布如图 10(a)所示。
第五步, 用 ABDDVDSS 法测量到时差。截取每个台站 Ss 波实际到时之后 4～14s 时窗内的数据, 进行两两互相关, 记录台站对的距离、ΔTvdss,ij 和互相关系数。选取互相关系数高于 0.8 且台间距小于1°的 ΔTvdss,ij 进行方程组联立求解(垂向分量的方程总量为 2059, 径向分量的方程总量为 2202), 得到每个台站的相对到时解
。
第六步, 对相对到时解进行质量控制。选取互相关系数大于 0.5 和 0.7 的 ΔTvdss,ij, 按照与 Ss 波到时异常质量控制相同的步骤进行方程组联立求解, 得到另两个的估计值。当同一台站 3 个
中任意两个的差异大于 1s 时, 将该台站的数据排除。对于剩余台站, 排除到时差在 1°以内其他台站到时差均值加减两倍标准差以外的数据。结果质量筛选后, 保留 402 个台站的测量结果, 将垂向分量和径向分量的到时差取均值, 作为该台站的平均相对到时差。
第七步, 计算偏移量 dt。选取波形拟合法和ABDDVDSS 法都保留的台站到时差, 根据式(3)计算 ABDDVDSS 法参考解的整体偏移量 dt, 将其与平均相对到时差相加, 得到绝对到时解 Tvdss, 地理分布如图 9(b)所示。Tvdss,VR 的地理分布如图 9(d)所示, Tvdss,VR 的概率分布如图 10(b)所示, 同一台站Tvdss,VR 与方程数量的分布如图 11 所示, 波形拟合法和 ABDDVDSS 法得到的 Tvdss 的差异 D(Tvdss)如图10(c)所示。
横坐标和纵坐标分别为两个台站的震中距; 黑色实线为 1°台间距等值线, 黑色虚线为 2°台间距等值线; 右上方为 V 分量的误差, 左下方为R分量的误差
图5 不同震中距的记录的ΔTvdss,ij误差
Fig. 5 Error of ΔTvdss,ij between records of different epicentral distances
粗黑线表示地质块体分界线, 细黑线表示断层, 下同; 灰色三角形为台站位置。CAOB: 中亚造山带; HB: 河套盆地; AB: 阿拉善块体; OB: 鄂尔多斯块体; QOB: 祁连造山带; WQOB: 西秦岭造山带; SGB: 松潘甘孜块体; SCC: 华南克拉通
图6 研究区域的地质构造分区
Fig. 6 Geographic structural zones of the study region
(a)每个台站的 Ss 波到时异常, 圆点为台站位置, 下同; (b)线性插值后的 Ss 波到时异常分布
图7 Ss波到时异常分布
Fig. 7 Distribution of Ss arrival time anomaly
第八步, 根据 Ss 波到时异常和 ABDDVDSS 法得到的 Tvdss进行深度偏移。首先根据式(2), 假设ΔtSs = 0, 根据台站的 Tvdss 计算莫霍面反射点及地表转换点的位置和深度; 然后根据当前地表转换点的位置更新 ΔtSs, 根据式(2)重新计算莫霍面反射点及地表转换点的位置和深度; 重复上述过程, 直至收敛。校正 Ss 波到时异常后的地壳厚度分布如图 12 (b)和(d)所示。作为对比, 忽略 Ss 波的到时异常, 根据式(1)推导得出的地壳厚度分布如图 12(a)和(c)所示。
本文以测量得到的到时差与真值的误差作为测量精度的指标, 以垂向分量与径向分量的到时差之差 ΔTvdss,VR 作为测量准确度的指标, 以保留的台站数量作为空间分辨率指标。
虚线框示意Tvdss,VR小于 1 s, 互相关系数大于 0.7 的部分
图8 波形拟合法得到的同一台站 Tvdss,VR 与互相关系数的分布
Fig. 8 Distribution of Tvdss,VR and coherence in wavefitting method
在含噪模拟数据中, 波形拟合法得到的 Tvdss较为弥散, 容易在部分点位出现较大的测量误差(图4(a)和表 3), ABDDVDSS法得到的 ΔTvdss,ij 则更接近真值 0s (表 2), 整体上小于台站对用波形拟合方法分别得到的 Tvdss 的差值, 说明 ABDDVDSS 法得到的到时差比单台拟合到时具有更高的精度。同时, ΔTvdss,ij 相关方程的数量多于台站数量。冗余的高精度双差数据有效地降低了 ABDDVDSS 法得到的Tvdss 的误差(图 4(a))。在测量的准确度方面, 虽然两种方法对两个分量的 Tvdss 都是独立测量, 但是ABDDVDSS 法得到的 Tvdss,VR 具有更好的准确度(图4(b))。这是因为, 波形拟合法对单一分量的测量容易出现大的误差, 导致两个分量的差异比 ABDD VDSS 法大。
(a)波形拟合法得到的Tvdss; (b)ABDDVDSS法得到的Tvdss; (c)波形拟合法得到的Tvdss,VR; (d)ABDDVDSS法得到的Tvdss,VR
图9 Tvdss和 Tvdss,VR的分布
Fig. 9 Distribution of Tvdss and Tvdss,VR
(a)波形拟合法的Tvdss,VR; (b)ABDDVDSS法的Tvdss,VR; (c)波形拟合法与ABDDVDSS法的Tvdss之差D(Tvdss)
图10 波形拟合法与 ABDDVDSS 法 Tvdss 差异的概率分布
Fig. 10 Probability distribution of difference of Tvdss of wave fitting method and ABDDVDSS method
取每个台站两个分量中双差方程数量较小的
图11 同一台站Tvdss,VR与方程数量的分布
Fig. 11 Distribution of Tvdss,VR and number of relative equations on same station
在实际数据应用中, 波形拟合法得到的两个分量的互相关系数均超过 0.8 的台站仅有 80 个, 大部分台站(301 个)的互相关系数在 0.7～0.8 之间, 且准确度与互相关系数之间没有显著的关系(图 8)。在ABDDVDSS 法的测量中, 两个分量互相关系数大于 0.8 且台间距小于 1°的台站对数量分别是 2059 和 2202, 有 364 个台站至少有两个双差方程约束。同时, 台站相关的双差方程越多, 到时测量的准确度越高(图 11)。与波形拟合法相比, ABDDVDSS 法的 Tvdss,VR概率分布更向 0 s 集中(图 10(a)和(b))。
经过数据质量控制后, 波形拟合法保留了 264个台站的到时差结果, 对阿拉善块体北部和鄂尔多斯块体西南部覆盖度不佳; ABDDVDSS 法保留了402 个台站的到时差结果, 对各个块体均有良好的覆盖度(图 9)。在相同的台站, 两种方法的差异与其各自两个分量的差异 ΔTvdss,VR 相当(图 10(c))。ABDDVDSS 法能在保持精度的前提下, 获得更多高质量到时差结果, 提高了空间分辨率。
ABDDVDSS 法使用更多的更高精度的双差数据 ΔTvdss,ij 进行到时差求解, 结果中保留更多的更高精度的到时差数据 Tvdss, 准确度也随双差方程数量的增加而提高。ABDDVDSS 法有效地提升了结果的精度、准确度和空间分辨率。
ABDDVDSS 法要求台站对之间相移之差不能太大。台站间距过大时, 台站对间相移之差大, 互相关双差误差较大; 台间距过小时, 台站对数量少, 双差方程少, 即每个台站的约束少, ABDDVDSS 法精度下降。
ΔTvdss,ij 的误差近似地随着震中距的增加而增加, 当震中距接近全反射震中距时, 误差急剧增加(图 5)。随着震中距增加, 台站间距的上限降低, 可以有效地控制互相关双差误差。对于本文的实际数据, 震中距集中在 29°～43°之间, 设置 1°的最大台间距, ΔTvdss,ij 误差可以控制在 0.1s 以内, 同时台站对间的双差方程数量也较多。所以, 1°是合理的最大台间距取值。
(a)为不进行Ss波到时异常校正的地壳厚度, (b)为Ss波到时异常校正后的地壳厚度, (c)为(a)中点位线性插值后的地壳厚度, (d)为(c)中点位线性插值后的地壳厚度; 圆圈为莫霍面反射点的位置, 使用的到时为ABDDVDSS法得到的Tvdss
图12 地壳厚度分布
Fig. 12 Variation of crustal thickness
“中国地震科学台阵探测计划”二期台阵的布设范围(图 6)覆盖中亚造山带、祁连造山带、西秦岭造山带、松潘甘孜块体、阿拉善块体和鄂尔多斯块体等区域, 地形起伏巨大, 地质结构复杂。该区域Ss 波到时异常整体上为东北早、西南晚, 与地质块体及高程的分布吻合。中亚造山带和鄂尔多斯地块的 Ss 波实际到时早于理论到时, 祁连造山带、西秦岭造山带和松潘甘孜块体的 Ss 波实际到时晚于理论到时(图 7)。研究区内, Ss 波到时异常差异可达10 s, 不可忽视。
在鄂尔多斯块体及其周缘裂谷带, Ss 波到时异常校正前, 鄂尔多斯块体内部地壳厚度是中间大、南北小, 变化幅度较大。Ss 波到时异常校正之后, 鄂尔多斯块体核部地壳厚度减小, 周缘裂谷带地壳厚度增大, 块体内部地壳厚度变化幅度显著下降, 与接收函数法的结果[5‒6]更加一致。在祁连造山带与阿拉善地块交界处, 校正前地壳增厚区跨过分界线, 史克旭等[15]据此认为青藏高原东北缘的扩张前缘已到达阿拉善块体南缘。然而, 校正后阿拉善块体处地壳厚度为 50km, 越过与祁连造山带分界线后迅速增至 55km, 地壳加厚区止于两个块体的分界线, 表明青藏高原东北缘的扩张前缘仍处于祁连山北缘断裂处。西秦岭造山带内部, 校正前 100°E以西地壳厚度骤增, 校正后则 102.6°E 以西地壳厚度骤增, 校正后的分界线与地形(图 6)更加吻合。
鄂尔多斯块体的地壳呈现中间厚、南北薄的特点, 核心区域地壳横向变化小, 厚度为 45～47km, 环鄂尔多斯裂谷带地壳厚度下降至 39～43km。这与多条测线的人工震源反射探测结果[1‒4]一致, 两者的差异可能来自预设的 P 波速度的差异。接收函数法的结果[5‒7]显示, 鄂尔多斯块体内部地壳厚度南厚北薄, 以 36.5°N 为界, 南北厚度差异达 3～5km, 本文结果与之有较大的差异。Wang 等[5]指出, 在鄂尔多斯块体北部, 由于沉积层的影响, 接收函数波形出现异常, 接收函数法估计的地壳厚度比基于地形和布格重力异常估计的地壳厚度小。这或许是其差别的一种解释, 需要进一步证实。
阿拉善块体地壳厚度自东向西在 45～52km 之间缓慢增加, 在西南侧越过与祁连造山带的交界线后骤增至 55km。祁连造山带的地壳厚度较大, 自东向西厚度逐渐增加, 在 102.6°E 西侧地壳厚度骤增, 东侧为 46～54km, 西侧为 56～66km, 与基于地形和接收函数法估计的地壳厚度[5,7]一致。阿拉善块体与祁连造山带地壳厚度的变化与青藏高原的东北向挤压相关。
本文提出虚拟地震台阵双差测深法来测量SsPmp 与 Ss 波的到时差。该方法利用相邻台站间信号的相似性, 通过多道互相关测量台站间到时差的差, 然后联立双差方程组求解计算每个台站的到时差, 主要结论如下。
1)台间距保持不变时, 互相关到时差的误差随着震中距的增加而增加, 选用随着震中距增加而减小的最大台间距阈值, 可以控制互相关到时差的误差。
2)在含噪模拟数据和实际数据上的应用结果表明, 本文方法得到的双差到时误差小, 精度高, 能够降低地壳厚度估计误差。
3)如果忽略 Ss 波到时异常, 尤其是在地质背景复杂的区域, 会增加地壳厚度估计误差。本文的例子中, 忽略Ss波到时异常时, 会加大鄂尔多斯块体核部以及阿拉善块体西南部的地壳厚度, 减小祁连造山带中部的地壳厚度。
4)鄂尔多斯块体核部平均地壳厚度为 43～45 km, 横向变化小, 阿拉善地块地壳厚度为 45～52 km, 祁连造山带地壳厚度为 46～66 km, 地壳厚度东北方向的递减与青藏高原的东北向挤压相关。
致谢 地震数据来自“中国地震科学台阵探测计划”二期的观测资料, 研究工作得到北京大学高性能计算校级公共平台支持, 斯坦福大学博士研究生刘天泽、密歇根州立大学李嘉琪博士、南方科技大学俞春泉副教授和建信基金管理有限责任公司亢豆博士在论文完成过程中给予技术方面的指导, 在此一并表示感谢。
参考文献
[1]李松林, 张先康, 张成科, 等. 玛沁‒兰州‒靖边地震测深剖面地壳速度结构的初步研究. 地球物理学报, 2002, 45(2): 210–217
[2]滕吉文, 王夫运, 赵文智, 等. 阴山造山带‒鄂尔多斯盆地岩石圈层、块速度结构与深层动力过程. 地球物理学报, 2010, 53(1): 67–85
[3]滕吉文, 李松岭, 张永谦, 等. 秦岭造山带与邻域华北克拉通和扬子克拉通的壳、幔精细速度结构与深层过程. 地球物理学报, 2014, 57(10): 3154–3175
[4]Wang S, Wang F, Zhang J, et al. The P-wave velocity structure of the lithosphere of the North China Craton — results from the Wendeng-Alxa Left Banner deep seismic sounding profile. Science China Earth Scien-ces, 2014, 57(9): 2053–2063
[5]Wang C Y, Sandvol E, Zhu L, et al. Lateral variation of crustal structure in the Ordos block and surroun-ding regions, North China, and its tectonic implica-tions. Earth and Planetary Science Letters, 2014, 387: 198–211
[6]Wang W, Wu J, Fang L, et al. Sedimentary and crustal thicknesses and Poisson’s ratios for the NE Tibetan Plateau and its adjacent regions based on dense seismic arrays — science direct. Earth and Planetary Science Letters, 2017, 462: 76–85
[7]Xu X, Niu F, Ding Z, et al. Complicated crustal deformation beneath the NE margin of the Tibetan plateau and its adjacent areas revealed by multista-tion receiver-function gathering. Earth and Planetary Science Letters, 2018, 497: 204–216
[8]Zhu L, Kanamori H. Moho depth variation in southern California from teleseismic RFs. Journal of Geophy-sical Research Atmospheres, 2000, 105(B2): 2969–2980
[9]危自根, 储日升, 陈凌, 等. 复杂地壳接收函数 H-κ叠加——以安纳托利亚板块为例. 地球物理学报, 2016, 59(11): 4048–4062
[10]Luo S, Zhu L, Huang R, et al. Determination of crus-tal thickness and velocities by using receiver func-tions and PmP travel times. Geophysical Journal In-ternational, 2019, 216(2): 1304–1312
[11]Tseng T L, Chen W P, Nowack R L. Northward thinning of Tibetan crust revealed by virtual seismic profiles. Geophysical Research Letters, 2009, 36(24): doi: 10.1029/2009GL040457
[12]Tian X, Chen Y, Tseng T L, et al. Weakly coupled lithospheric extension in southern Tibet. Earth and Planetary Science Letters, 2015, 430: 171–177
[13]亢豆, 俞春泉, 陈九辉, 等. 川滇地区地壳结构的虚拟地表震源反射测深法研究. 北京大学学报(自然科学版), 2017, 53(5): 825–832
[14]魏晓拙, 姜明明, 陈凌, 等. 一种基于密集线性台阵的虚拟地震测深新方法及其应用. 地球物理学进展, 2018, 33(3): 986–992
[15]史克旭, 张瑞青, 肖勇. 利用虚拟地震测深法约束青藏高原东北缘及周边地区地壳厚度. 地球物理学报, 2020, 63(12): 4369–4381
[16]Chen W P, Jiang Y. Undulating Moho beneath a near-uniform surface of central Tibet. Earth and Planetary Science Letters, 2020, 543: 116343
[17]Yu C Q, Chen W P, Ning J Y, et al. Thick crust beneath the Ordos plateau: implications for instability of the North China craton. Earth & Planetary Science Letters, 2012, 357/358: 366–375
[18]Liu T, Klemperer S L, Yu C, et al. Post-critical SsPmp and its applications to virtual deep seismic sounding (VDSS) — 1: sensitivity to lithospheric 1-D and 2-D structure. Geophysical Journal International, 2018, 215(2): 880–894
[19]Liu T, Klemperer S L, Yu C, et al. Post-critical SsPmp and its applications to virtual deep seismic sounding (VDSS) — 3: back-projection imaging of the crust–mantle boundary in a heterogeneous lithosphere, the-ory and application. Geophysical Journal Internatio-nal, 2020, 223(3): 2166–2187
[20]Yu C Q, Chen W P, van der Hilst R D. Constraints on residual topography and crustal properties in the wes-tern United States from virtual deep seismic soun-ding. Journal of Geophysical Research, 2016, 121: 5917‒5930
[21]Parker Jr E H, Hawman R B, Fischer K M. Estima-ting crustal thickness using SsPmp in regions covered by low-velocity sediments: imaging the Moho beneath the Southeastern Suture of the Appalachian margin experiment (SESAME) array, SE Atlantic Coastal Plain. Geophysical Research Letters, 2016, 43(18): 9627–9635
[22]Kang D, Yu C Q, Ning J Y, et al. Simultaneous deter-mination of crustal thickness andpwavespeed by vir-tual deep seismic sounding (VDSS). Seismological Re-search Letters, 2016, 87(5): 1104–1111
[23]Thompson D, Rawlinson N, Tkalčić H. Testing the limits of virtual deep seismic sounding via new crus-tal thickness estimates of the Australian continent. Geophysical Journal International, 2019, 218: 787–800
[24]Yu C Q, Chen W P, van der Hilst R D. Removing source-side scattering for virtual deep seismic soun-ding (VDSS). Geophysical Journal International, 2013, 195: 1932‒1941
[25]刘震, 田小波, 朱高华, 等. SsPmp震相地壳探测方法. 地球物理学报, 2015, 58(10): 3571–3582
[26]Randall G E. Efficient calculation of differential seis-mograms for lithospheric receiver functions. Geophy-sical Journal International, 1989, 99(3): 469‒481
[27]Buland R, Chapman C H. The computation of seismic travel times. Bulletin of the Seismological Society of America, 1983, 73(5): 1271–1302
[28]Crotwell H P, Owens T J, Ritsema J. The TauP toolkit: flexible seismic travel-time and ray-path utilities. Seismological Research Letters, 1999, 70(2): 154–160
[29]BondÁR I, Storchak D. Improved location procedures at the international seismological centre. Geophysical Journal International, 2011, 186(3): 1220‒1244
Array-Based Double Difference Virtual Deep Seismic Sounding Method and Its Application
Abstract Through multi-channel cross-correlation, the differences of arrival time differences of SsPmp and Ss between adjacent stations in dense broadband seismic array are measured, and then arrival time differences of SsPmp and Ss of each station are obtained by means of least squares. This method is an advanced virtual deep seismic sounding (VDSS) method and named Array-Based Double Difference Virtual Deep Seismic Sounding (ABDDVDSS) method, which effectively reduces arrival time difference error and estimates crustal thickness precisely by taking advantages of array data. This method is applied to synthetic data and ChinArray II data. Results show that the proposed method can significantly improve the precision of measurement.
Key words virtual deep seismic sounding (VDSS); array; multi-channel cross-correlation; double difference; Ss arrival time anomaly
doi: 10.13209/j.0479-8023.2022.044
国家自然科学基金(41874071)和地震行业科研专项经费(201408013)资助
收稿日期: 2021-05-28;
修回日期: 2021-06-25