地球物理学报  2020, Vol. 63 Issue (2): 562-572   PDF    
应力张量反演的遗传算法及其在青藏高原东北缘的应用
李振月1, 万永革2, 胡晓辉2, 李泽潇2, 杨帆2, 闫睿3     
1. 中国地震局地球物理研究所, 北京 100081;
2. 防灾科技学院, 河北 三河 065201;
3. 河南省地震局, 郑州 450016
摘要:为准确而快速地求解区域构造应力张量,构建了基于震源机制解反演应力张量的遗传算法策略,详细介绍了计算原理,分析了以断层面上剪应力方向和滑动方向的偏差为约束反演应力张量而忽略剪应力大小偏差对反演结果的影响,表明仅考虑剪应力方向和滑动方向的偏差就可以得到正确结果.利用不同应力状态下人工合成的包含不同噪声水平的震源机制数据对该方法进行检验,并与网格搜索法得到的结果比较.表明本文方法所得结果的拟合差均小于网格搜索法结果的拟合差,并且二者相差较小,体现了本文方法的稳健性.将该方法应用于青藏高原东北缘(103-106°E,34.5-37.5°N)应力张量的估计,结果显示,该区域主要受控于青藏高原近东西向的挤压,从而导致阿拉善块体以及华南块体方向的拉张,并且向华南块体方向的拉张作用强于向阿拉善块体的拉张作用.
关键词: 震源机制解      应力张量      遗传算法      青藏高原东北缘     
A genetic algorithm for stress tensor inversion and its application to the northeast margin of the Tibetan Plateau
LI ZhenYue1, WAN YongGe2, HU XiaoHui2, LI ZeXiao2, YANG Fan2, YAN Rui3     
1. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
2. Institute of Disaster Prevention, Sanhe Hebei 065201, China;
3. Henan Earthquake Agency, Zhengzhou 450016, China
Abstract: In order to obtain regional tectonic stress tensors accurately and quickly, the strategy of stress tensor inversion using the genetic algorithm based on focal mechanisms is constructed. Its principle is introduced in detail. We analyze the influence on the inversion results under constraint of direction deviation between shear stress and slippage on the fault plane while neglecting the shear stress magnitude deviation. The result indicates that we can obtain the correct result only under constraint of direction deviation between shear stress and slippage. The method is tested using synthetic data produced by different stress states with different noisy. We compare the result with another result obtained by a grid search algorithm. It indicates that the misfit of results obtained by the genetic algorithm are all smaller than those of results obtained by the grid search algorithm. And the difference of misfit obtained by the two method is small, which reflects the robustness of the method in this paper. We apply this method to the northeast margin of the Tibetan Plateau (103-106°E, 34.5-37.5°N). The result shows that this region is mainly controlled by nearly east-west compression from the Tibetan Plateau, leading to extension towards the Alxa block and South China block. And the extension towards South China block is stronger than that towards Alxa block.
Keywords: Focal mechanisms    Stress tensor    Genetic algorithm    Northeast margin of Tibetan Plateau    
0 引言

随着地震台站的加密和监测能力的提高,积累了越来越多的地震记录,利用地震记录得到的震源机制解描述了震源断层的错动方式,为研究地壳构造演化的动力来源提供了很好的数据基础(许忠淮,2001Wan,2010Christova,2015).目前,利用震源机制解反演应力张量的方法可以分为两类,一类是力轴张量法(钟继茂和程万正,2006盛书中等,2013);另一类是滑动方向拟合法(Gephart and Forsyth, 1984许忠淮,1985Michael,1987Wan et al., 2016).力轴张量法通过把地震释放的矩张量叠加得到应力张量;滑动方向拟合法是拟合断层面的滑动方向和最大剪应力的方向,使其偏差最小作为约束,求取相应应力张量.不同的是,力轴张量法是将地震释放的矩张量作为应力张量,滑动方向拟合法求取的是导致地震发生的应力张量.本文研究的主要目标是改进滑动方向拟合法.滑动方向拟合法反演应力张量需要明确知道真实的发震断层面,然而震源机制解给出了震源断层两个等可能的节面,究竟选择哪个节面作为“真实断层面”来反演应力张量?不同的学者给出了不同的方法.Gephart和Forsyth(1984)提出网格搜索法,在全空间范围内搜索上述约束下的最优解,该方法通过使两个节面绕任意轴旋转,当给定应力状态下节面的剪应力方向与滑动方向一致时,旋转角度较小的节面作为“真实断层面”;许忠淮(1985)通过穷举所有断层面的组合,将每个组合都进行试算,以确定最优的应力张量;Michael(1987)将反演问题线性化,利用bootstrap方法随机抽样随机选择地震和节面进行反演,统计反演结果在一定置信区间下的集中范围,从而确定最优解和不确定范围;Wan等(2016)提出考虑震源机制权重的网格搜索法,给定应力状态下,将两个节面中滑动方向与剪应力方向差别较小的节面作为“真实断层面”.

有两个问题值得思考:(1)假定研究区内有m次地震,若对所有可能的断层面组合一一试算,得到最优应力张量共需要计算2m次,当m较大时,由于计算量太大,计算过程缓慢.若不对所有的组合一一计算的话又可能由于“真实断层面”的选择不当导致反演结果不可信,能否兼顾反演结果的可靠性与计算耗时?(2)反演应力张量都是拟合断层面的滑动方向和剪应力的方向,未考虑剪应力大小偏差对反演结果的影响.断层面上剪应力大小的偏差在反演应力张量问题中起什么样的作用,它能否被忽略?

本文研究基于这两个问题展开.针对问题(1),为保证应力张量反演结果的可靠性,采用Wan等(2016)网格搜索法中的方法:给定应力状态下,对两个节面都进行试算,选择偏差较小的节面作为“真实断层面”,此方法与计算2m次断层面的组合是等效的.至此,把问题聚焦于减少计算耗时.本文将近来风靡全球的人工智能技术引入——使用遗传算法反演应力张量加快计算速度.将人工智能与地球物理相结合,前人已经做了一些探索(Kobayashi and Nakanishi, 1994石耀霖和金文,1995万永革和李鸿吉,1995Wu et al., 2008),而今,随着人工智能技术的飞速发展,将其优势应用于地球物理势在必行,这也是本文的一大尝试.对于问题(2),下文(第2节)将进行详细分析.

1 方法

震源机制解给出了震源断层的滑动方向,仅根据滑动方向不能确定应力张量各分量的绝对大小,只能得到各分量的相对大小.因此,必须对应力张量进行简化.简化过程见附录,简化后的应力张量包含四个参数:ϕδλR,取值范围分别为:ϕ∈[0, 2π]、δ∈[0, π/2]、λ∈[0, 2π]、R∈[0, 1]. , 称为应力形因子(Angelier, 1979, 1984Gephart and Forsyth, 1984Michael,1984Wan et al., 2016),σ1σ2σ3分别为最大、中间、最小主应力.符号规定与弹性力学(徐芝纶,2006)一致:拉张为正,挤压为负.R描述了构造应力的状态:当R=0,σ2=σ3,对应于单轴拉张的应力状态;当R=0.5,中间应力轴(σ2)不表现构造作用;当R=1,σ1=σ2,对应于单轴挤压的应力状态(万永革等,2011).

下面首先阐述根据应力张量计算的理论滑动方向和剪应力大小(此处得到的是相对剪应力,方便起见,本文称相对剪应力为剪应力)与震源机制解给出的滑动方向和剪应力大小之间的偏差得到应力张量的拟合差,其次介绍利用遗传算法全空间搜索与震源机制解最匹配,即拟合差最小的应力张量.

1.1 应力张量的拟合差

断层面上剪应力方向与滑动方向之间的角度偏差:

(1)

式中σs为根据应力张量计算的剪应力矢量(理论滑动方向矢量),s为观测滑动方向矢量.

断层面上理论剪应力大小与实际剪应力大小的偏差:

(2)

式中p为断层面的应力矢量.

引言中叙及的方法(Gephart and Forsyth, 1984许忠淮,1985Michael,1987Wan et al., 2016)都是对(1)式角度偏差拟合,略去了(2)式剪应力大小的偏差,Thakur等(2017)考虑了剪应力大小的偏差,将偏差函数定义为

(3)

其中λ∈[0, 1],是剪应力方向与滑动方向角度偏差所占的权重.遗憾的是,文中未给出λ的取值会如何影响反演结果,角度偏差和剪应力大小偏差在应力张量反演中扮演的角色究竟孰轻孰重?针对这一问题,下文(第2节)将进行进一步的分析.

双力偶震源机制模型给出两个等可能的地震断层节面,这两个节面的偏差(f1f2f3)不同.分别计算这两个节面的偏差,将偏差较小的节面作为发震断层面,其偏差作为应力张量对于该震源机制解的偏差,这样可以避免由于断层面选择不当导致应力张量反演结果不可信.简化的应力张量共有4个待求参数,所以至少需要4个震源机制解数据求解某一区域的应力张量,震源机制解数据中包含的噪声越大,就需要越多的震源机制解数据拟合应力张量的四个参数.将某应力张量对所有震源机制解的平均偏差作为该应力张量的拟合差(4式),下一步的目标是利用遗传算法寻求拟合差最小的应力张量:

(4)

式中,Fi是第i种偏差约束(1、2或3式)下应力张量的拟合差,fij表示第i种偏差约束计算得到的第j个震源机制解的偏差,N是震源机制解个数.

1.2 遗传算法反演应力张量

同其他应力张量反演方法(Gephart and Forsyth, 1984Michael,1984Michael,1987许忠淮,1985Wan et al., 2016)一样,本文应力张量反演基于以下假设:

(1) 研究区域内应力是均匀的;

(2) 地震发生在先存断层上;

(3) 断层沿着剪应力的方向滑动;

(4) 断层之间的滑动互不影响.

遗传算法反演应力张量基于一定规模的种群,循环对其进行选择、交叉、变异操作,直到种群收敛到一定程度或者进化达到一定代数终止循环,将最优个体作为结果输出.

(1) 生成初始种群

随机生成一定规模的初始种群,种群中每个个体都是由第1节所述的四个参数(ϕδλR)组成,每个参数由8个基因(可根据精度要求调整)二进制编码,单个个体如图 1所示.

图 1 二进制编码个体 Fig. 1 Binary coded individual

将二进制的个体转换到相应的区间:

(5)

其中,pop10是二进制参数转化为十进制的值,parametermax是参数的最大值.

(2) 选择

个体是以一定的概率被选择的,被选择的概率取决于其适应度,个体的适应度定义如下:

(6)

其中,Fi是(4)式得到的拟合差.

个体被选择的概率为

(7)

结合(7)式可以看出,在(6)式分母加0.1,是为了避免拟合差很小的个体适应度太大,使得种群过早收敛,这样极有可能使搜索陷入局部最优.

(3) 交叉和变异

交叉和变异是种群中引入新个体不可或缺的手段,期望通过交叉和变异得到更优秀的个体,这方面有一定的相似性.交叉和变异都是以一定概率进行的,若概率太大,可能会导致已有的优秀个体被破坏,若概率太小,种群则不能有效进化,过早收敛,使得搜索陷入局部最优.不同之处在于,交叉是通过两个个体交换一部分基因产生新个体,变异是通过改变个体的某个基因产生新个体(如图 2所示),所以交叉产生新个体是隐含方向性的,变异产生新个体则是随机的.本文交叉方式为单点交叉,交叉点范围为第二个基因到倒数第二个基因;变异为单点变异,变异点范围为所有基因.

图 2 交叉(a)和变异(b)示意图 加下划线的基因表示染色体变化的部分. Fig. 2 Schematic diagram of crossover (a) and mutation (b) Underlined genes are variable part of the chromosome.

(4) 种群大小、交叉和变异概率以及进化终止条件的设置

遗传算法寻优也不是绝对安全的,这取决于各参数(种群大小、交叉和变异概率以及进化终止的条件)的配合.为了节省计算时间,设置最大进化代数为60.一般设置交叉概率为较大常数,变异概率为很小的常数.在前几代种群内个体多样性丰富,交叉在创造新个体中起主要作用,随着进化的进行,种群多样性会逐渐减少,此时变异的作用就会凸显出来.本文效仿万永革和李鸿吉(1995)提出的变异概率随进化代数指数衰减的模型,增加变异在前几代的作用,相当于种群大小不变的情况下引入更多的个体,增加了搜索到最优个体的概率,同时也降低了搜索陷入局部最优的风险.变异概率随进化代数衰减情况见(8)式及图 3,式中参数取值需考虑到变异概率衰减的范围及速度,意在使进化过程中种群收敛的同时时刻保持种群多样性.这里使变异概率大约由0.1衰减至将近0.01,在最大进化代数的一半(30代),变异概率衰减至较小水平.经多次试算,第n代变异概率设定为

(8)

图 3 变异概率随进化代数指数衰减图 Fig. 3 Exponential reduction of mutation probability with generation

给定应力张量和震源断层的空间取向正演震源机制,基于这些震源机制,通过改变种群大小和交叉概率反演应力张量,综合考虑反演结果的拟合差、种群的收敛情况和计算速度,以确定这两个参数合适的值.选择操作不一定每次都能选出当代最优个体,交叉和变异也有可能破坏最优秀的个体.因此,在算法中需要对最优个体进行保护:保存当前最优个体,若下一代进化结果中出现更好的个体,将其替换为最优个体,把最后一代的最优个体作为结果输出.多次进行上述实验,最终确定种群大小为1500,交叉概率为0.7.

2 滑动方向和剪应力大小的偏差在反演应力张量中的表现

利用单个震源机制(震源机制的个数为1)多次实验,研究上一节叙及以滑动方向偏差(f11)和剪应力大小偏差(f21)为约束反演应力张量,哪种方案的表现更好,或者滑动方向偏差与剪应力大小偏差以什么样的组合(f31)会得到更理想的效果.这里任意给定一个应力状态(竖直方向挤压,东西向拉张,R=0.5)和断层取向(走向0°,倾角45°),计算得到断层面理论滑动方向(最大剪应力方向),在此方向基础上顺时针偏转0°~180°作为实际滑动方向,见图 4.根据(4)式得到不同约束方案f11f21f31(λ=0.5)和f41=sin(f1/2)的拟合差F1F2F3F4,观察随着最大剪应力方向与实际滑动方向角度偏差的增大(图 5横坐标),不同拟合差的变化情况.显然,在偏转角度为0°时所有的拟合差均为零,结果如图 5所示.

图 4 由理论滑动方向构造实际滑动方向示意图 α(0°~180°)是偏转角度. Fig. 4 Schematic diagram of actual sliding direction constructed by theoretical sliding direction α (0°~180°) is rotation angle.
图 5 不同拟合差随角度偏差的变化情况 Fig. 5 Variation of different misfit with angle deviation

图 5可以看出,随着剪应力方向与实际滑动方向角度偏差的增大,F1线性增加,F4正弦方式增加,F3介于F2F4之间,这些是可以预期的,并且四种拟合差都在单调增加,意味着,不论使用这四种约束方式中的哪个反演应力张量,均能得到正确结果.不同的是,这四种拟合差的增加速度不同,角度偏差较小时(这是我们关心的),增加速度的大小关系为:F1>F4>F3>F2,反映了四种约束方式对最大剪应力方向与实际滑动方向角度偏差的敏感程度,增加速度越快的对角度偏差越敏感.当角度偏差小于20°时,剪应力大小偏差(F2)对角度偏差极不敏感,如果用剪应力大小偏差作为约束反演应力张量,会增加遗传算法寻优陷入局部最优的风险.F1的斜率为1,因此,本文选择以滑动方向偏差为约束反演应力张量.

3 方法检验

本节方法检验有两个目的:(1)利用不同应力状态下人工合成的震源机制数据反演应力张量以检验不同应力状态下本文方法的有效性;(2)实际应力张量研究工作中,震源机制数据总是存在一定程度的噪声,这里在人工合成的震源机制数据中加入不同程度噪声来反演应力张量,并和Wan等(2016)四个参数(ϕδλR)搜索间隔分别为1°×1°×1°×0.1的网格搜索法反演的结果比较,以检验本文方法是否稳健.

理论上,应在人工合成震源机制数据的走向、倾角、滑动角中均添加噪声,但f1的表达式仅是滑动角的函数,走向和倾角的误差最终会传递到滑动角,因此,仅在人工合成震源机制的滑动角中添加不同程度的噪声,添加方式如下:

(9)

其中,rake为添加噪声的滑动角,rak为人工合成震源机制数据的滑动角,r为-1到1的随机数,noisy为添加的噪声大小.

本文选择和绝对安全的全局寻优算法——网格搜索法得到的结果相比较,一方面是因为网格搜索法可以给出一定精度下的最优解,与其结果比较,可以衡量本文方法所得结果是否正确;另一方面,Wan等(2016)在其文章中已证明网格搜索法优于Michael(1987)许忠淮(1985)的方法,若本文方法所得结果优于网格搜索法,也说明本文方法优于Michael(1987)许忠淮(1985)的方法.不同应力状态和噪声水平下本文方法与网格搜索法的反演结果见表 1,从表 1可以看到:不同应力状态,不同噪声水平的情况下,两种方法所得结果的拟合差差别较小,说明本文方法是正确的;本研究方法所得结果的拟合差均小于网格搜索法结果的拟合差,说明本文方法所得结果优于网格搜索法.该结果体现了本文方法的稳健性,同时也说明1.2节中遗传算法各参数设置是合理的.

表 1 不同应力状态和噪声水平下本文方法与网格搜索法的反演结果 Table 1 Inversion results by using the method in this study and the grid search algorithm under different stress states with different noisy

利用遗传算法反演应力张量的另一个优点是耗时短.本文选择种群大小为1500,最大进化代数为60,得出结果共需尝试9×104(1500×60)次;四个参数(ϕδλR)搜索间隔分别为1°×1°×1°×0.1的网格搜索法得出结果共需尝试1.1664×108(360×90×360×10)次.遗传算法计算的次数远远少于网格搜索法,是耗时短的原因.综上,本文方法在时间和精度方面均优于网格搜索法.

4 青藏高原东北缘应力张量的研究

将本文方法应用于青藏高原东北缘应力张量的研究,研究区域(103°E—106°E,34.5°N—37.5°N)是青藏高原、阿拉善块体、鄂尔多斯块体、华南块体的交界地区(图 6),该区域在这四个块体的共同作用下呈现一种什么样的应力状态,哪个块体对该区域的作用最大?这是本节亟待解决的问题.

图 6 震源机制分布图 图中沙滩球为震源机制的Schmidt投影,填充象限为压缩象限. Fig. 6 Distribution of focal mechanisms The beach balls in the figure are Schmidt projection of focal mechanisms, in which filled quadrant is compressed.
4.1 数据来源

研究区域内震源机制解数据来源于两部分,一是从赵知军和刘秀景(1990)利用P波初动得到宁夏及其邻区的震源机制解中挑选本文研究区域内的数据,并且将以下情况的数据剔除:(1)震源机制的PBT轴或两个节面明显不垂直(偏差大于10°);(2)震源机制两个节面和应力主轴明显不匹配(偏差大于10°);(3)P波初动所得震源机制解矛盾比大于20%.文献中给出的是震源机制解两个节面的走向、倾角和PBT轴的走向、仰角.双力偶震源机制解模型下,一个节面的法向和另一个节面的滑动方向一致或相反,因此将两个节面的法向矢量相加,若和矢量与P轴矢量一致,则一个节面的法向就是另外一个节面的滑动方向,否则一个节面法向的反方向是另外一个节面的滑动方向,据此得到两个节面的滑动角.二是来源于全球矩心矩张量目录(Global CMT).共得到震源机制解50个,时间跨度为1970—2015年,见表 2图 6.

表 2 研究区域内地震的震源机制解 Table 2 Focal mechanism solutions in the study region
4.2 应力张量反演结果

基于表 2数据,利用本文方法得到青藏高原东北缘的构造应力结果如图 7所示,σ1(张轴)走向152°,倾伏角50°;σ2(中间轴)走向350°,倾伏角39°;σ3(压轴)走向253°,倾伏角9°;R=0.64,根据第1节可知,中间应力轴表现为拉张作用,并且该拉张作用比σ1要弱很多.该结果说明,研究区域主要受控于青藏高原近东西向的挤压,从而导致在垂直于主压应力轴的平面上产生拉张作用,向华南块体方向的拉张作用强于向阿拉善块体的拉张作用,这与长期的动力作用导致研究区域南侧地形高于北侧是一致的!青藏高原对研究区域近东西向的挤压作用被认为来源于青藏高原受到印度板块北东向的挤压和青藏高原北部亚洲块体的阻挡,导致青藏高原物质不断隆升,以及侧向扩张(Shen et al., 2017李煜航,2018).另外,本文结果与其他学者对此区域的认识是一致的(许忠淮等,1987郭祥云等,2017),验证了本文方法和结果的正确性.

图 7 青藏高原东北缘构造应力张量反演结果 “+”表示应力主轴下半球投影的位置,黑色的弧线表示震源断层的下半球投影,玫红色小箭头表示断层的理论滑动方向,浅蓝色小箭头表示断层的实际滑动方向. Fig. 7 The inverted stress tensor in northeast edge of Tibetan Plateau "+" represents the position of stress axis in lower hemisphere projection. The black arc represents the fault in lower hemisphere projection. The small magenta arrow indicates the theoretical slip direction of fault. The small turquoise arrow indicates the actual slip direction of fault.
5 结论与讨论

本文给出了基于震源机制数据反演应力张量的遗传算法原理,分析了以断层面上滑动方向偏差为约束反演应力张量而忽略剪应力大小偏差对反演结果的影响,利用不同应力状态下人工合成的包含不同噪声大小的震源机制数据对本文方法进行了检验,最后研究了青藏高原东北缘的应力张量,主要结论如下:

(1) 随着断层面上滑动方向偏差的增大,剪应力大小的偏差也在单调增加,并且二者以不同的权重组合,其偏差也是单调增加的,意味着不论以哪种约束反演应力张量,均能得到正确的结果.然而,由于断层面上滑动方向偏差较小的时候,相对剪应力大小的偏差对其反映不感敏,容易使遗传算法搜索陷入局部最优,所以本文仅以断层面上滑动方向的偏差作为约束反演应力张量.

(2) 利用不同应力状态下人工合成的包含不同噪声大小的震源机制数据对本文方法进行检验,并与网格搜索法所得结果比较,结果显示,本文方法所得结果的拟合差均小于网格搜索法结果的拟合差,并且二者差别较小,体现了本文方法的稳健性.通过该实验,证明了本文方法在时间和精度上的优势.

(3) 对青藏高原东北缘应力张量的研究得到:青藏高原东北缘主要受控于青藏高原近东西向的挤压作用,从而导致在垂直于主压应力轴方向的平面上产生拉张作用,并且R值显示,向华南块体方向的拉张作用强于向阿拉善块体方向的拉张作用,这与其地形是一致的.

将人工智能技术应用于应力张量的研究,本文结果表明了该算法的高效性.遗传算法能否高效地搜索到全局最优解,取决于种群大小、适应度函数、交叉、变异的概率、进化终止条件之间的相互配合.最需要注意的是,任何时候都要保持种群个体的多样性,才能有效避免搜索陷入局部最优.应该指出的是,本文未给出一套完善的评估反演结果不确定范围的方案.作者曾尝试将进化最后一代个体中拟合差小于平均拟合差的个体作为不确定范围的评估,经过多次试验发现,对于同样的震源机制数据,由于遗传算法每次进化的轨迹不同,最后一代个体的分布可能不同,造成每次计算得到的不确定范围不尽相同.解决该问题是改进本文方法以后发展的一个方向.

本文研究开发的基于震源机制解反演应力张量的遗传算法MATLAB程序可向感兴趣的读者提供.

附录

应力张量的简化

遵照弹性力学(徐芝纶,2006)中的符号规定:拉张为正,挤压为负.在主轴坐标系(PBT)下,应力张量表达为

(A1)

σPBT做运算:,得到

(A2)

R∈[0, 1],称为应力形因子(Angelier, 1979, 1984Gephart et al., 1984Michael,1984Wan et al., 2016).将(A2)式球应力减去,得到构造应力张量为

(A3)

R描述了构造应力的状态:当R=0,σ2=σ3,对应于单轴拉张的应力状态;当R=0.5,中间应力轴(σ2)不表现构造作用;当R=1,σ1=σ2,对应于单轴挤压的应力状态(万永革等,2011).球应力不会影响断层上的剪应力,本文选择σ*作为主轴坐标系下简化的应力张量.

通过三次坐标旋转,可以将北东下坐标架旋转至以断层面的滑动方向为X轴,断层面上沿滑动方向顺时针旋转90°为Y轴,法向的反向为Z轴的任意方位.具体过程参考《地震学导论》(万永革,2016).将该过程的逆过程进行,就可以将任意坐标架旋转至北东下坐标架,从而得到北东下坐标系(NED)三个坐标轴的方向在任意坐标系(XYZ)下的表达(A4)式.

(A4)

其中,ϕδλ分别为断层的走向、倾角、滑动角.

主轴坐标系(PBT)下,三个坐标轴的方向为

(A5)

将式(A5)代入式(A4)得到北东下坐标系(NED)的三个坐标轴方向在主轴坐标系(PBT)下的表达:

(A6)

主轴坐标系(PBT)下的应力张量σ*在北东下坐标系(NED)表达为

(A7)

因此,应力张量简化为四个参数:ϕδλR,取值范围分别为:ϕ∈[0, 2π]、δ∈[0, π/2]、λ∈[0, 2π]、R∈[0, 1].

致谢  两位审稿专家对本文提供了有益的建议,特此致谢!
References
Angelier J. 1979. Determination of the mean principal directions of stresses for a given fault population. Tectonophysics, 56(3-4): T17-T26. DOI:10.1016/0040-1951(79)90081-7
Angelier J. 1984. Tectonic analysis of fault slip data sets. Journal of Geophysical Research: Solid Earth, 89(B7): 5835-5848. DOI:10.1029/JB089iB07p05835
Christova C V. 2015. Spatial distribution of the contemporary stress field in the Kurile Wadati-Benioff zone by inversion of earthquake focal mechanisms. Journal of Geodynamics, 83: 1-17. DOI:10.1016/j.jog.2014.11.001
Gephart J W, Forsyth D W. 1984. An improved method for determining the regional stress tensor using earthquake focal mechanism data: Application to the San Fernando earthquake sequence. Journal of Geophysical Research: Solid Earth, 89(B11): 9305-9320. DOI:10.1029/JB089iB11p09305
Guo X Y, Jiang C S, Wang X S, et al. 2017. Characteristics of small to moderate focal mechanism solutions stress field of the circum-Ordos block. Journal of Geodesy and Geodynamics (in Chinese), 37(7): 675-685.
Kobayashi R, Nakanishi I. 1994. Application of genetic algorithms to focal mechanism determination. Geophysical Research Letters, 21(8): 729-732. DOI:10.1029/94GL00593
Li Y H. 2018. Study on the lateral motion of Northeastern Tibetan Plateau. Recent Developments in World Seismology (in Chinese), (5): 45-47.
Michael A J. 1984. Determination of stress from slip data: faults and folds. Journal of Geophysical Research: Solid Earth, 89(B13): 11517-11526. DOI:10.1029/JB089iB13p11517
Michael A J. 1987. Use of focal mechanisms to determine stress: A control study. Journal of Geophysical Research: Solid Earth, 92(B1): 357-368. DOI:10.1029/JB092iB01p00357
Shen X Z, Liu M, Gao Y, et al. 2017. Lithospheric structure across the northeastern margin of the Tibetan Plateau: Implications for the plateau′s lateral growth. Earth and Planetary Science Letters, 459: 80-92. DOI:10.1016/j.epsl.2016.11.027
Sheng S Z, Wan Y G, Xu Z G, et al. 2013. Mean stress field inferred from the total seismic moment released by earthquakes. Seismology and Geology (in Chinese), 35(1): 92-100.
Shi Y L, Jin W. 1995. Genetic algorithms inversion of lithospheric structure from surface wave dispersion. Acta Geophysica Sinica (in Chinese), 38(2): 189-198.
Thakur P, Srivastava D C, Gupta P K. 2017. The genetic algorithm: A robust method for stress inversion. Journal of Structural Geology, 94: 227-239. DOI:10.1016/j.jsg.2016.11.015
Wan Y G, Li H J. 1995. The preliminary study on the seismic hypocenter location using genetic algorithms. Seismological and Geomangnetic Observation and Research (in Chinese), 16(6): 1-7.
Wan Y G. 2010. Contemporary tectonic stress field in China. Earthquake Science, 23(4): 377-386. DOI:10.1007/s11589-010-0735-5
Wan Y G, Sheng S Z, Xu Y R, et al. 2011. Effect of stress ratio and friction coefficient on composite P wave radiation patterns. Chinese Journal of Geophysics (in Chinese), 54(4): 994-1001. DOI:10.3969/j.issn.0001-5733.2011.04.014
Wan Y G, Sheng S Z, Huang J C, et al. 2016. The grid search algorithm of tectonic stress tensor based on focal mechanism data and its application in the boundary zone of China, Vietnam and Laos. Journal of Earth Science, 27(5): 777-785. DOI:10.1007/s12583-015-0649-1
Wan Y G. 2016. Introduction to Seismology (in Chinese). Beijing: Science Press: 393-394.
Wu Y M, Zhao L, Chang C H, et al. 2008. Focal-mechanism determination in Taiwan by genetic algorithm. Bulletin of the Seismological Society of America, 98(2): 651-661. DOI:10.1785/0120070115
Xu Z H. 1985. Mean stress field in Tangshan aftershock area obtained from focal mechanism data by fitting slip direction. Acta Seismologica Sinica (in Chinese), 7(4): 349-362.
Xu Z H, Wang S Y, Huang Y R, et al. 1987. Directions of mean stress axes in southwestern China deduced from microearthquake data. Acta Geophysica Sinica (in Chinese), 30(5): 476-486.
Xu Z H. 2001. A present-day tectonic stress map for eastern Asia region. Acta Seismologica Sinica (in Chinese), 23(5): 492-501.
Xu Z L. 2006. Elasticity (in Chinese). Beijing: Higher Education Press: 4-5.
Zhao Z J, Liu J X. 1990. Seismic activity and local tectonic stress field in Ningxia and nearby regions. Seismology and Geology (in Chinese), 12(1): 31-46.
Zhong J M, Cheng W Z. 2006. Determination of directions of the mean stress field in Sichuan-Yunnan region from a number of focal mechanism solutions. Acta Seismologica Sinica (in Chinese), 28(4): 337-346.
郭祥云, 蒋长胜, 王晓山, 等. 2017. 鄂尔多斯块体周缘中小地震震源机制及应力场特征. 大地测量与地球动力学, 37(7): 675-685.
李煜航. 2018. 青藏高原东北横向扩展运动研究. 国际地震动态, (5): 45-47. DOI:10.3969/j.issn.0253-4975.2018.05.010
盛书中, 万永革, 徐志国, 等. 2013. 由地震释放的地震矩叠加推导平均应力场. 地震地质, 35(1): 92-100. DOI:10.3969/j.issn.0253-4967.2013.01.008
石耀霖, 金文. 1995. 面波频散反演地球内部构造的遗传算法. 地球物理学报, 38(2): 189-198. DOI:10.3321/j.issn:0001-5733.1995.02.007
万永革, 李鸿吉. 1995. 遗传算法在确定震源位置中的应用. 地震地磁观测与研究, 16(6): 1-7.
万永革, 盛书中, 许雅儒, 等. 2011. 不同应力状态和摩擦系数对综合P波辐射花样影响的模拟研究. 地球物理学报, 54(4): 994-1001. DOI:10.3969/j.issn.0001-5733.2011.04.014
万永革. 2016. 地震学导论. 北京: 科学出版社: 393-374.
许忠淮. 1985. 用滑动方向拟合法反演唐山余震区的平均应力场. 地震学报, 7(4): 349-362.
许忠淮, 汪素云, 黄雨蕊, 等. 1987. 由多个小震推断的青甘和川滇地区地壳应力场的方向特征. 地球物理学报, 30(5): 476-486. DOI:10.3321/j.issn:0001-5733.1987.05.005
许忠淮. 2001. 东亚地区现今构造应力图的编制. 地震学报, 23(5): 492-501. DOI:10.3321/j.issn:0253-3782.2001.05.005
徐芝纶. 2006. 弹性力学. 北京: 高等教育出版社: 4-5.
赵知军, 刘秀景. 1990. 宁夏及其邻区地震活动带与小区域构造应力场. 地震地质, 12(1): 31-46.
钟继茂, 程万正. 2006. 由多个地震震源机制解求川滇地区平均应力场方向. 地震学报, 28(4): 337-346. DOI:10.3321/j.issn:0253-3782.2006.04.001