用矩阵分解方法识别地下核爆炸

赵克常 张献兵

北京大学地球与空间科学学院, 北京 100871; †通信作者, E-mail: xbzhang@pku.edu.cn

摘要 利用 3 次地下核试验及其附近 3 次天然地震的地震资料, 反演获得震源地震矩张量, 通过矩阵分解方法识别地下核爆炸。结果表明, 地下核爆炸地震的震源中有较显著的爆炸源(EXP), 也有补偿线性矢量偶极(CLVD)源和双力偶(DC)源。CLVD 的物理机理是爆炸引起介质破裂, 在地下核爆炸震源中占较大的比例。与地下核爆炸相比, 天然地震一般为剪切位错模式, 其震源中DC占较大的比例。

关键词 地下核爆炸; 矩张量; 爆炸源(EXP); 补偿线性矢量偶极(CLVD)源; 双力偶(DC)源

随着全球核试验不断增加, 核查工作显得越来越重要。地下核爆炸是一种人为的地震事件, 与天然地震的本质区别在于其震源特性的差异。通常, 核试验是秘密进行的, 其地震资料往往难以获得。因此, 如何利用较少的资料来识别地下核爆炸, 是地震学界面对的重要课题。

在众多地震学方法中, 地震震源矩张量是识别核爆炸地震、陷落地震和构造地震的有力工具之一[1]。Liu 等[2]和 Wang 等[3]认为, 核爆炸震源主要为各向同性源。Ford 等[4]和 Wang 等[5]认为, 固体矿产资源开采过程中易发生的陷落地震主要表现为闭合源。绝大多数天然地震是构造地震, 由断层错动形成, 主要为双力偶(double couple, DC)震源。利用区域地震波形资料, 可以反演得到震源矩张量, 从而推断震源的性质。Dreger 等[6]对 3 次地下核爆炸进行震源矩张量反演, 结果显示, 在震源机制中DC 分量占比很小, 不具备天然地震的特性。Ford等[4]根据矩张量, 识别出数次核爆炸地震、陷落地震及天然地震。Han 等[7]和 Liu 等[2]甚至识别出一个地震序列的“主震”和两个“余震”的震源性质, “主震”是地下核爆炸, “余震”则是由核爆炸引发的陷落地震和构造地震。Bowers 等[8]将震源矩张量分解成各向同性源、双力偶源和补偿线性矢量偶极(com-pensate linear vector dipole, CLVD)源 3 个矩张量, 用于确定震源的标量矩。Ford 等[4]用分解地震矩张量特征值的方法, 分析震源的属性。

用震源矩张量表示震源特征, 能够将震源效应与传播路径的效应区分开来, 将资料、震源和传播路径三者之间归结为一种线性关系。如果已知震源位置和地球介质的结构, 由给定的矩张量就可以线性地正演出位移场的分布。反之, 如果已知震源位置和相应介质模型下的格林函数, 就可以由记录资料线性地反演出地震矩张量。位移场采用矩张量中各分量作为权重系数的基本格林函数的线性组合来表示, 就可以利用面波记录资料反演震源矩张量。

反演天然地震的研究以 Jost 等[9]的工作为代表。Herrmann 等[10]在前人工作的基础上, 给出任意一个位错点源和爆炸源(explosion, EXP)产生的地震波场的表达式及所需的 10 个格林函数, 可以从理论上反演含有爆炸源成分的全矩张量(full moment tensor)。

本文采用 Chen[11]和 Yao 等[12]基于广义反射系数和透射系数的水平层状地球模型中理论地震图的计算方法, 反演前苏联在东哈萨克斯坦地下核试验场进行的地下核爆炸的震源矩张量。在求得震源矩张量的基础上, 利用矩阵分解方法, 探索识别地下核爆炸的依据。

1 方法原理

1.1 矩张量的定义

地震震源矩张量通常指在点源近似的条件下, 地震震源等效力完整描述的一级近似。Gilbert[13]首先将矩张量的概念应用到地球表面位移的计算中, 将位移表示为矩张量与描述震源与接收器之间传播路径脉冲响应的格林函数的卷积。根据震源矩张量、格林函数和观测资料之间的线性关系, 由台站观测资料可以计算得到矩张量解。本文即采用这种反演方法, 求解地震震源矩张量。

矩张量 Mpq 是二阶张量, 具有应力张量的所有性质, 点矩张量在 t 时刻及 x 场点产生的位移场可以表示成下列卷积形式:

width=80.25,height=16.5。(1)

Gip 是震源处 p 方向单位脉冲力引起的格林函数, 代表在场点处产生的位移的 i 方向分量, Gip,q 则是 Gipq 方向震源点坐标的微商(偏导数)。若已知格林函数的表达式, 则位移场可用矩张量 Mpq 表示。

1.2 地震矩张量的表述方法

地震矩张量 Mij 是一个对称张量, 因此有 6 个独立分量。利用天然地震断层面的倾角(δ)、滑动角(λ)和走向(ϕf)计算各分量的公式[14]如下:

width=210.7,height=17.25(2)

width=216.8,height=17.25(3)

width=162,height=16.5(4)

width=214.5,height=27(5)

width=201,height=16.5(6)

width=201,height=16.5(7)

其中, Mij = Mji (i, j = x, y, z; ij)。

爆炸与天然地震的震源具有完全不同的矩张量形式。通过求解未知地震事件的震源矩张量, 并基于结果进行差异分析, 就可以进行震源特性的研究, 达到识别地下核爆炸的目的。

1.3 主要震源模式的矩阵表述

从理论上讲, 地下核试验的震源机制是球对称的爆炸机制, 爆炸能量会朝各个方向辐射。但是, 实际情况并非如此。由于地壳结构是非均匀的, 并且各个方向的构造不同, 因此地下核爆炸的震源不会出现理想的对称形式。Day 等[15]提出的层裂辅助震源模型, 在研究地下核爆炸震源的物理机制时发挥了重要作用。

如图 1 所示, 在垂直方向, 地下核爆炸产生的爆炸冲击波从爆炸中心向上传播, 经地表面反射后, 由压缩波转变成拉伸波。该拉伸波与继续向上传播的压缩波相互叠加, 导致爆炸源上方的地表岩层破裂。该过程可以用张裂源(tension crack, TC)描述。由于爆炸源上方的地表岩层破裂, 引起应力场重新分布, 导致介质的构造应力释放。此释放过程可以用逆倾滑断层模式描述。此震源形状为顶点在爆炸点的倒立的圆锥体, 锥体底面位于层裂面至地表面。此应力释放过程可以用 CLVD 源[15‒17]来描述。

对一次地震事件而言, 可以把震源分解成为爆炸(EXP)源、补偿线性矢量偶极(CLVD)源和双力偶(DC)源。3 种形式的震源矩张量[18‒19]表示如下。

width=141.7,height=85.05

图1 层裂模型示意图

Fig. 1 Diagram of spalling model

EXP 源的震源矩张量形式为

width=88.1,height=46.2(8)

其中, M0 = ΔV(λ + 2/3μ)为爆炸源矩, ΔV 为球对称的爆炸源引起的体积变化, λμ为拉梅常数。

CLVD 源的震源矩张量形式为

width=127.9,height=46.2(9)

其中, MCLVD 为补偿线性矢量偶极源矩。

DC 源的震源矩张量形式为

width=101,height=46.2(10)

其中, MDC 为双力偶源矩。双力偶(DC)描述的是各项同性介质中的沿断层面发生的剪切滑动, 与剪切位错是等效的。通常认为, 双力偶是天然地震震源的力学模型。

1.4 用矩阵分解方法识别地下核爆炸

地震震源矩张量是实对称矩阵, 可以求出其对角化的相似矩阵:

width=145.05,height=46.2。(11)

按照 EXP, DC 和 CLVD 的基本震源模式, 将式(11)右边的对角化矩阵展开如下:

width=228.9,height=93.5

width=207.4,height=46.2(12)

于是, 从地震矩张量的角度, 将地下核爆炸的鉴定问题归结为式(12)中 3 个系数width=50.5,height=27.95, −λ1 +λ2width=67.7,height=27.95相对比例的计算问题, 据此可以得到震源中 EXP, CLVD 和 DC 的准确比例, 从而达到识别天然地震和地下核爆炸的目的。

从理论上讲, 对于地下核爆炸地震, EXP 的占比应很大, 甚至接近 100%, CLVD 也占一定的比例, 与其产生的物理机理有关, DC 的占比应为 0; 对于天然地震, DC 几乎占 100%, EXP 和 CLVD 的占比应接近 0。当然, 介质模型不可能与实际情况完全符合, 加上计算方法和数据处理会有误差。所以, 无论是核爆炸地震还是天然地震中, EXP, CLVD 和DC 这 3 种震源都有, 识别依据如下: 核爆炸震源中EXP 比较明显, 且 CLVD 占比较大; 天然地震震源中则 EXP 占比很小, DC 占比很大。

2 实际应用

本文所用前苏联在东哈萨克斯坦核试验场进行的 3 次地下核试验以及试验场附近发生的 3 次天然地震事件的台站资料, 来自美国地震学研究联合会(Incorporated Research Institutions for Seismology,IRIS)发布的乌鲁木齐台站(WMQ station)的记录。采用适合东哈萨克斯坦地区的 Steven 模型[20‒21]计算理论格林函数。把式(1)离散化, 用以矩张量分量为权重系数的 10 个格林函数的线性组合表示位移。首先通过反演得到 6 次事件的震源矩张量, 然后对二类震源矩张量进行矩阵分解, 并对结果进行分析。为了易于区分, 把 3 次地下核试验记为事件1、事件 2 和事件 3, 把 3 次天然地震记为事件 4、事件 5 和事件 6。

2.1 地下核爆炸震源矩张量的矩阵分解及识别

首先分析事件 1。经过计算, 该地下核试验地震的震源矩张量如式(13)所示:

width=133.25,height=46.2(13)

其中, M0 是事件 1 的震源矩(M0 = 1017Nm)。对与式(13)右侧矩阵相似的对角阵进行矩阵分解, 如式(14)所示:

width=168.7,height=46.2

width=167.65,height=46.2(14)

接着分析事件 2。经过计算, 该地下核试验地震的震源矩张量如式(15)所示:

width=133.25,height=46.2(15)

其中, M0 是事件 2 的震源矩(M0 = 1017 Nm)。对与式(15)右侧矩阵相似的对角阵进行矩阵分解, 如式(16)所示:

width=173,height=93.5(16)

最后分析事件 3。经过计算, 该地下核试验地震的震源矩张量如式(17)所示:

width=135.4,height=46.2,(17)

其中, M0 是事件 3 的震源矩(M0 = 1017 Nm)。对与式(17)右侧矩阵相似的对角阵进行矩阵分解, 如式(18)所示:

width=171.95,height=93.5(18)

事件 1、事件 2 和事件 3 震源矩张量的矩阵分解结果如表 1 所示。可以看出, 地下核试验的震源矩张量地矩阵分解结果中, 同时包含 EXP 源、DC源和 CLVD 源。虽然 EXP 占比较大, 但 CLVD 的占比远大于 EXP 源, DC 源占比最小。这种现象是由震源机理决定的, 震源附近区域地下介质的复杂性、爆炸导致的层裂[22]等因素, 使得核爆炸不会显示理想的对称性。然而, 天然地震的震源机制主要为 DC 源, 所以完全可以将 DC 源占比小作为识别核爆炸的标志。图 2 是 3 次地下核爆炸的震源机制解, 明显具有地下核爆炸地震的特征。

表1 地下核试验的震源矩张量的矩阵分解结果

Table 1 Matrix decomposition of source moment sensor for an test of underground nuclear

事件编号占比/% EXPDCCLVD 124.8910.2264.89 222.7113.5663.73 322.9712.9864.05

2.2 天然地震震源矩张量的矩阵分解及识别

同样, 可以得到上述地核爆炸附近 3 次天然地震(相关信息见表 2)的震源矩张量分解结果。

首先分析事件 4。该天然地震的震源矩张量如式(19)所示:

width=133.25,height=46.2(19)

其中, M0 是事件 4 的震源矩(M0 = 6.91 × 1016 Nm)。对与式(19)右侧矩阵相似的对角阵进行矩阵分解, 如式(20)所示:

width=169.8,height=93.5(20)

width=221.15,height=68

(a)事件1; (b)事件2; (c)事件3

图2 3次地下核爆炸震源机制解

Fig. 2 Focal mechanism solution of three underground nuclear tests

表2 3次天然地震的相关信息

Table 2 Information about three nature earthquakes

事件编号发生日期震中位置震源深度/km震级/mb 41990‒10‒2443.79°N 83.99°E15.05.4 51995‒05‒0243.49°N 84.54°E18.05.5 62011‒05‒0143.72°N 77.73°E24.65.3

接着分析事件 5。该天然地震的震源矩张量如式(21)所示:

width=133.25,height=46.2(21)

其中, M0是事件5的震源矩(M0 = 2.24 × 1017 Nm)。对与式(21)右侧矩阵相似的对角阵进行矩阵分解, 如式(22)所示:

width=169.8,height=93.5(22)

最后分析事件6。该天然地震的震源矩张量如式(23)所示:

width=138.65,height=46.2(23)

其中, M0是事件 6 的震源矩(M0 = 6.47 × 1016 Nm)。对与式(23)右侧矩阵相似的对角阵进行矩阵分解, 如式(24)所示:

width=167.65,height=93.5(24)

事件 4、事件 5 和事件 6 震源矩张量的矩阵分解结果如表 3 所示。可以看出, 与 2.1 节中地下核试验的震源矩张量的分解结果截然不同, 天然地震的震源中 DC 占绝对优势(80%以上), 这是因为天然地震多为断层地震, 其震源一般是位错源; EXP 源的比例小到可以忽略; CLVD 源占有一定的比例, 这是可以理解的, 地球介质的各向非均匀性、实施核试验时的刻意干扰以及数据处理方法和计算带来的误差等因素都会导致 CLVD 的存在, 但占比不大。因此, 依然可以清晰地区分地下核爆炸和天然地震。图 3 是 3 次天然地震的震源机制解, 与图 2所示地下核爆炸的震源机制解截然不同, 具有鲜明的断层地震特征。

表3 3次天然地震的震源矩张量的矩阵分解结果

Table 3 Matrix decomposition of source moment sensor for three nature earthquakes

事件编号占比/% EXPDCCLVD 40.0193.62 6.37 50.0099.90 0.10 60.0482.6617.30

width=221.15,height=68

(a)事件4; (b)事件5; (c)事件6

图3 3次天然地震震源机制解

Fig. 3 Focal mechanism solution of three nature earthquakes

3 结论与讨论

本文对 3 次地下核试验及其附近 3 次天然地震的震源矩张量进行矩阵分解, 结果表明, 地下核爆炸的震源中有较明显的爆炸源, 补偿线性矢量偶极源占比大于 60%; 天然地震震源中双力偶源占 80%以上, 爆炸源和补偿线性矢量偶极源占比较小。这些特征可以用来清晰地区分地下核爆炸与天然地震事件, 从而达到识别地下核爆炸的目的。本文的研究结果对于从震源角度了解地下核爆炸的物理机制具有参考意义, 同时也为利用区域少量甚至是单站地震数据来识别地下核试验提供了理论依据。

然而, 由于诸多因素的制约, 本文得出的结论在识别地下核爆炸时会存在精度瓶颈。随着核爆炸鉴别技术的发展, 逃避核查的手段也在不断进步, 目前主要有以下 3 种方法。1)解耦: 如果把核爆炸放在一个足够大的空腔中进行, 爆炸产生的地震波就会小很多。这种方法非常有效, 可以使地震波的振幅减少 1~2 个数量级。然而, 解耦的方法对小当量的核爆炸适用, 对大当量的核爆炸则无法实施。2)多次爆炸: 将大当量的核爆炸试验紧接在较小当量的核爆炸试验之后进行, 就容易被人误读成体波和面波, 并且, 多次爆炸的波形叠加, 导致波形复杂化, 使得多次核爆炸的地震记录更像一次天然地震的记录, 也容易使人们误解成一个地震序列。3)利用天然地震来掩盖核爆炸事件: 在一个中等地震发生后, 立即实施核爆炸试验, 把核爆炸信号掩藏到地震信号当中。

识别核爆炸与逃避核查的方法是矛盾的统一体, 互相促进, 交替上升。逃避核查技术的效果与核爆炸鉴别方法有着紧密的联系, 既不存在一种确凿的核爆炸识别方法, 也不存在一种不留痕迹的逃避核查方法。随着包含横向不均匀性的地球三维模型在格林函数计算中的引入, 对地下核爆炸的识别会更加精准。

参考文献

[1]林鑫, 王向腾, 赵连锋, 等. 核试验监测的地震学研究综述.地球物理学报, 2019, 62(11): 4047‒4066

[2]Liu J Q, Li L, Zahradnik J, et al. North Korea’s 2017 test and its nontectonic aftershock. Geophysical Re-search Letters, 2018, 45(7): 3017‒3025

[3]Wang T, Shi Q B, Nikkhoo M, et al. The rise, collapse, and compaction of Mt. Mantap from the 3 September 2017, North Korean nuclear test. Science, 2018, 361: 166‒170

[4]Ford S R, Dreger D S, Walter W R. Source analysis of the memorial day explosion, Kimchaek, North Korea. Geophysical Research Letters, 2009, 36(21): L21304

[5]Wang X T, Wang S F, Li Z W, et al. Source charac-terization of some collapse earthquakes due mining activities in Shandong and Beijing, North China. Seis-mological Research Letters, 2018, 90(1): 183‒193

[6]Dreger D, Woods B. Regional distance seismic mo-ment tensors of nuclear explosions. Tectonophy-sics. 2002, 356: 139‒156

[7]Han L B, Wu Z L, Jiang C S, et al. Properties of three seismic events in September 2017 in the northern Korean Peninsula from moment tensor inversion. Sci Bull, 2017, 62(22): 1569‒1571

[8]Bowers D, Hudson J A. Defining the scalar moment of a seismic source with a general moment tensor. Bull Seism Soc Am, 1999, 89(5): 1390‒1394

[9]Jost M L, Herrmann R B. A student’s guide to and review of moment tensor. Seism Res Lett, 1989, 60 (2): 37‒57

[10]Herrmann R B, Wang C Y. A comparison of synthetic seismograms. Bull Seism Soc Am, 1985, 75: 41‒56

[11]Chen X F. A systematic and efficient method of computing normal modes for multilayered half-space. J Geophys Int, 1993, 115: 391‒409

[12]Yao Z X, Harkrider D G. A generalized reflection-transmission coefficient matrix and discrete wave-number method for synthetic seismograms. Bull Seism Soc Am, 1983, 73: 1685‒1699

[13]Gilbert F. Excitation of the normal modes of the earth by earthquake source. J Geophys Res, 1971, 22: 223‒ 226

[14]Aki K, Richards P G. Quantitative seismology, theo-ry and methods. San Francisco: Freeman W H and Company, 1980

[15]Day S M, McLaughlin K L. Seismic source repre-sentations for spall. Bull Seism Soc Am, 1991, 81: 191‒201

[16]何永锋, 赵克常, 张献兵. 地下核爆炸的主要非爆炸源机制. 地球物理学进展, 2010, 25(3): 789‒794

[17]Knopoff L, Randall M J. The compensated linear-vector: a possible mechanism for deep earthquakes. J Geophys Res, 1970, 75: 4957‒4963

[18]何永锋, 陈晓非, 何耀峰, 等. 地下爆炸 Rg 波低谷点激发机理. 地球物理学报, 2005, 48(3): 643‒648

[19]Bowers D, Hudson J A. Defining the scalar moment of a seismic source with a general moment tensor. Bull Seism Soc Am, 1999, 89(5): 1390‒1394

[20]何永锋, 赵克常, 姚国政, 等. 基于 Lg 波识别判据的特性分析. 地震学报, 2011, 33(6): 715‒722

[21]何永锋, 陈晓非, 张海明. 地下核爆炸Lg波的激发机制. 地球物理学报, 2005, 48(2): 367‒372

[22]何永锋, 陈晓非, 张海明. 层裂对区域震相Lg波的影响. 地震学报, 2005, 27(3): 309‒316

Distinguishing Underground Nuclear Test by Matrix Decomposition

ZHAO Kechang, ZHANG Xianbing

School of Earth and Space Sciences, Peking University, Beijing 100871; † Corresponding author, E-mail: xbzhang@pku.edu.cn

Abstract Using the seismic data of three underground nuclear tests and three nearby natural earthquakes, the focal seismic moment tensor is obtained by inversion, and the underground nuclear explosion is identified by matrix decomposition method. The results show that there are obvious explosion sources (EXP), compensated linear vector dipole (CLVD) source and double couple (DC) source in the source of underground nuclear explosion earthquake. The physical mechanism of CLVD is medium rupture caused by explosion, which accounts for a large proportion of underground nuclear explosion sources. Compared with underground nuclear explosion, natural earthquake is generally shear dislocation mode, and DC accounts for a large proportion of its source.

Key words underground nuclear explosion; seismic moment tensor; EXP source; CLVD source; DC source

doi: 10.13209/j.0479-8023.2022.042

国家自然科学基金(41374068)资助

收稿日期: 2021-05-21;

修回日期: 2021-06-29