张 鑫 谭 凯
（中国地震局地震研究所，湖北 武汉 430071）
【摘 要】提出一种可以同时采用近/远场地震波反演地震破裂过程的方法，以便可以更好地利用近场高频GNSS数据及远场地震波形数据。在一个给定的区间内以一定步长搜索破裂传播速度，对于每一个破裂速度，采用遗传算法得到归一化的子断层震源时间函数，通过线性反演快速得到各子断层的滑动分布。采用2013年芦山6.6级地震的源-场模型，进行了数值模拟实验，验证了本方法的可行性与稳定性。
教育期刊网 http://www.jyqkw.com
关键词 破裂过程反演；遗传算法；高频GNSS
※基金项目：中国地震局地震研究所所长基金（IS201116013）。
０ 引言
高频GNSS数据已经被越来越广泛地应用于地震学研究。同时，也有越来越多的连续运行GNSS跟踪站被布置到地震带附近，使得我们可以得到许多同震的近场高频GNSS数据。如果仅仅使用这些数据来获得静态同震形变，就会损失数据中的时域信息。本文提出的方法可以将同震高频GNSS时间序列携带了地震震源破裂过程的信息提取出来，以服务于震源物理研究。
本文所采用的方法基于遗传算法，将计算过程并行化，并对计算过程进行了重要的分解，提前计算出一些关键的矩阵，能够快速计算每一个“个体”的适应度，从而大大减少了计算所需时间。
1 正演
将断层面划分为m个子断层。假定纯剪切破裂起始于震中，并以一个（待求的）恒定速度传播，同时假定每个子断层都有统一的震源时间函数（除了破裂开始时间和幅值），每个子断层的开始破裂时间由且仅由破裂传播速度决定，而各子断层震源时间函数的幅值由各子断层的滑动量决定。这个假定也意味着各子断层在破裂过程中不会改变滑动方向。每个子断层的滑动矢量可分解为两个标量：Xm和Ym，其中Xm为滑动矢量在断层走向上的投影，而Ym为滑动矢量在断层倾向上的投影。
采用n个上升时间为Δt的斜坡函数来表示各子断层正则化的震源时间函数，其中第一个斜坡函数的幅值a1为1，接下来各斜坡函数的幅值为a2，a3，…，an。
台站xj处的理论位移波形由下式表示：
Ｆk（ti，xj）=∑mmanXmfmk（ti-（n-1）Δt-Tm，xj）+∑mmanXmgmk（ti-（n-1）Δt-Tm，xj）（１）
其中fmk和gmk为第m个子断层上滑动角为0°和90°的单位滑动对应格林函数的第k分量，Tm为第m个子断层的起始破裂时间，它由破裂速度决定。
2 反演
其中F0为观测波形向量，D为观测波形的协方差矩阵，?茁为经验阻尼系数，Gn为某破裂速度下的格林函数矩阵。
这是一个非线性问题。它可以被线性化，通过迭代求解，但这种方法对迭代初始值的选取非常敏感，很容易终止在局部极值处。
固定（a2，…，an），将此问题变成一个典型的二次型问题，其最小值由下式得到：
上式为一个关于（a2，…，an）的函数，通过遗传算法来得到它的最小值以及对应的参数（a2，…，an）。
在遗传算法中，每个候选的矢量（a2，…，an）就是一个“个体”，一定数量的个体组成一个“种群”，而整个遗传、变异、选择、交叉的过程是以种群而非个体为单位。对每个个体，根据（3）式来计算它的适应度，并以此来确定每个个体被选择到下一代的概率，使得“更优”的个体有更大的概率被选择保存下来。
由于矩阵Gn过于庞大，若每次将它代入（3）式计算适应度，会使得计算十分缓慢。因此，在遗传算法开始之前，先计算出
和，前者是一个小得多的方阵，后者则是一个向量。它们的大小不会随着观测值的增加而变化，而只与子断层数量有关。这一环节至关重要，它大大减小了遗传算法中的计算量，使得本方法有能力使用大量观测数据进行反演。
通过遗传算法找到使得（3）式最小的（a2，a3，…）后，便可据下式求得滑动量的最优估计：
3 数值实验
制作一个检测板，假定一个单位斜坡函数为各子断层统一的震源时间函数，并设破裂从震中开始以3.0km/s的速度传播，采用2013年芦山6.6级地震的源-场模型，选取了26个GPS台站和35个地震台站作为位移/速度波形的接收场点。在合成的波形中，加入RMS为波形平均振幅5%的高斯误差。以此为伪观测波形，对震源破裂过程进行了反演。
数值实验对破裂速度从2.6km/s到3.5km/s以步长0.1km/s进行了搜索。在上图中可以看到，波形拟合不符度在“真实”破裂速度3.0km/s时取得最小值。在搜索3.0km/s的破裂速度时，对应于各个搜索破裂速度的滑动分布几乎完全恢复了检测板给定的滑动分布。在遗传算法的迭代过程中，虽然最大演化代数被设为4000代，但适应度收敛得非常快，在前100代就完全收敛到了最小值上。
数值实验表明，本文所采用的方法对破裂速度的探测是灵敏的，对误差具有很好的抵抗能力，是一种稳定的反演方法。
4 结论与讨论
本研究是对GNSS数据时域信息进行利用的一次理论尝试。数值实验显示，我们的方法对破裂速度敏感，而且对误差也具有稳定性。
对于该方法来说，台站的分布十分重要。在一些数值实验中某些场点分布下，即使是采用未加入噪声的理论合成波形进行反演，也无法得到“真实的”滑动分布和子断层震源时间函数。因此，对每一个震例来说，对于源-场位置关系进行数值实验都是必要的。
该方法固定了断层面的走向和倾角，这限定了子断层的位置必需位于这个平面上，并且对子断层滑动矢量作了限制，即滑动必需平行于我们设定的断层面，这种假定虽然在物理上是合理的，实际上却增加了先验信息对结果的影响。本文的方法中，对子断层位置的限制和对滑动方向的限制在数学上是并不直接相关的，子断层的位置要被用于格林函数和地震波走时差的计算，所以需要假定已知，但是对于滑动方向的限制并不是必需的。应此可以将各子断层的地震矩张量设为参数，以取代子断层滑动的走滑、倾滑分量，并计算地震矩张量各分量对应的格林函数来组成观测方程。
教育期刊网 http://www.jyqkw.com
参考文献
［１］Hirahara K， Nakano T， Hoso Y， et al. An experiment for GPS strain seismometer[C]//Japanese symposium on GPS. 1994：67-75.
［２］Ji C， Larson K M， Tan Y， et al. Slip history of the 2003 San Simeon earthquake constrained by combining 1‐Hz GPS， strong motion， and teleseismic data[J]. Geophysical research letters， 2004，31（17）.
［３］Larson K M. GPS seismology[J]. Journal of Geodesy， 2009， 83（3-4）：227-233.
［４］Hartzell S H， Heaton T H. Inversion of strong ground motion and teleseismic waveform data for the fault rupture history of the 1979 Imperial Valley， California， earthquake[J]. Bulletin of the Seismological Society of America， 1983，73（6A）：1553-1583.
［责任编辑：邓丽丽］