0 引言

据中国地震台网中心测定，北京时间2023年9月9日6时10分，摩洛哥发生M_W6.9地震，震中位于8.55°W，31.00°N，震源深度10km。截至北京时间9月13日10时，本次地震已造成2901人遇难，2530人受伤，震中区域建筑和道路损毁严重。本次地震位于摩洛哥马拉喀什市西南方向，摩洛哥首都拉巴特和马拉喀什等多个城市震感明显，该地震是2023年1—9月继土耳其7.8级地震之后发生的第二大地震灾难。

摩洛哥位于非洲西北部，与欧洲接壤，属于非洲板块和欧亚板块挤压碰撞影响区域。摩洛哥西部横亘着的亚特拉斯(Atlas)山脉一直延伸到突尼斯，全长约2000km，是非洲板块和欧亚板块交汇形成的大陆内造山带，从晚新生代至今经历了复杂而漫长的地质构造历史(Timoulali et al，2015)。亚特拉斯山脉东部地表较为平缓，西部逐渐升高(Teixell et al，2005)，本次地震就发生在西南段山脉。长周期地震折射剖面在该区域的地球物理调查结果显示，亚特拉斯地区地壳厚度在EW向上变化剧烈，在25~39km之间变化(Wigger et al，1992)，西部区域莫霍面深度可达35~44km(Ghomsi et al，2020)。

1920年以来，摩洛哥及其附近区域共发生3.0级以上地震1900多次，其中5.0级以上地震42次(图 1)。1960年2月29日摩洛哥阿加迪尔(Agadir)M_W5.9地震由于震源深度较浅(2~3km)，虽然地震强度为中等，但仍造成约1.2万人丧生，75%的建筑损毁，地表最大烈度达到Ⅸ度(Timoulali et al，2011)。本次地震发生在1960年地震的EN向约140km处，是百年来摩洛哥历史最大地震，本次地震释放能量约为1960年地震的三十多倍。

破坏性强震发生后，快速且准确地确定震源破裂能量时空分布特征和矩震级大小，对于精准评估受灾区域大小及灾情严重程度、指导应急救援至关重要。区别于广泛应用的有限断层反演方法(许力生等，2002；张勇等，2010)，基于远场P波数据的台阵反投影方法不涉及震源参数，较少依赖先验经验，无需计算格林函数等，具有计算效率高、结果直接而稳定的优点(Wang et al，2017；曹博男等，2021；Chen et al，2023)。2004年12月26日苏门答腊—安达曼M_W9.1地震发生后，Ishii等(2007)和Krüger等(2005)开始将台阵反投影技术应用于远震震源过程研究，得到了2004年苏门答腊—安达曼地震能量最大值的位置和对应的破裂时刻。此后，反投影方法逐步广泛应用于破坏性大震的破裂过程成像和减灾中(Koper et al，2011；刘志鹏等，2018；Yao et al，2019；Chen et al，2022)。矩震级是完全由地震矩决定的震级标度，是评估地震潜在破坏的最佳震级(Hanks et al，1979)。在Hara(2007)研究的基础上，Wang等(2017)逐步发展出基于反投影计算的矩震级估算方法，将震源破裂持续时间和P波振幅与震级联系起来，具有计算速度快且结果稳定的优势。

本文基于美国阿拉斯台阵记录到的宽频带远场P波数据，利用远场台阵反投影方法分析本次摩洛哥地震震源破裂空间和时间分布特征，并结合全球GSN台网波形数据获取P波最大振幅，估算地震的矩震级。

1 研究方法与数据处理

本文利用Wang等(2017)发展的基于远场P波数据的台阵反投影方法估算破坏性强震震源破裂过程和矩震级，计算过程可分为三步。

首先，根据方法适用震中距范围(30°~85°)和台站分布特征，选取美国阿拉斯加台阵震中距74°~87°远场P波垂向分量数据(图 2)，进行反投影成像计算。对数据进行去均值、去线性趋势等预处理。为尽可能减小地下结构的横向不均匀性对成像位置的影响，选取台阵中距离几何中心最近的台站作为参考台站，将台阵波形按照参考台站P波到时对齐并利用两步互相关计算获取台阵校正量。对原始数据先进行0.05~1.0Hz滤波处理，提取波形中的低频成分，进行波形粗对齐。然后在低频波形粗对齐的基础上，对原始数据进行0.5~2.0Hz带通滤波处理，对波形的高频成分进行精确对齐，对齐时窗为10s。仅保留与参考台站波形相关系数大于0.7的波形数据参与后续计算。

经过波形相似度筛选后，保留122条波形，以参考台站AK.BPAW的直达P波理论到时122.5s为标准，将高频波形数据经互相关对齐后可以看出(图 3)，不同台站之间波形的相似性极高，存在几乎完全对齐的同向轴，可以开展后续反投影叠加计算。

然后，根据震中位置和震源深度，以及震级和破裂尺度之间的经验关系(Wells et al，1994)，在震源深度处将震源区划分为2km×2km等间隔的网格，共45×45个格点。使用0.5~2Hz高频数据进行反投影叠加成像，选取10s长的滑动时窗，步长为1s，对数据进行分段截取和叠加成像。当叠加能量达到总能量的90%，或叠加能量达到总能量的80%且后续能量振幅小于最大能量的0.1倍时，选取两个持续时间中较短的时间定义为地震破裂持续时间。

最后，从IRIS下载震中距85°以内的GSN台站垂向分量波形数据。利用在全球分布较为均匀的GSN台网垂向分量波形数据，共取得61个台站P波波列最大位移数据。对数据进行去均值、去线性趋势和去仪器响应处理，根据PREM模型估算理论P和S波到达时间，获取P波最大位移数据。结合反投影计算过程得到的震源破裂持续时间，估算地震矩震级M_dt，即

$ M_{\mathrm{d} t}=\frac{2}{3} \frac{1}{N} \sum\limits_{i=1}^N\left[\lg \left(|u|_{(\max)}^P \cdot R \cdot T\right)\right]+\delta(h, k) $

(1)

其中，|u|_(max)^P为第i个台站P波最大位移，R为对应震中距，T为第二步反投影成像得到的震源破裂持续时间，N为总的台站数。δ(h，k)为常数项，与震源时间函数的形状和震源深度有关，即

$ \delta(h, k)=\frac{2}{3}\left[\lg \left(k \cdot 4 {\rm{ \mathsf{ π} }} \rho_s^{\frac{1}{2}} \rho_r^{\frac{1}{2}} \alpha_s^{\frac{5}{2}} \alpha_r^{\frac{1}{2}}\right)-9.1\right] $

(2)

其中，k值与震源时间函数的形状有关，其最佳拟合值为0.48，ρ_s、ρ_r、α_s、α_r分别为震源深度处和台站位置处的密度和P波速度(基于PREM模型)。

2 结果与讨论

通过上述计算，对0.5~2Hz高频数据进行反投影成像分析，得到一系列能量辐射子源的时空分布特征。反投影空间成像结果中不同颜色的圆圈代表不同时刻能量极大值的空间分布位置，圆圈越大表示释放的能量越多(图 4(a))。从空间上看，能量辐射子源集中在震中约10km范围内，约22s之后在震中的西北部区域存在一些小的能量释放点，考虑到地震震级与破裂时间之间的经验关系(Yao et al，2019)，以及破裂方向不大可能出现突然的转向，筛去传播距离较大且分布散乱的破裂子源。从时间上看，震中EN方向能量释放整体较早，WS方向较晚，空间位置上具有较好的连续性。

中国地震台网中心利用W震相反演方法(Kanamori et al，2008)，获得了本次地震的震源机制解，节面Ⅰ：128°/25°/136°，节面Ⅱ：260°/73°/72°，显示出逆冲叠加走滑的性质。USGS(2023)利用长周期远震数据波形拟合方法得到的震源破裂过程显示破裂整体长度约15km，断层面的走向和倾角选择了节面Ⅱ结果(249°/65°)，破裂中心位于震中偏EN方向(图 4)。震中区域地质构造图显示，本次地震发生在亚特拉斯山脉的西南段(Sébrier et al，2006)，震中所在区域邻近山脉北部一个向南倾斜的逆冲断层，余震数量较少但基本分布在主震的WS方向(图 1)。综合上述结果，本次地震整体表现出NE-SW走向的破裂特征。

地震破裂过程中能量随时间的变化规律较为简单(图 5)，仅在10~20s内存在一个归一化能量大于0.8的峰值，表明地震能量主要在破裂开始后的10~20s内释放。根据能量释放比例计算得到的震源破裂持续时间约28s，与能量空间分布结果基本一致，这一数值将参与后续矩震级估算。

震级是衡量地震强度大小的度量，震级的大小直接与灾害评估和应急处置相关(杨晶琼等，2019)。利用在全球分布较为均匀的GSN台网垂向分量波形数据，共取得61个台站P波波列最大位移数据(图 6)。随震中距增大，P波最大位移数值逐渐减小。基于反投影叠加成像计算得到的震源破裂持续时间以及P波最大位置数据，计算得到的矩震级M_dt为7.0。该数值与中国地震台网中心发布的矩震级(M_W6.9)、USGS测定矩震级(M_W6.85)和GCMT测定的矩震级(M_W6.9)较为吻合，证实了该震级测定方法针对中强震的有效性。

计算量小、自动化程度高是基于反投影震源破裂过程和矩震级估算方法的最明显优势。利用该方法可在震后约15min快速获取地震破裂的时空分布和震级信息，且地震震级越大，破裂过程越复杂，该方法的优势越明显(Wang et al，2017)。本研究结果表明，对于震级低于7级的中强震，该方法能取得稳定可靠的结果，其应用范围未来可以进一步扩大。在时效性上，本次地震采用了震中距74°~87°的阿拉斯加台阵数据，如果存在距离更近的区域密集地震台阵数据，则可在约10min获得稳定可靠的矩震级结果。

3 结论

本文基于全球一维速度模型(PREM)，使用美国阿拉斯加台阵记录到的远场P波数据，对2023年9月9日摩洛哥M_W6.9地震震源破裂过程时空分布特征进行反投影分析，并结合GSN台网P波最大位移数据估算了该地震的矩震级。结果显示，本次地震整体呈现出NE-SW走向的破裂特征，在10~20s地震能量集中释放，且集中在震中约10km范围内，矩震级M_dt为7.0。结合震源机制解和震中区域地质构造特征，地震能量辐射子源主要分布在亚特拉斯山脉西南段北部南向倾斜逆的冲断层上，与主震后余震的空间分布特征基本一致，表明本次地震可能发生在高倾角WS向的逆冲走滑断层上。

致谢: 感谢王墩教授为本文提供台阵反投影和矩震级计算程序，感谢审稿专家提出的宝贵修改意见。