0 引言

天然地震监测，相比于人工地震勘探，更容易受复杂的地质条件以及采集环境的影响，通常野外地震台站布设稀疏并且呈不规则分布，这使得对地震的监测能力以及对地下介质成像的分辨率受到地震台站数量和记录的地震数据质量的影响.近年来，随着监测仪器性能的提升和轻便化，天然地震(包括背景噪声)的观测向着密集台阵观测的趋势发展.例如南加州长滩的密集台阵观测采用了5200台，台间距大约为100 m(Lin et al., 2013).美国San Helen火山的密集台阵观测用到了904台，台间距大约为10 m (Hansen and Schmandt, 2015).美国南加州San Jacinto断层的密集台阵观测布设了1108台，台间距为10m(Ben-Zion et al., 2015).随着密集台阵观测的增多，原来由于台站稀疏而无法采用的地震处理方法，例如地震偏移成像，就可以用于天然地震数据处理.但是天然地震的密集台阵观测依然受可布设台站的数量和监测区域地形的影响，而呈现非规则分布并且对于某些地震数据处理方法来说依旧过于稀疏.为此，本文针对天然地震数据的重建或者规整化展开研究，并且以实际布设在南加州San Jacinto断层的密集台阵为例.

人工地震勘探数据的重建或者规整化，在近十几年来得到了快速发展.这是因为随着勘探工作的深入推进，连片偏移成像处理(Gierse et al., 1949)、宽方位地震数据OVT(Offset Vector Tilt，偏移矢量片)处理(Ling et al., 2016)、逆时偏移(Reverse Time Migration，RTM)(Zhou et al., 2011; Li et al., 2016)、全波形反演(Full Waveform Inversion，FWI)(Anderson et al., 2012; Chi et al., 2015)、多分量地震(赵波等, 2012)等面向地下构造的成像和参数估计方法技术，以及AVO(Amplitude Variation with Offset)反演(Downton et al., 2010)、横波分裂分析(Stanton and Sacchi, 2013)等面向目标储层的反演和属性分析方法技术都对地震勘探数据的规则性和完整性提出了更高的要求.在理论方面，基于一定模型的地震数据重建是经典的反问题.在实际应用方面，地震数据作为后续处理和反演的先天之本，对其进行实现高信噪比、高分辨率和高保真度的重建具有十分重要的现实意义(Fu et al., 2018).

经典的地震数据重建方法中，基于稀疏约束的信号重构算法利用各种变换方法寻求信号的稀疏表达，从而重建缺失的地震信号.基于此类算法的信号重构最关键问题是寻找信号的最稀疏表达方法，已有众多的文章研究适合各种不同特征信号的稀疏变换方法，如小波变换、曲波变换、自适应字典等(曹静杰等, 2012; Mansour et al., 2013; Chen et al., 2013).运用降秩的方法对地震数据重建的思想与基于稀疏性约束的信号重构思想较为一致，是在秩空间中寻求地震数据的稀疏表达.

近几年，基于低秩理论的矩阵补全方法受到广泛关注.杨等人发现了地震数据重建和低秩补全问题之间的联系(Yang et al., 2014)，已经成功地将矩阵补全问题应用于人工地震数据重建和噪声衰减.降秩方法通常假设地震数据在经过一些预转换后具有低秩结构，比如频率切片上构造Hankel矩阵(Chen and Sacchi, 2014; Naghizadeh and Sacchi, 2013; Oropeza and Sacchi, 2011)，然而地震道的缺失和随机噪声的存在会增加预变换数据矩阵的秩.因此，地震数据重建的问题归结为矩阵降秩问题，可以使用降秩方法来重建缺失的地震道数据.由于秩函数最小化本身是一个NP难问题，比较难求解，因此凸松弛化为核范数(奇异值的l₁-范数)最小化求解.Trickett等人将Cadzow方法用于地震数据插值，该方法首先将数据变换到频率(f-x)域，再对每个频率切片生成Hankel矩阵，最后对Hankel矩阵进行降秩处理，该方法保证Hankel矩阵对于包含线性同相轴的地震数据是低秩的(Trickett et al., 2010).Oropeza等提出多道奇异谱分析(Multichannel Singular Spectrum Analysis，MSSA)方法，利用奇异值分解方法进行矩阵降秩处理(Oropeza and Sacchi, 2011). Naghizadeh等提出了去假频的Cadzow/MSSA方法，主要将低频成分的Hankel矩阵分解结果用于恢复高频Hankel矩阵分解中，从而达到去假频、重建地震数据的目的(Naghizadeh and Sacchi, 2009).然而，这些降秩方法都需要基于奇异值分解(SVD)方法，比较耗时.Gao等利用块Toeplitz矩阵替代Hankel矩阵，提出了一种快速的Cadzow/MSSA方法(Gao et al., 2013).正交秩-1矩阵追踪(OR1MP)算法是一种基于矩阵基追踪的快速降秩方法，只有左右顶部的奇异值向量需要进行计算，从而避免了SVD分解所需的计算复杂度(Wang et al., 2014)，因此采用OR1MP算法解决核范数最小化问题.

上述的地震数据重建方法主要用于人工地震勘探数据重建.与人工地震数据重建相比，天然地震数据的重建研究较少，处于刚起步阶段.Zhang等人将图像处理中常用的凸集投影算法(Projection Onto Convex Sets，POCS)运用到天然地震波形的修复，并对芦山地震发生在近场的剪切波形进行了重建(Zhang et al., 2016).相比较来说，天然地震与人工地震的震源机制是不一样的，天然地震事件多为剪切破裂源，而人工地震多为爆炸源，能量辐射方式不一样.通常天然地震发生地点距离观测台站较远，考虑到震源辐射花样，相对于震源不同方位的检波器所接收到的波形振幅与相位有较大差别.在实际应用中需分方位进行重建，即将数据体按照震源辐射区域进行分割，对分割后的数据体重建后再整合成完整的数据体.在相同的辐射区内，可以认为是平面波并且起跳的相位相同.因此对于天然地震波形，尽管不是爆炸源，也可以得到较为一致的排列，可以使用人工地震重建的方法对天然地震数据进行重建.在本文中，由于检波器范围仅为600m，而天然地震震源在检波器范围之外几十公里.故认为检波器在天然地震震源辐射的相同区域，进而直接将整个数据体进行数据重建.

本文用基于低秩理论的正交秩-1矩阵追踪(OR1MP)图像处理技术重建天然地震数据.这种方法在图像补全和空间欠采样的人工地震数据重建中已经得到成功应用(Naghizadeh and Sacchi, 2009; Jia et al., 2016).本文首先介绍基于低秩理论的天然地震数据重建的基本原理，然后基于美国南加州San Jacinto断层的密集台阵数据测试2D和3D的天然地震数据重建.理论分析和数值实验表明，OR1MP方法可以很好地对天然地震数据进行重建.

1 基于降秩理论的天然地震数据重建基本原理

降秩方法通常假设地震数据在经过一些预变换后具有低秩结构，比如对频率切片构造Hankel矩阵，由于地震道的缺失和随机噪声的存在会增加预变换数据矩阵的秩，因此，地震数据重建的问题可归结为矩阵降秩问题.3D的天然地震数据D由时间-空间域中的两个空间维度和一个时间维度表示，假设D的维数为n_t×n_y×n_x，其中n_t表示时间维度的大小，n_y和n_x表示两个空间维度的大小.固定D的一个空间维度n_y即可得到2D的天然地震数据d，维数为n_t×n_x.由于Hankel变换是对频率切片而言的，首先对观测到的时间域的天然地震信号D和d沿时间轴做一维傅里叶变换.对于2D的天然地震数据d，向量m_ω ∈ R ^n_x表示在f-x域对于给定的频率ω的频率切片，对m_ω构造一个Hankel矩阵.对于3D的天然地震数据D，矩阵 M_ω ∈ R ^n_y×n_x表示在f-x-y域对于给定的频率ω的频率切片，需要对M_ω构造一个块Hankel矩阵.因此，我们首先讨论Hankel矩阵和块Hankel矩阵的构造方法.

1.1 2D天然地震数据频率切片Hankel矩阵的构造

定义M_ω=[m₁, m₂, …, m_nx]是给定频率ω的2D的天然地震数据的频率切片，其中m_i(i=1, 2, …, n_x)是经过傅里叶变换后给定频率ω的地震数据，由M_ω生成的Hankel矩阵定义如下：

(1)

其中，定义为向下取整，即取不大于x的最大整数.Hankel矩阵近似为方阵，其反对角线上的元素是相等的.

1.2 3D天然地震数据频率切片块Hankel矩阵的构造

3D天然地震数据给定频率切片ω上的数据 M_ω ∈ R ^n_y×n_x，n_x和n_y表示空间x和y方向上的地震道数，M_ω给定如下：

(2)

其中，M_ω的每一列元素都是一个大小为n_y的向量，应用上述频率切片m_ω的Hankel操作，M_ω的每一列元素都可以构造一个Hankel矩阵.定义M_ω的第j列生成的Hankel矩阵 H_j(j=1, 2, …, n_x)，这些Hankel矩阵可以看作一个新的向量P：

(3)

再对 P 作上述Hankel操作，可以得到块Hankel矩阵 Y ，表示如下：

(4)

其中，同样定义为向下取整.上面构造的 Y 矩阵是一个块Hankel矩阵，它的大小为(L_x×L_y)×(K_x×K_y)，每一个元素都是一个小的Hankel矩阵，Hankel矩阵的构造过程如图 1所示.图 1描述了一个5×5的频率切片，对第3行的向量构造一个3×3的Hankel矩阵和把5×5的频率切片转换为9×9的块Hankel矩阵的过程.

2 OR1MP算法

对于由2D天然地震数据频率切片构造的Hankel矩阵S和3D天然地震数据的频率切片构造的块Hankel矩阵 Y来说，在理想条件下S和Y的秩和分析窗口中平面波的数量是对应的.即当数据中含有k个平面波时，对应的Hankel矩阵的秩为k(Hua, 1992).缺失的天然地震数据增加了矩阵S和Y的秩，所以天然地震数据重建问题可以看做矩阵的降秩补全问题.Hankel矩阵S和块Hankel矩阵Y的大小不一样，但都是低秩矩阵补全问题，因此这部分的论述只针对S，对于矩阵Y，只需要把补全问题模型的S换成Y即可.

定义 S ∈ R ^n×m是缺失的地震数据，低秩矩阵补全问题是找到一个低秩矩阵 L ∈ R ^n×m较好地近似矩阵S.可以通过求解以下秩最小化模型来重建随机缺失的数据：

(5)

其中，L ∈ R ^n×m是矩阵S的估计；rank(L)是矩阵L的秩，即矩阵L非零奇异值的个数；Ω{1, …, n}×{1, …, m}是采样数据的指标集；P_Ω(·)定义为对指标集Ω的采样操作，定义如下：

(6)

(5) 式所示的秩最小化模型是NP难问题，现已知的有限时间算法在理论和实践中都至少有双倍的指数运行时间.因此，针对问题(5)提出了许多近似解(Candès and Recht, 2009; Lai et al., 2009)，Candes和Recht提出了用核范数代替上述模型中的rank(·) (Candès and Recht, 2009)，核范数作为秩函数的松弛函数，定义为，因此可以得到如下凸松弛模型：

(7)

其中，是矩阵L的奇异值，k是矩阵L的秩.Candes和Recht从理论上证明在一般的条件下上述凸松弛模型可以以高概率准确地恢复原始低秩矩阵(Candès and Recht, 2009).并且发现核范数是奇异值的和，秩是非零奇异值的个数，因此用矩阵的核范数替代矩阵秩的这种做法类似于基于压缩感知的稀疏信号恢复中用l₁ -范数替代l₀ -范数的做法.通常噪声在观测数据中是不可避免的，所以上述模型(7)经常被放宽到核范数正则化无约束优化问题，对应的优化问题的目标函数如下式所示：

(8)

其中，定义为矩阵的Frobenius范数，即矩阵每个元素的平方和的平方根，λ是一个可以调节的参数.

用正交秩-1矩阵追踪算法(OR1MP)(Wang et al., 2014)求解上述优化问题，任何一个秩为k的矩阵 L ∈ R ^n×m都可以写成秩1矩阵的线性组合，即

(9)

其中，B_i(i=1, 2, …, k)是秩-1基矩阵，并且满足条件是基矩阵的权重系数向量.

假设在第(l－1)次迭代之后，秩为1的基础矩阵 B₁, …, B_l－1和它们当前的权重系数θ^l－1已经完成了计算，那么在第l次迭代过程中，将寻找一个新的满足的秩-1基矩阵 B_l.这主要与当前观察到的回归残留有关，残差计算如下：

(10)

其中，

(11)

因此，B_l可以通过求解下面的优化问题得到

(12)

其中，〈 L, S 〉表示两个矩阵 L ∈ R ^n×m和S ∈ R ^n×m的内积，设 L =(l₁, …, l_m)，ver(L)=(l₁^T, …, l_m^T)^T是矩阵 L 的向量化，则〈 L , S 〉=〈vec(L), vec(S)〉.因为每个满足的秩-1矩阵 B 都可以写成两个单位向量的乘积，B = uv ^T(u ∈ R ⁿ, v ∈ R ^m，问题(12)可以改写为以下优化问题：

(13)

优化问题(13)的最优解〈 u_*, v_*〉是残差 E_l的一对左上和右上奇异值向量，因此新的秩-1基矩阵更新公式如下：

(14)

在找到新的秩-1基矩阵 B_l后，通过求解以下最小二乘回归问题来更新所有当前可用基矩阵{ B₁, …, B_l}的权重系数θ^l：

(15)

θ^l的最优解为

(16)

其中，向量s和 b是由矩阵( S )_Ω和(B)_Ω转换得到的，设矩阵( S )_Ω=(s₁, s₂…, s_m)∈ R ^n×m，其中向量 s_i ∈ R ^n×1 (i=1, 2, …m)，则s=vec((S)_Ω)=(s₁^T, …, s_m^T)^T. 是由所有重建的基向量形成的矩阵.

反复运行上述两个更新步骤，直到满足所设置的迭代终止条件，就可以得到优化问题(8)中的L矩阵.对我们的优化问题来说，有两种常用的迭代终止条件，第一种零迭代次数l=k(地震数据的秩)，或者使得残差E_k小于一个预先设定的参数.算法如下：

算法：基于OR1MP的低秩矩阵补全算法
输入：S，秩k.
初始化：L0=0，θ0=0，l=1.
for l=1:k
步骤一：找到公式(10)残差的一对左上和右上的奇异向量(ul, vl)，令Bl=ulvlT.
步骤二：用公式更新权重θl.
步骤三:
end
输出：

3 基于Hankel矩阵降秩的天然地震数据插值

本文分别对2D天然地震数据和3D天然地震数据进行重建，两种重建方法的区别在于Hankel矩阵的构造和反Hankel矩阵的处理，下面分别介绍2D天然地震数据和3D天然地震数据重建的技术路线.

3.1 2D天然地震数据重建的技术路线

对于空间欠采样的2D天然地震数据，沿着时间轴做傅里叶变换后得到频率域的地震数据，每个频率切片构造的Hankel矩阵如公式(1)所示，Hankel矩阵反对角线上的元素是相等的.首先介绍反Hankel过程，该过程是构造Hankel矩阵的逆过程，对公式(1)构造的Hankel矩阵反对角线取平均得到重建的频率切片数据，这个过程是在最优化矩阵 L 中实现的，反Hankel操作可以得到重建后的频率切片，完成频率切片的重建过程.基于Hankel矩阵的二维地震数据重建的技术路线如下示：

输入：时间域缺失道的2D天然地震数据d∈ R ^n_t×n_x.

第一步：沿时间轴对2D天然地震数据d的每一道做傅里叶变换，得到频率域的2D天然地震数据.

第二步：对每个频率切片m_ω

(1) 按公式(1)对频率切片构造Hankel矩阵S，Hankel矩阵S具有低秩结构，地震的缺失增加了矩阵S的秩;

(2) 用OR1MP算法对矩阵S做降秩重建，得到矩阵L；

(3) 对降秩重建后得到的Hankel矩阵L做反Hankel，得出重建后的频率切片的列元素，完成该频率切片的重建过程，得到频率域数据.

第三步：沿频率轴做逆傅立叶变换.

输出：重建的2D天然地震数据.

3.2 3D天然地震数据重建的技术路线

对于空间欠采样的3D天然地震数据，与之前的2D天然地震处理相比，三维数据体处理包含同样的六步，不同之处在于Hankel矩阵的构造.空间欠采样的3D天然地震数据沿着时间轴做傅里叶变换后得到频率域的地震数据，每个频率切片是一个矩阵，不同于(1)式构造的Hankel矩阵，3D天然地震数据的频率切片需要按公式(4)构造块Hankel矩阵.反块Hankel矩阵处理方法是，首先按块Hankel反对角线取平均，得到n_y个小矩阵，然后再将每个小矩阵按照上述2D天然地震处Hankel矩阵的处理方法反向对角线平均，最终得到重建的频率切片的矩阵元素.基于Hankel矩阵的三维地震数据重建的技术路线如下示：

输入：时间域缺失道的3D天然地震数据 D ∈ R ^{n_t×n_x×n_y}.

第一步：沿时间轴对3D天然地震数据 D 的每一道做傅里叶变换，得到频率域的3D天然地震数据.

第二步：对每个频率切片M_ω

(1) 按公式(4)对频率切片构造块Hankel矩阵，构造的块Hankel矩阵具有低秩结构，地震道的缺失增加了矩阵的秩;

(2) 用OR1MP算法对块Hankel矩阵做降秩重建，得到矩阵重建的块Hankel；

(3) 对降秩重建后得到的Hankel矩阵做逆块Hankel操作，得出重建后的频率切片的矩阵元素，完成该频率切片的重建过程.

第三步：沿频率轴做逆傅立叶变换.

输出：重建的3D天然地震数据.

4 圣哈辛托断层带天然地震数据重建

本节分别对2D和3D的圣哈辛托断层带(San Jacinto Fault Zone，SJFZ)天然地震数据进行重建.为了衡量重建效果，引入信噪比(SNR)定义如下：

(17)

其中，x₀表示原始的数据，x_rec表示重建的数据，SNR值越高，表示重建效果越好.

我们讨论了加州San Jacinto断层密集地震台阵中记录的大于4周的地震记录.其中有1108个垂直(10 Hz)ZLand地震检波器，检波器布置的范围长度约600 m×600 m.该阵列于2014年5月7日至2014年6月13日在Anza南部复杂的三叉区内Sage Brush Flat (SGB)站点部署，核心网格由垂直于克拉克断层轨迹并居中的20排组成，每行有50个传感器，标称距离为10 m，行间标称间距为30 m，剩余的108个检波器被部署为多行的扩展.该阵列以每秒500个采样点连续记录地震和噪声数据，即时间采样率为0.002 s，可用频率高达200 Hz.SGB站点1108个垂直分量ZLand检波器空间上密集的节点阵列如图 2所示，其中，黑色三角形表示检波器，x轴表示北向，y轴表示东向.

4.1 2D天然地震数据重建实验结果

2D天然地震数据实验选取的数据是核心网格第10排的50个检波器采集到的圣哈辛托断层的天然地震数据，如图 2中黑色虚线矩形区域所示，检波器间距为10 m，数据采样间隔为0.002 s，实验数据集都包含512个时间采样和50个空间采样，即实验数据大小是512×50.本实验是对2D天然地震数据进行重建，随机欠采样率从10%增加到50%.利用OR1MP算法对缺失的数据体进行重建，最大迭代次数设置为500，参数秩定义为k=20.在处理实际数据时，不同倾角的数量并不是先验已知的.用小参数k可能会导致不切实际的解决方案，另一方面，增加k会导致复杂的重建.对于实际的地震数据还没有能够制定选择k的客观策略，实验时使用k=18, 19, 20, 21, 22获得了非常相似的结果.选择秩k类似于选择f-x解卷积中预测滤波器的长度(Cadzow, 1988).

图 3和图 4显示了将OR1MP算法应用到San Jacinto断层带采集到的天然地震数据后的结果图，图 3a和图 4a展示的是2D原始天然地震数据体，采样率为20%的圣哈辛托断层带的地震道集如图 3b所示，采样率为40%的圣哈辛托断层带的地震道集如图 4b所示.图 3c和图 4c显示的是用OR1MP算法对欠采样率为20%和40%的San Jacinto断层带天然地震数据体重建后的结果.从图 3c和图 4c可以看出，OR1MP算法对缺失的天然地震数据重建结果与原始的天然地震数据体的一致性均较好.

图 5图 6显示的是将OR1MP算法应用到San Jacinto断层带采集到的天然地震数据后的频谱图.图 5a和图 6a是2D原始天然地震数据体的频谱图，图 5b和图 6b是缺失数据的频谱图，表明地震数据的缺失会产生空间频散现象.而重建后的地震数据频谱图 5c和6c表明重建后频散现象得到了抑制.

为了进一步说明重建效果，一方面我们展示了从原始数据和重建结果得到的第36道地震数据的单道对比图，如图 7所示.从单道结果对比图可以看出重建后的天然地震数据和原始的数据相位基本一致，说明重建结果的有效性.另一方面我们分别记录奇异谱分析(Singular Spectrum Analysis, SSA)算法(Hassani and Zhigljavsky, 2009)和OR1MP算法重建后信噪比(SNR)与采样率之间的关系.如表 1所示，在采样率为20%的情况下，OR1MP算法信噪比由重建之前的6.9727 dB提高到23.3876 dB，而SSA算法的重建后的信噪比是17.7292 dB.采样率为40%的情况下，OR1MP算法重建信噪比(19.9332 dB)比SSA算法重建信噪比(14.7027 dB)高5 dB左右.由图 8信噪比(SNR)的变化曲线图可以看出，相对于经典的SSA(奇异谱分析)算法，OR1MP算法能得到更高精度的重建结果.

4.2 3D天然地震数据重建实验结果

3D天然地震数据实验选取的数据是圣哈辛托断层带核心网格检波器采集到的地震数据，核心网格有20排，行间距为50 m，每一排有50个检波器，行间距为10 m，如图 2中黑色实线矩形区域所示.采样频率为0.002 s，实验数据大小是512×50×20，其中，第一个维度表示时间采样(t轴)，另外两个维度表示空间采样(x-y轴).

本实验对3D圣哈辛托断层带不同欠采样率的地震数据进行重建.秩k=20的选择同上述2D地震数据重建中k的选择，实验结果只显示3D天然地震数据体的2D切片图(t-x切片)，图 9、图 10和图 11分别展示了随机缺失50%的3D圣哈辛托断层带天然地震数据y=5，y=10和y=16的t-x切片图.图 9a、图 10a和图 11a表示原始天然地震数据的切片图，图 9b、图 10b和图 11b是缺失50%的天然地震数据的切片图，对应的重建后天然地震数据的切片图如图 9c、图 10c和图 11c所示.图 9d、图 10d和图 11d表示重建后天然地震数据和原始天然地震数据的残差图.其中横轴表示空间轨迹采样，纵轴表示时间采样.从图 9c、图 11c和图 12c的切片图可以看出ORIMP算法对3D天然地震数据的重建效果很好.

另一方面，多道奇异谱分析(Multichannel Singular Spectrum Analysis, MSSA)算法和OR1MP算法重建信噪比可以进一步说明重建效果差异.表 2记录了不同采样率下信噪比(SNR)数值，采样率为30%的情况下，ORIMP算法信噪比由重建之前的5.2274 dB提高到21.8079 dB，MSSA算法只提高到18.1214 dB.采样率为50%的情况下，ORIMP算法(19.2680 dB)比MSSA算法(14.7548 dB)重建精度高，说明ORIMP算法对3D天然地震数据的重建能取得不错的效果.图 12实验结果进一步表明，OR1MP算法相对于经典的MSSA算法有较高的重建精度.

5 结论

虽然人工地震的数据重建已经有了大量的研究工作，但是天然地震数据的重建目前国际上研究较少.因为天然地震观测台站较为稀疏，而天然地震的数据重建对于地震信号的检测和地下结构的成像具有非常大的意义，所以天然地震数据的重建研究更有紧迫性和重要性.本文把将基于低秩理论的矩阵补全算法(OR1MP)运用到天然地震数据的重建问题，并且基于一个实际的密集台站观测数据对算法进行了测试.实际数据测试表明，对于2D天然地震数据，在数据缺失达到50%的情况下，基于低秩理论的矩阵补全算法(OR1MP)依然能够很好地恢复缺失的信号.在3D天然地震数据观测的情况下，在数据缺失达到70%的情况下，重建算法依然可以高保真地恢复出缺失的信号.也就是说，对于San Jacinto断层的密集台阵天然地震观测，只需要布设300个传感器就几乎可以达到1000个传感器能够达到的效果.可以看出通过天然地震数据重构算法，在某种程度上可以大大减少台站布设的数量，降低台站布设的成本.

致谢  作者感谢南加州大学洛杉矶分校地球科学系的Yehuda Ben-Zion教授提供的数据，以及感谢中国科学技术大学地球与空间科学学院的常凯同学在数据前期准备中的协助和论文修改中给予的建议.