CN109143360B - 一种高分辨率确定地震事件p波反方位角和慢度的方法 - Google Patents

一种高分辨率确定地震事件p波反方位角和慢度的方法 Download PDF

Info

Publication number
CN109143360B
CN109143360B CN201811101677.9A CN201811101677A CN109143360B CN 109143360 B CN109143360 B CN 109143360B CN 201811101677 A CN201811101677 A CN 201811101677A CN 109143360 B CN109143360 B CN 109143360B
Authority
CN
China
Prior art keywords
slowness
station
azimuth angle
stations
wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811101677.9A
Other languages
English (en)
Other versions
CN109143360A (zh
Inventor
张海江
胡景
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
ANHUI WANTAI GEOPHYSICAL TECHNOLOGY Co Ltd
Original Assignee
ANHUI WANTAI GEOPHYSICAL TECHNOLOGY Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by ANHUI WANTAI GEOPHYSICAL TECHNOLOGY Co Ltd filed Critical ANHUI WANTAI GEOPHYSICAL TECHNOLOGY Co Ltd
Priority to CN201811101677.9A priority Critical patent/CN109143360B/zh
Publication of CN109143360A publication Critical patent/CN109143360A/zh
Application granted granted Critical
Publication of CN109143360B publication Critical patent/CN109143360B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种高分辨率确定地震事件P波反方位角和慢度的方法，该方法的工作步骤：设台阵由N个台组成，取离台阵中心最近的台作为参考台，设为原点，读取N个台的垂直分量地震位移波形记录，并对记录进行去仪器响应，去均值处理，重采样到相同采样率，对步骤二中的所有台的波形记录按发阵时刻对齐，选取时间窗口，对选取的地震记录加窗后进行傅立叶变换，得到记录谱即观测向量。本发明的有益效果是：本发明相比传统的台阵FK算法可以提供高分辨率的地震事件P波反方位角和慢度参数。利用源在空间分布的稀疏特性的约束，能够快速反演出地震事件的P波反方位角和慢度参数。

Description

一种高分辨率确定地震事件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＝[xi yi]。对于所有子台接收到的数据定义为d＝[d1(t)，…，di(t)，…，dN(t)]矩阵，矩阵大小为ns×N，其中ns是在时间维度的点数，di(t) 表示为第i个子台接收到的时间序列，为一向量，长度为ns。
在频率域中，对于给定的频率f0，聚束能量E可以用公式（1）计算：
Figure DEST_PATH_BDA0001806832070000031
其中j表示虚部单位,Di(f0)为数据di(t)在f0的傅立叶变换，相移项中的τi表示在第i个台站的与原点的到时差，可以用公式（2）计算 (Schweitzer et al. 2002),
τi＝-(xisinΦ+yicosΦ)s (2)
在实际情况中，给定慢度和反方位角的搜寻范围，我们将慢度和反方位角网格划分为M＝Sslw×Bbaz,其中Sslw和Bbaz分别表示慢度网格点数和反方位角网格点数，对于FK方法，我们可以搜寻网格点来使得聚束能量最大，该网格点对应的慢度和反方位角为我们的解。但这种方法给出的结果具有低分辨率的特点，并且不能够轻易的区分多个源混合的现象。
假设每隔网格点为一个激发源，具备相应的慢度和反方位角。那么这个问题就转化为在慢度-反方位角空间定位的问题（Yao et al., 2011），每隔网格点的激发源的谱可以简单表示为X(f)＝[x1(f)，x2(f)，…，xM(f)]T，第i个激发源到达第n个台站的延迟时间为τni，这个延迟时间是与反方位角和慢度的函数。可以用（2）公式计算。第n个台站记录到的位移谱dn(f)，可以认为是所有潜在的激发源谱经过相位校正后在该台站的线形叠加，这里不考虑衰减，如公式（3）所示：
Figure DEST_PATH_BDA0001806832070000041
对于N个子台，我们可以建立带有E(f)高斯噪声谱的方程组如（4）式
Figure 1952DEST_PATH_IMAGE029
其中，G(f) 是测量矩阵，其中第n行，第m列相应元素可以表示为
Figure DEST_PATH_BDA0001806832070000042
Figure DEST_PATH_BDA0001806832070000043
τnm是第m个网格点到第n个台站的延时，可以用公式（2）计算。
通常情况下，慢度-反方位角空间上的网格点数（或者说潜在的源个数）M是大于我们的测量数N，那么方程组（4）这个线形系统是欠定的，有无穷多个解。我们解向量X在慢度-反方位角空间是稀疏的，利用这个约束可以大大减少解的不唯一性。解稀疏可以用L1范数测量：
Xsol(f)＝argmin||X(f)||1s.t.||d(f)-G(f)X(f)||2＜∈ （5）
其中，
Figure DEST_PATH_BDA0001806832070000044
这里我们采用正交匹配追踪（OrthogonalMatching Pursuit，OMP）来解决公式（5） (Tropp et al., 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＝[xi yi]。对于所有子台接收到的数据定义为d＝[d1(t)，…，di(t)，…，dN(t)]矩阵，矩阵大小为ns×N，其中ns是在时间维度的点数，di(t) 表示为第i个子台接收到的时间序列，为一向量，长度为ns。
在频率域中，对于给定的频率f0，聚束能量E可以用公式（1）计算：
Figure DEST_PATH_BDA0001806832070000061
其中j表示虚部单位,Di(f0)为数据di(t)在f0的傅立叶变换，相移项
Figure DEST_PATH_BDA0001806832070000062
中的τi表示在第i个台站的与原点的到时差，可以用公式（2）计算 (Schweitzer et al. 2002),
τi＝-(xisinΦ+yicosΦ)s (2)
在实际情况中，给定慢度和反方位角的搜寻范围，我们将慢度和反方位角网格划分为M＝Sslw×Bbaz,其中Sslw和Bbaz分别表示慢度网格点数和反方位角网格点数，对于FK方法，我们可以搜寻网格点来使得聚束能量最大，该网格点对应的慢度和反方位角为我们的解。但这种方法给出的结果具有低分辨率的特点，并且不能够轻易的区分多个源混合的现象。
假设每隔网格点为一个激发源，具备相应的慢度和反方位角。那么这个问题就转化为在慢度-反方位角空间定位的问题（Yao et al., 2011），每隔网格点的激发源的谱可以简单表示为X(f)＝[x1(f)，x2(f)，…，xM(f)]T，第i个激发源到达第n个台站的延迟时间为τni，这个延迟时间是与反方位角和慢度的函数。可以用（2）公式计算。第n个台站记录到的位移谱dn(f)，可以认为是所有潜在的激发源谱经过相位校正后在该台站的线形叠加，这里不考虑衰减，如公式（3）所示：
对于N个子台，我们可以建立带有E(f)高斯噪声谱的方程组如（4）式
其中，G(f) 是测量矩阵，其中第n行，第m列相应元素可以表示为
Figure DEST_PATH_BDA0001806832070000072
τnm是第m个网格点到第n个台站的延时，可以用公式（2）计算。
通常情况下，慢度-反方位角空间上的网格点数（或者说潜在的源个数）M是大于我们的测量数N，那么方程组（4）这个线形系统是欠定的，有无穷多个解。我们解向量X在慢度-反方位角空间是稀疏的，利用这个约束可以大大减少解的不唯一性。解稀疏可以用L1范数测量：
Xsol(f)＝argmin||X(f)||1s.t.||d(f)-G(f)X(f)||2＜∈ （5）
其中这里我们采用正交匹配追踪（OrthogonalMatching Pursuit，OMP）来解决公式（5） (Tropp et al., 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列的元素为下式：
Figure DEST_PATH_BDA0001806832070000081
其中，τ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所示。
以上内容仅仅是对本发明结构所作的举例和说明，所属本技术领域的技术人员对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离发明的结构或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。

Claims (2)

1.一种高分辨率确定地震事件P波反方位角和慢度的方法，其特征在于，该方法具体包括以下步骤：
步骤一：设台阵由N个台组成，取离台阵中心最近的台作为参考台，设为原点；
步骤二：读取N个台的垂直分量地震位移波形记录，并对记录进行去仪器响应，去均值处理，重采样到相同采样率；
步骤三：对步骤二中的所有台的波形记录按发阵时刻对齐，选取时间窗口；
步骤四：对选取的地震记录加窗后进行傅立叶变换，得到记录谱即观测向量d(f)；
步骤五：将慢度和反方位角参数化，选取慢度以及反方位角范围进行均匀网格划分；
步骤六：根据台位置，以及慢度和反方位角参数化的信息，构建一个测量矩阵G(f)；
步骤七：通过测量矩阵G(f)以及步骤四得到的观测向量d(f)，应用正交匹配追踪方法求出解向量Xsol；
步骤八：通过查寻解向量Xsol中非零元素对应的索引值，反推到慢度-反方位角空间中，得到最终的解。
2.根据权利要求1所述的一种高分辨率确定地震事件P波反方位角和慢度的方法，其特征在于，步骤三中所述时间窗口将每道的P波包含在内。
CN201811101677.9A 2018-09-20 2018-09-20 一种高分辨率确定地震事件p波反方位角和慢度的方法 Active CN109143360B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811101677.9A CN109143360B (zh) 2018-09-20 2018-09-20 一种高分辨率确定地震事件p波反方位角和慢度的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811101677.9A CN109143360B (zh) 2018-09-20 2018-09-20 一种高分辨率确定地震事件p波反方位角和慢度的方法

Publications (2)

Publication Number Publication Date
CN109143360A CN109143360A (zh) 2019-01-04
CN109143360B true CN109143360B (zh) 2020-01-10

Family

ID=64815441

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811101677.9A Active CN109143360B (zh) 2018-09-20 2018-09-20 一种高分辨率确定地震事件p波反方位角和慢度的方法

Country Status (1)

Country Link
CN (1) CN109143360B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109856671A (zh) * 2019-03-06 2019-06-07 合肥国为电子有限公司 一种基于无线通讯的地震探测方法及系统
CN112462415B (zh) * 2020-11-02 2023-07-21 中国电子科技集团公司第三研究所 一种对多振动源进行定位的方法及装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2444954B (en) * 2006-12-20 2009-05-20 Westerngeco Seismic Holdings Method of monitoring microseismic events
CN107180512B (zh) * 2017-06-22 2018-12-11 禁核试北京国家数据中心 一种特定地区地震事件的报警方法

Also Published As

Publication number Publication date
CN109143360A (zh) 2019-01-04

Similar Documents

Publication Publication Date Title
Gubbins Time series analysis and inverse theory for geophysicists
Zha et al. Determining the orientations of ocean bottom seismometers using ambient noise correlation
US20040243312A1 (en) Method of correcting for time shifts in seismic data resulting from azimuthal variation
CN109143360B (zh) 一种高分辨率确定地震事件p波反方位角和慢度的方法
Van Dalen et al. Retrieving surface waves from ambient seismic noise using seismic interferometry by multidimensional deconvolution
Yang et al. Determination of the local magnitudes of small earthquakes using a dense seismic array in the Changning–Zhaotong Shale Gas Field, Southern Sichuan Basin
Lokmer et al. Properties of the near-field term and its effect on polarisation analysis and source locations of long-period (LP) and very-long-period (VLP) seismic events at volcanoes
CN112444773A (zh) 基于空域融合的压缩感知二维doa估计方法
Vinnik Receiver function seismology
Poppeliers et al. Three‐dimensional seismic‐wave gradiometry for scalar waves
CA2331712C (en) Method of and system for processing multicomponent seismic data
AU2002317974B2 (en) A method of processing geophysical data
Haldorsen Spatial aliasing and 3C seismic sensors
Lehujeur et al. Eikonal Tomography Using Coherent Surface Waves Extracted From Ambient Noise by Iterative Matched Filtering—Application to the Large‐N Maupasacq Array
Muyzert et al. A five component land seismic sensor for measuring lateral gradients of the wavefield
US6826485B1 (en) Determination of the fast and slow shear wave polarisation directions
NIU et al. The Polynomial Radon Transform
Poppeliers et al. Three‐dimensional wave gradiometry for polarized seismic waves
CN111830560B (zh) 一种基于降秩算法的地震数据重建方法
CN110824555B (zh) 地震能量均衡方法及装置、计算机可读存储介质
Kabir et al. Migration velocity analysis using the common focus point technology
Urban et al. Rethinking the polar cap: Eccentric dipole structuring of ULF power at the highest corrected geomagnetic latitudes
WO2001036999A2 (en) Determination of the fast and slow shear wave polarisation directions
Dellinger et al. Alford rotation after tensor migration
CN111538084B (zh) Ovt域数据转换成方位角度域成像道集的方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant