一种高分辨率确定地震事件P波反方位角和慢度的方法与流程

本发明属于地震事件p波处理领域，涉及一种地震台阵数据处理方法，具体是一种高分辨率确定地震事件p波反方位角和慢度的方法。

背景技术：

地震事件中的p波即为纵波，p波为从震源发出，一路不断折射传播过来的波，它是远震测量中的首波。

通常确定地震事件的p波反方位角和慢度参数的方法有频率域聚束(fk)和时间域聚束(beamforming)。对于频率域聚束(fk)或时间域聚束(beamforming)的方法确定地震事件的p波反方位角和慢度参数的分辨率低，不便于工作人员分析、察看，且通过fk或beamforming的方法不能够区分多个地震事件的p波反方位角和慢度参数，为解决上述缺陷，现提供一种解决方案。

技术实现要素：

本发明的目的是为了解决对于频率域聚束(fk)或时间域聚束(beamforming)的方法确定地震事件的p波反方位角和慢度参数的分辨率低，不便于工作人员分析、察看，且通过fk或beamforming的方法不能够区分多个地震事件的p波反方位角和慢度参数的问题，而提出一种高分辨率确定地震事件p波反方位角和慢度的方法。

本发明的目的可以通过以下技术方案实现；

一种高分辨率确定地震事件p波反方位角和慢度的方法，该方法具体包括以下步骤：

步骤一：设台阵由n个台组成，取离台阵中心最近的台作为参考台，设为原点；

步骤二：读取n个台的垂直分量地震位移波形记录，并对记录进行去仪器响应，去均值处理，重采样到相同采样率；

步骤三：对步骤二中的所有台的波形记录按发阵时刻对齐，选取时间窗口；

步骤四：对选取的地震记录加窗后进行傅立叶变换，得到记录谱即观测向量d(f)；

步骤五：将慢度和反方位角参数化，选取慢度以及反方位角范围进行均匀网格划分；

步骤六：根据子台位置，以及慢度和反方位角参数化的信息，构建一个测量矩阵g(f)；

步骤七：通过测量矩阵g(f)以及步骤四得到的观测向量d(f)，应用正交匹配追踪方法求出解向量xsol；

步骤八：通过查寻解向量xsol中非零元素对应的索引值，反推到慢度-反方位角空间中，得到最终的解。

进一步地，步骤二中所述时间窗口将每道的p波包含在内。

进一步地，该方法的工作原理如下：

地震远离台阵的情况下，假设地震波为平面波的情况下，以反方位角φ和慢度s(s/km)到达台阵下方。假设2维地震台阵由n子台构成，这里不考虑子台间的垂直方向的高程差。以台阵的中心最近的台站为参考台为坐标原点建立笛卡尔坐标系，第i个子台的位置向量为ri＝[xiyi]。对于所有子台接收到的数据定义为d＝[d1(t)，…，di(t)，…，dn(t)]矩阵，矩阵大小为ns×n，其中ns是在时间维度的点数，di(t)表示为第i个子台接收到的时间序列，为一向量，长度为ns。

在频率域中，对于给定的频率f0，聚束能量e可以用公式(1)计算：

其中j表示虚部单位,di(f0)为数据di(t)在f0的傅立叶变换，相移项中的τi表示在第i个台站的与原点的到时差，可以用公式(2)计算(schweitzeretal.2002),

τi＝-(xisinφ+yicosφ)s(2)

在实际情况中，给定慢度和反方位角的搜寻范围，我们将慢度和反方位角网格划分为m＝sslw×bbaz,其中sslw和bbaz分别表示慢度网格点数和反方位角网格点数，对于fk方法，我们可以搜寻网格点来使得聚束能量最大，该网格点对应的慢度和反方位角为我们的解。但这种方法给出的结果具有低分辨率的特点，并且不能够轻易的区分多个源混合的现象。

假设每隔网格点为一个激发源，具备相应的慢度和反方位角。那么这个问题就转化为在慢度-反方位角空间定位的问题(yaoetal.,2011)，每隔网格点的激发源的谱可以简单表示为x(f)＝[x1(f)，x2(f)，…，xm(f)]^t，第i个激发源到达第n个台站的延迟时间为τni，这个延迟时间是与反方位角和慢度的函数。可以用(2)公式计算。第n个台站记录到的位移谱dn(f)，可以认为是所有潜在的激发源谱经过相位校正后在该台站的线形叠加，这里不考虑衰减，如公式(3)所示：

对于n个子台，我们可以建立带有e(f)高斯噪声谱的方程组如(4)式

d(f)＝g(f)x(f)+e(f)(4)

其中，g(fh是测量矩阵，其中第n行，第m列相应元素可以表示为τnm是第m个网格点到第n个台站的延时，可以用公式(2)计算。

通常情况下，慢度-反方位角空间上的网格点数(或者说潜在的源个数)m是大于我们的测量数n，那么方程组(4)这个线形系统是欠定的，有无穷多个解。我们解向量x在慢度-反方位角空间是稀疏的，利用这个约束可以大大较少解的不唯一性。解稀疏可以用l1范数测量：

xsol(f)＝argmin||x(f)||1s.t.||d(f)-g(f)x(f)||2＜∈(5)

其中，这里我们采用正交匹配追踪(orthogonalmatchingpursuit，omp)来解决公式(5)(troppetal.,2007)，该方法快速，并且容易实现。

与现有技术相比，本发明的有益效果是：本发明相比传统的台阵fk算法可以提供高分辨率的地震事件p波反方位角和慢度参数。利用源在空间分布的稀疏特性的约束，能够快速反演出地震事件的p波反方位角和慢度参数。

附图说明

为了便于本领域技术人员理解，下面结合附图对本发明作进一步的说明。

图1为地震事件远离2维台阵时的示意图，可以假设地震波为平面波入射到台阵。

图2为上海台阵的布设示意图，其中三角形为台站，s09为参考台站，即笛卡尔坐标系的原点。

图3上海台阵记录到的汶川地震波形示意图，灰色区域为截取的时间窗口，作为本发明的数据输入。

图4为图3中截取的波形对应的傅立叶分析示意图，主频在0.4hz，为本发明所用的频率。

图5为fk和本发明最终结果在慢度-反方位角空间的对比示意图，左图为传统fk方法得到的结果，右图为本发明得到的结果。

图6为fk和本发明最终结果在慢度方向和反方位角方向的结果对比示意图，虚线为根据汶川地震震中计算的实际反方位角。

具体实施方式

如图1所示，一种高分辨率确定地震事件p波反方位角和慢度的方法，该方法具体包括以下步骤：

步骤一：设台阵由n个台组成，取离台阵中心最近的台作为参考台，设为原点；

步骤二：读取n个台的垂直分量地震位移波形记录，并对记录进行去仪器响应，去均值处理，重采样到相同采样率；

步骤三：对步骤二中的所有台的波形记录按发阵时刻对齐，选取时间窗口；

步骤四：对选取的地震记录加窗后进行傅立叶变换，得到记录谱即观测向量d(f)；

步骤五：将慢度和反方位角参数化，选取慢度以及反方位角范围进行均匀网格划分；

步骤六：根据子台位置，以及慢度和反方位角参数化的信息，构建一个测量矩阵g(f)；

步骤七：通过测量矩阵g(f)以及步骤四得到的观测向量d(f)，应用正交匹配追踪方法求出解向量xsol；

步骤八：通过查寻解向量xsol中非零元素对应的索引值，反推到慢度-反方位角空间中，得到最终的解。

进一步地，步骤二中所述时间窗口将每道的p波包含在内。

进一步地，该方法的工作原理如下：

地震远离台阵的情况下，假设地震波为平面波的情况下，以反方位角φ和慢度s(s/km)到达台阵下方。假设2维地震台阵由n子台构成，这里不考虑子台间的垂直方向的高程差。以台阵的中心最近的台站为参考台为坐标原点建立笛卡尔坐标系，第i个子台的位置向量为ri＝[xiyi]。对于所有子台接收到的数据定义为d＝[d1(t)，…，di(t)，…，dn(t)]矩阵，矩阵大小为ns×n，其中ns是在时间维度的点数，di(t)表示为第i个子台接收到的时间序列，为一向量，长度为ns。

在频率域中，对于给定的频率f0，聚束能量e可以用公式(1)计算：

其中j表示虚部单位,di(f0)为数据di(t)在f0的傅立叶变换，相移项中的τi表示在第i个台站的与原点的到时差，可以用公式(2)计算(schweitzeretal.2002),

τi＝-(xisinφ+yicosφ)s(2)

在实际情况中，给定慢度和反方位角的搜寻范围，我们将慢度和反方位角网格划分为m＝sslw×bbaz,其中sslw和bbaz分别表示慢度网格点数和反方位角网格点数，对于fk方法，我们可以搜寻网格点来使得聚束能量最大，该网格点对应的慢度和反方位角为我们的解。但这种方法给出的结果具有低分辨率的特点，并且不能够轻易的区分多个源混合的现象。

假设每隔网格点为一个激发源，具备相应的慢度和反方位角。那么这个问题就转化为在慢度-反方位角空间定位的问题(yaoetal.,2011)，每隔网格点的激发源的谱可以简单表示为x(f)＝[x1(f)，x2(f)，…，xm(f)]^t，第i个激发源到达第n个台站的延迟时间为τni，这个延迟时间是与反方位角和慢度的函数。可以用(2)公式计算。第n个台站记录到的位移谱dn(f)，可以认为是所有潜在的激发源谱经过相位校正后在该台站的线形叠加，这里不考虑衰减，如公式(3)所示：

对于n个子台，我们可以建立带有e(f)高斯噪声谱的方程组如(4)式

d(f)＝g(f)x(f)+e(f)(4)

其中，g(f)是测量矩阵，其中第n行，第m列相应元素可以表示为τnm是第m个网格点到第n个台站的延时，可以用公式(2)计算。

通常情况下，慢度-反方位角空间上的网格点数(或者说潜在的源个数)m是大于我们的测量数n，那么方程组(4)这个线形系统是欠定的，有无穷多个解。我们解向量x在慢度-反方位角空间是稀疏的，利用这个约束可以大大较少解的不唯一性。解稀疏可以用l1范数测量：

xsol(f)＝argmin||x(f)||1s.t.||d(f)-g(f)x(f)||2＜∈(5)

其中这里我们采用正交匹配追踪(orthogonalmatchingpursuit，omp)来解决公式(5)(troppetal.,2007)，该方法快速，并且容易实现。

如图2-6所示，以上海台阵接收到的汶川地震为例，设台站s09为坐标原点，该方法具体包括以下步骤：

步骤一：读取16个台站垂直分量的地震位移波形记录，并对记录进行去仪器响应，去均值处理。重采样到相同采样率。

步骤二：对步骤一的所有台的波形记录按发震时刻对齐，选取时间窗口，该时间窗口不宜过大，也不宜过小，恰好将每道的p波包含在内，选取7s(图3中的灰色阴影部分)。

步骤三：对步骤二中的选取的记录加窗进行傅立叶变换，得到记录谱即观测向量d(f)，在这里根据图4中的频谱，频率f选择0.4hz。

步骤四：将慢度和反方位角参数化，慢度范围为0-0.2，间隔为0.01；反方位角范围为0-360度，间隔为2度。

步骤五：根据子台坐标，以及慢度和反方位角参数化的信息，构建一个测量矩阵g(f)，该测量矩阵的第n行第m列的元素为下式：

其中，τnm是第n个慢度-反方位角空间上的网格点到第m个台站的延时，可以根据公式：τi＝-(xisinφ+yicosφ)s求得。

步骤六：通过测量矩阵g(f)以及观测向量d(f)，应用正交匹配追踪方法解下列l1范数问题，求出解向量xsol

xsol(f)＝argmin||x(f)||1s.t.||d(f)-g(f)x(f)||2＜∈

步骤七：通过查寻解向量xsol中非零元素对应的索引值，反推到慢度-反方位角空间中，得到最终的解在慢度和反方位角值，如图5中的右图以及图6所示。

以上内容仅仅是对本发明结构所作的举例和说明，所属本技术领域的技术人员对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离发明的结构或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。