CN115755183A - 一种基于尾波干涉的震源位置追踪方法 - Google Patents

Abstract

本发明的目的在于提供一种基于尾波干涉的震源位置追踪方法，包括如下步骤：获取震源位置变化前、后的时域信号；信号分段；震源位置变化的时延估计；距离差计算；震源位置估计；震源位置追踪。从理论角度出发，本发明结合尾波特性，利用尾波对介质变化敏感，远距离接收下尾波信号能量强等特点，从尾波相位变化展开分析，区分介质变化和震源位置变化的时延，尤其是非均匀介质变化时，去除介质变化导致的时延，准确可靠地求取震源位置变化导致的时延，实现了高精度的震源位置追踪；从工程应用角度，本发明对于台站分布要求更低，对介质速度结构的依赖性不强，实用性更高。

Description

一种基于尾波干涉的震源位置追踪方法

技术领域

本发明涉及的是一种地震探测方法，具体地说是震源位置追踪方法。

背景技术

中小地震的震源位置及其变化对于研究地震构造、地震活动性、断层性质以及地震应力场的演化等问题都有极为重要的意义。常用的震源位置追踪方法主要是基于直达波的处理应用，直达波虽然能够估计地震事件间时延，但无法确定时延是由震源位置变化引起，还是由介质变化引起。同时追踪结果的精度依赖于台站分布和介质速度结构，当对于介质速度精细结构了解不足时，不同台站间接收信号的路径不同，可能存在直达波不是最快到达的，将导致基于直达波的震源位置追踪方法效果欠佳。在复杂地质环境下，传播路径或介质速度变化发生微弱扰动时，地震事件间直达波部分保持高度一致，很难进行区分。同时直达波信号会随着散射转换为尾波信号，其幅值会快速衰减，尾波能量逐渐增强，甚至在远程传播的场景下，直达波的能量可能会弱于尾波能量。例如，文献《Time—lapseMonitoring with Coda Wave Interferometry》和《Monitoring Velocity Variations inthe Crust Using Earthquake Doublets'An Application to the Calaveras Fault,California》中存在台站接收信号的尾波能量强于直达波。因此，现有的震源位置追踪方法难以区分震源位置变化及介质结构变化，同时难以在远距离接收下实现对震源位置变化的高精度追踪。

综上所述，现有震源位置追踪方法难以实现远距离接收，震源位置扰动以及存在非均匀速度变化场时的震源位置高精度追踪，无法满足地震构造等研究的高精度需求。

发明内容

本发明的目的在于提供能满足地震活动性、断层性质等研究需求的一种基于尾波干涉的震源位置追踪方法。

本发明的目的是这样实现的：

本发明一种基于尾波干涉的震源位置追踪方法，其特征是：

(1)获取震源位置变化前、后的时域信号；

(2)信号分段；

(3)震源位置变化的时延估计；

(4)距离差计算；

(5)震源位置估计；

(6)震源位置追踪。

本发明还可以包括：

1、步骤(1)获取震源位置变化前、后的时域信号具体为：

确定初始震源位置S₁(x₁,y₁)、接收点位置R_i＝(X_i,Y_i),1≤i≤n及介质传播速度v，获取接收点R_i处震源位置变化前、后的时域信号A_i(t)和B_i(t)。

2、步骤(2)信号分段具体方法为：

基于总时间窗口t_all、移动时间窗的窗长t_mw、移动步距t_step，获取分段信号A_i,j(t),1≤j≤k和B_i,j(t)，其中k指分段个数。

3、步骤(3)震源位置变化的时延估计具体为：

基于尾波干涉法获取分段信号间偏移量α_i,j，由偏移量α_i,j和流逝时间t_j拟合，确定多项式方程，对多项式方程做差分与反差分运算，求取速度变化导致的时延β_i,j，去除该时延，获取震源位置变化导致的时延τ_i：

a、计算偏移量：基于尾波干涉法计算接收点R_i的震源位置变化前、后的分段信号A_i,j(t)和B_i,j(t)间的偏移量α_i,j；

b、速度变化的时延估计：对偏移量α_i,j和流逝时间t_j做多项式拟合，由曲线变化趋势自行确定多项式阶数n，基于残差最小原则确定多项式系数a_n,a_n-1,...,a₀，以确定多项式方程；对多项式方程做差分与反差分运算，求取速度变化导致的时延β_i,j，将β_i,j与t_j做线性拟合，其斜率为接收点R_i处估计得到的速度扰动Δv_i；

c、求取震源位置改变的时延：偏移量α_i,j减去时延β_i,j，记为τ_i,j，再对τ_i,j取平均，即为接收点R_i处信号A_i(t)和B_i(t)间震源位置变化导致的时延τ_i。

4、步骤(4)距离差计算具体为：

基于步骤(3)所得的时延τ_i与传播速度v相乘获得震源点S₁与震源点S₂到接收点R_i距离差Δr_i。

5、步骤(5)震源位置估计具体为：

基于步骤(4)求取的距离差Δr_i，获取未知坐标震源点S₂到接收点R_i的距离，结合接收点坐标信息，使用双曲线定位，估计出震源点S₂的坐标位置；

A、震源位置变化后与接收点的距离估计：首先计算初始震源位置S₁(x₁,y₁)到接收点R_i＝(X_i,Y_i),i≥1的距离d_i，结合距离差Δr_i，获得震源点S₂到接收点R_i的距离D_i；

B、双曲线定位：分别以接收点R_i和R_j,j≠i为圆心，D_i和D_j为半径，绘制双曲线并求取双曲线焦点，选择靠近初始震源点S₁的焦点作为震源点S₂由接收点R_i和R_j估计的一个坐标位置S_i,j,i≠j；

C、位置平均：任意两个接收点估计一个震源位置S_i,j,i≠j，结合S_i,j和初始震源点S₁，将十分偏离S₁或与其余S_i,j的坐标值剔除掉，对剩余的估计震源位置取平均，即为震源S₂的坐标位置。

6、步骤(6)震源位置追踪具体为：

基于速度扰动Δv_i,1≤i≤n取均值获得平均速度扰动Δv，再由Δv补偿介质传播速度v；再将步骤5定位出的震源S₂的位置作为初始震源点，重复步骤(1)到步骤(5)定位出下一个震源移动位置；以此类推，实现震源位置追踪。

本发明的优势在于：

1、从理论角度出发，本发明结合尾波特性，利用尾波对介质变化敏感，远距离接收下尾波信号能量强等特点，从尾波相位变化展开分析，区分介质变化和震源位置变化的时延，尤其是非均匀介质变化时，去除介质变化导致的时延，准确可靠地求取震源位置变化导致的时延，实现了高精度的震源位置追踪；

2、从工程应用角度，本发明对于台站分布要求更低，对介质速度结构的依赖性不强，实用性更高。

附图说明

图1为本发明的流程图；

图2为本发明的强散射介质速度模型图；

图3是本发明的弱散射介质速度模型图；

图4是本发明的震源位置移动图；

图5是本发明的强散射介质接收点R1处信号图；

图6是本发明的弱散射介质接收点R1处信号图；

图7是本发明的非线性时延估计结果图；

图8是本发明的追踪结果误差棒图；

图9是本发明的累积分布函数CDF图。

具体实施方式

下面结合附图举例对本发明做更详细地描述：

结合图1-9，本发明具体包括如下步骤：

步骤1获取震源位置变化前、后的时域信号的详细过程是本发明选择构建仿真信号，首先构建散射介质速度模型，如附图2和附图3所示，确定介质传播速度v，设置初始震源点位置S₁(x₁,y₁)和接收点R_i＝(X_i,Y_i),i≥1坐标位置。然后激发震源，接收点R_i接收记录震源点位于S₁处的时域信号A_i(t)。随后，改变震源点位置至S₂(x,y)，激发震源，接收点R_i接收记录震源点位于S₂的时域信号B_i(t)；

步骤2信号分段：信号分段用于估计局部时延，观察其局部时延随流逝时间的变化趋势，用于后续区分时延变化的产生原因。首先确定总时间窗口t_all，确保总时间窗口内尾波信号占主导，且信噪比较高，并对接收点R_i处时域信号A_i(t)和B_i(t)重新赋值，即A_i(t)＝A_i(t_all)和B_i(t)＝B_i(t_all)，使其只包含总时间窗口内的信号。设置移动时间窗的窗长t_mw，移动时间窗指具有指定长度t_mw的时间窗口，包含在总时间窗口内，按照设定的移动步距t_step沿着时间轴滑动该移动窗，实现信号分段。需要注意的是移动时间窗口内需要包含足够多的信息量，即包含的波形要足够复杂，但移动时间窗的窗长过大时，计算速度将会降低，因此需要合理的设置窗长。同时，为了保证分析结果的连续性，移动步距需要较小，使移动时间窗口内信号有重叠。假定总时间窗口t_all＝[t₁-t_w,t_k+t_w]，移动时间窗的窗长为t_mw＝2t_w，其中滑动移动时间窗口后该窗口的中心位置为t_j,1≤j≤k，其中k为分段信号个数，本发明使移动时间窗内包含五个周期的波形变化，移动步距设置为窗长的1/15，即t_step＝|t_j-t_j-1|＝1/15×t_mw，获取分段信号A_i,j(t),1≤j≤k和B_i,j(t),1≤j≤k；

步骤3震源位置改变的时延估计过程包括：

3.1：计算偏移量。震源位置变化或介质变化时，将导致信号间偏移量呈现长周期的变化趋势，因此可通过获取信号间偏移量来检测介质变化，震源位置变化，或同时存在介质变化和震源位置变化。以Windowed Cross-Correlation(WCC)方法为例，属于一种尾波干涉法，计算接收点R_i震源位置变化前、后分段信号A_i,j(t)和B_i,j(t)间的偏移量。假定时延变量t_s在[-2t_w,2t_w]范围内，求取接收点R_i第j段分段信号A_i,j(t)和B_i,j(t)，t∈[t_j-t_w,t_j+t_w]的互相关系数R(t_s)，公式如下所示，互相关系数最大值所对应的时延变量t_s，即获取分段信号A_i,j(t)和B_i,j(t)的偏移量α_i,j,1≤j≤k，其中k指分段信号个数；

3.2：速度变化的时延估计。均匀速度变化使信号间呈现线性时延，而非均匀速度变化导致非线性时延，绘制根据步骤3.1计算得到的偏移量α_i,j随流逝时间t_j的变化曲线，由曲线变化趋势自行确定多项式阶数，推荐设置2阶或以上，避免总时间窗口靠后时，时延可能呈现线性变化而实际可能更倾向于非线性变化。通过对偏移量α_i,j和流逝时间t_j做多项式拟合，根据残差最小的原则确定合适的多项式系数a_n到a₀，即确定多项式方程y＝a_nxⁿ+a_n-1x^n-1+...+a₁x+a₀，其中x∈[1,t_k]，t_k表示流逝时间t_j的最后一个值，获取y变化值。基于变量y做差分运算，指后一个数减去前一个数，即X(m)＝y(m+1)-y(m)，且m∈[1,t_k-1]，X变量即为差分运算的结果。然后对X做反差分运算，Z(n)＝X(n)+Z(n-1)，n∈[1,t_k-1]且设置Z(0)＝0，Z变量表示反差分运算结果。最后，求取非均匀速度变化导致的时延β_i,j，β_i,j＝Z(t_j-1),1≤j≤k；将β_i,j与t_j做线性拟合，其斜率为接收点R_i处估计得到的速度扰动Δv_i；

3.3：求取震源位置改变的时延。将偏移量α_i,j减去时延β_i,j，记为τ_i,j，即去除速度变化导致的时延。再对τ_i,j取平均，即获取接收点R_i信号A_i(t)和B_i(t)间震源位置变化导致的时延τ_i，τ_i＝(τ_i,1+...+τ_i,k)/k；

步骤4距离差计算的过程是对每个接受点R_i计算得到的时延τ_i转换为震源位置变化前、后到接收点R_i的距离差Δr_i，Δr_i＝τ_i×v，v指介质传播速度；

步骤5震源位置估计的详细过程是：

5.1：震源位置变化后与接收点的距离估计。初始震源点S₁和接收点R_i的坐标位置已知，根据公式

求取初始震源点S₁到每个接收点R_i的距离d_i。再根据步骤4计算得到的距离差Δr_i，获得震源点S₂到接收点R_i的距离D_i，D_i＝d_i-Δr_i；

5.2：双曲线定位。从步骤5.1获得震源点S₂到接收点R_i的距离D_i，虽然获得S₂到接收点R_i的距离，但其S₂的坐标可以是以R_i为圆心，D_i为半径的弧线上任意点。因此，选择以接收点R_i为圆心，D_i为半径；再以接收点R_j,j≠i为圆心，D_j为半径，通过双曲线求焦点，即联立方程组，求解坐标S₂(x,y)。方程组如下所示，理论上存在两个解，即S₂的坐标存在两个，但由于尾波干涉法计算相对震源位置时，需要满足震源间距离小于一个波长，否则存在周期跳跃。因此，选择靠近初始震源点S₁的焦点作为震源点S₂由接收点R_i和R_j估计的一个坐标位置S_i,j,i≠j；

5.3：位置平均。每两个接收点能估计一个震源位置S_i,j,i≠j，结合S_i,j和初始震源点S₁，将十分偏离S₁或与其余S_i,j的坐标值剔除掉，对剩余的估计震源位置取平均，即为震源S₂的坐标位置；

步骤6对所有接收点获取的速度扰动Δv_i,1≤i≤n取均值，获取平均速度扰动Δv＝(Δv₁+,...,+Δv_n)/n，对介质传播速度v进行补偿，即v＝v+Δv。将步骤5定位出的震源S₂的位置作为初始震源点，重复步骤1到步骤5定位出下一个震源移动位置。以此类推，实现震源位置追踪。

下面通过仿真实例来对其进行进一步说明：

首先构建散射介质速度模型图，如图2和图3所示，其中图2为强散射介质速度模型，图3为弱散射介质速度模型，初始震源点S₁坐标为(1606m,1600m)，S₂坐标为(1606m,1602m)，逆时针近似圆形移动，如附图4所示，选择3个接收点R₁(3958m,3900m)，R₂(3958m,3958m)和R₃(1978m,3958m)接收记录信号。随后，使用Seismic Unix软件中SUFDMOD2模块获取仿真信号，本仿真的速度场边界设置为吸收边界，使用主频为200Hz的雷克子波作为震源，传播时间为5s，时间采样率为4000点，且介质传播速度v为速度模型取均值后的速度，例如，图2中强散射介质模型求取的速度v＝3322m/s。图5中A₁(t)和图6分别是接收点R₁在强散射介质和弱散射介质下接收到的初始震源点S₁处时域信号，图6中直达波能量强，持续时间短，信号能量在10^-2级别，而图5中直达波不明显，且信号能量在10^-3级别，说明在强散射介质下直达波不明显且信号整体能量偏弱，与弱散射介质中尾波的能量相近，也反映强散射介质下可能尾波信号占主导，传统基于直达波的定位方法将不适用。同样的方式在接收点R₁接收到的震源点S₂的时域信号，并对该信号进行非线性拉伸，构造既有非均匀速度变化又存在震源位置改变的仿真信号B₁(t)，如图5中B₁(t)所示，以获取可控的介质变化。

步骤2中以R₁接收到的信号A₁(t)为例，总时间窗口t_all设置为[3.75s，4.25s]，对A₁(t)重新赋值，即A₁(t)＝A₁(t_all)。为了方便处理，将时间以时间点数表示，转换关系为以秒为单位的时间乘以时间采样率，即t_all为[15001，17000]，移动时间窗的窗长t_mw为150，移动步距t_step设置为10，则分段信号个数k＝184，获得分段信号A_1,1,A_1,2,...,A_1,184，其中A_1,1＝A1(15001:15150)，A_1,2＝A₁(15151:15300)，以此类推。同理，参数设置一致的情况下，由B₁(t)信号也将获得分段信号B_1,1,B_1,2,...,B_1,184。以同样的方式可以处理接收点R₂、R₃的时域信号A₂(t)和B₂(t)以及A₃(t)和B₃(t)，获取分段信号。

由步骤3.1处理信号A_1,1(t)和B_1,1(t)，获取时延α_1,1＝-2.62，同理处理A_1,2(t)和B_1,2(t)直至A_1,184(t)和B_1,184(t)获取接收点R₁处的水平偏移量α_1,2,α_1,2,...,α_1,184。同样的方式处理接收点R2、R3的分段信号，获得不同接收点下的水平偏移量α_2,1,α_2,2,...和α_3,1,α_3,2,...。绘制α_1,1,α_1,2,...,α_1,184随流逝时间t_j,1≤j≤k的变化曲线，其中t_j为移动时间窗的中心位置所对应的时间，例如t₁＝(15001+15150)/2＝15076，如附图7中点划线α_1,j所示，选择二阶多项式对曲线进行拟合，确定多项式方程y＝-1.695e^-08x²+1.077e^-05x+0.7598，其中x∈[1,16906]。对y做差分运算，例如X(1)＝y(2)-y(1)＝1.0722e^-05，依次类推，获取X(2),X(3),...,X(16095)。然后依据X获取Z值，例如Z(1)＝X(1)+Z(0)＝1.0722e-05，Z(2)＝Z(1)+X(2)＝2.1411e^-05，其中Z(0)＝0，以此类推，直到获得Z(16095)。然后，对Z取与流逝时间t_j对应的值即为非均匀速度变化导致的时延β_1,j，β_1,j＝Z(t_j-1),1≤j≤k，如图7中黑色实线所示，发现与理论的非均匀速度变化的变化趋势相近，能较好地估计出非线性时延，但存在一定的误差。对β_1,j与流逝时间t_j做线性拟合，其斜率即为即速度扰动Δv₁＝-0.0005312。将分段信号求取的α_1,j与β_1,j相减，并取均值，即为接收点R₁估计的震源位置变化导致的时延，即τ₁＝(α_1,1-β_1,1)+,...+(α_1,k-β_1,k)/k。同理，估计接收点R₂、R₃的震源位置变化导致的时延τ₂和τ₃，以及速度扰动Δv₂、Δv₃。

由公式Δr_i＝τ_i×v求得距离差Δr₁、Δr₂以及Δr₃。由

计算初始震源点S₁到接收点R₁、R₂以及R₃的距离d₁、d₂和d₃。再由公式D_i＝d_i-Δr_i，估计震源点S₂到各接收点R_i的距离，具体值如表1所示。结合接收点位置，联立方程组，估计震源点S₂的坐标位置。接收点R₁和R₂双曲线定位的位置是S_1,2＝(3490.2402,645.1071)，同理，得到S_1,3＝(1605.7274,1602.1878)以及S_2,3＝(1605.7791,1601.8645)，其中S_1,3和S_2,3的定位结果与S₁位置相近，而S_1,2的结果与S₁位置相差甚远，因此，不考虑S_1,2的位置，对S_1,3和S_2,3的结果取均值，作为S₂的估计值，即S₂＝(1605.7533,1602.0261)与理论值S₂＝(1606,1602)基本吻合。

将所有接收点R_i估计得到的速度扰动Δv_i,1≤i≤3取均值，获取平均速度扰动Δv＝(Δv₁+Δv₂+Δv₃)/3，将介质传播速度进行补偿，即v＝v+Δv，以S₂的估计坐标信息作为初始震源点，估计S₃的坐标位置，依次类推，实现源点追踪。其定位结果绘制误差棒图，如图8所示，端线的长短表示误差的大小，上、下、左、右朝向分别表示Z坐标估计偏大，Z坐标估计偏小，X坐标估计偏小以及X坐标估计偏大，整体来看X坐标误差小，Z坐标的误差明显，但也能够较为准确地估计出震源位置变化的轨迹。为了更直观地描述定位误差，选择累积分布函数(Cumulative Distribution Function,CDF)评价指标来衡量，其中累积分布函数CDF如图9所示。累积分布函数CDF横坐标为定位测量误差，纵坐标是累积分布函数，表示在某个定位精度下的定位次数在总定位次数中的百分比，WCC方法的最大位置误差，即概率为1时对应的定位误差为1.58m。总的来说，定位结果误差不大，说明利用尾波干涉法能够高精度地实现震源位置追踪。

表1估计参数

Claims (7)

1.一种基于尾波干涉的震源位置追踪方法，其特征是：

(1)获取震源位置变化前、后的时域信号；

(2)信号分段；

(3)震源位置变化的时延估计；

(4)距离差计算；

(5)震源位置估计；

(6)震源位置追踪。

2.根据权利要求1所述的一种基于尾波干涉的震源位置追踪方法，其特征是：步骤(1)获取震源位置变化前、后的时域信号具体为：

确定初始震源位置S₁(x₁,y₁)、接收点位置R_i＝(X_i,Y_i),1≤i≤n及介质传播速度v，获取接收点R_i处震源位置变化前、后的时域信号A_i(t)和B_i(t)。

3.根据权利要求1所述的一种基于尾波干涉的震源位置追踪方法，其特征是：步骤(2)信号分段具体方法为：

基于总时间窗口t_all、移动时间窗的窗长t_mw、移动步距t_step，获取分段信号A_i,j(t),1≤j≤k和B_i,j(t)，其中k指分段个数。

4.根据权利要求1所述的一种基于尾波干涉的震源位置追踪方法，其特征是：步骤(3)震源位置变化的时延估计具体为：

基于尾波干涉法获取分段信号间偏移量α_i,j，由偏移量α_i,j和流逝时间t_j拟合，确定多项式方程，对多项式方程做差分与反差分运算，求取速度变化导致的时延β_i,j，去除该时延，获取震源位置变化导致的时延τ_i：

a、计算偏移量：基于尾波干涉法计算接收点R_i的震源位置变化前、后的分段信号A_i,j(t)和B_i,j(t)间的偏移量α_i,j；

b、速度变化的时延估计：对偏移量α_i,j和流逝时间t_j做多项式拟合，由曲线变化趋势自行确定多项式阶数n，基于残差最小原则确定多项式系数a_n,a_n-1,...,a₀，以确定多项式方程；对多项式方程做差分与反差分运算，求取速度变化导致的时延β_i,j，将β_i,j与t_j做线性拟合，其斜率为接收点R_i处估计得到的速度扰动Δv_i；

c、求取震源位置改变的时延：偏移量α_i,j减去时延β_i,j，记为τ_i,j，再对τ_i,j取平均，即为接收点R_i处信号A_i(t)和B_i(t)间震源位置变化导致的时延τ_i。

5.根据权利要求1所述的一种基于尾波干涉的震源位置追踪方法，其特征是：步骤(4)距离差计算具体为：

基于步骤(3)所得的时延τ_i与传播速度v相乘获得震源点S₁与震源点S₂到接收点R_i距离差Δr_i。

6.根据权利要求1所述的一种基于尾波干涉的震源位置追踪方法，其特征是：步骤(5)震源位置估计具体为：

基于步骤(4)求取的距离差Δr_i，获取未知坐标震源点S₂到接收点R_i的距离，结合接收点坐标信息，使用双曲线定位，估计出震源点S₂的坐标位置；

A、震源位置变化后与接收点的距离估计：首先计算初始震源位置S₁(x₁,y₁)到接收点R_i＝(X_i,Y_i),i≥1的距离d_i，结合距离差Δr_i，获得震源点S₂到接收点R_i的距离D_i；

B、双曲线定位：分别以接收点R_i和R_j,j≠i为圆心，D_i和D_j为半径，绘制双曲线并求取双曲线焦点，选择靠近初始震源点S₁的焦点作为震源点S₂由接收点R_i和R_j估计的一个坐标位置S_i,j,i≠j；

C、位置平均：任意两个接收点估计一个震源位置S_i,j,i≠j，结合S_i,j和初始震源点S₁，将十分偏离S₁或与其余S_i,j的坐标值剔除掉，对剩余的估计震源位置取平均，即为震源S₂的坐标位置。

7.根据权利要求1所述的一种基于尾波干涉的震源位置追踪方法，其特征是：步骤(6)震源位置追踪具体为：

基于速度扰动Δv_i,1≤i≤n取均值获得平均速度扰动Δv，再由Δv补偿介质传播速度v；再将步骤5定位出的震源S₂的位置作为初始震源点，重复步骤(1)到步骤(5)定位出下一个震源移动位置；以此类推，实现震源位置追踪。

Images