基于Cesium的震源机制解三维可视化方法

何俊强，杨 钦,2

(1.北京网格天地软件技术股份有限公司 研发部，北京 100088；2.北京航空航天大学 计算机学院，北京 100191)

1 引 言

近年来，基于Web的三维可视化技术有了较大的发展，其中Cesium作为一款开源的三维GIS(Geographic Information Systems,GIS)引擎深受广大学者的欢迎，基于Cesium的各类数据可视化方法研究也是越来越广泛[1-4]，如王晓艳等[5]实现了将3ds Max格式的三维模型转成glTF和3D Tiles两种格式后，在Cesium进行三维可视化的方法；郭晓非等[6]使用Cesium实现了海表面温度场和海面风场的可视化，使得海洋环境的模拟更加形象直观。另一方面，国家在大力发展数字地球的建设，积极推动地上地下一体化的可视化技术研究。2022年，全国两会提出了实施“透明国土”和建设“全息数字地球基础框架”的提案，预示着以三维地质模型为基础的地震、地质、水文、物探、化探等信息的存储、共享、可视化与分析应用将是未来的一大发展与研究方向[7]。

震源机制解是地震领域的重要参数之一，在研究地震断层类型、断层稳定性分析、发震机制、应力分布和区域地质构造等工作中起着重要作用[8-12]。震源机制解的可视化结果是沙滩球，其本质是震源球的下半球在赤平面的投影，所以地震领域中的沙滩球是二维的。目前震源机制解常用的绘制工具有GMT(Generic Mapping Tools,GMT)等[13]，其绘制出来的结果是二维的沙滩球图片，将其放到三维信息化平台中的显示效果一般，且与使用三维地质模型等数据构建出来的场景也不相协调，考虑到三维信息化建设的发展方向以及人们对绘制模拟更形象直观的需求，本文基于Cesium对震源机制解的三维可视化方法进行了研究。

2 方法介绍

2.1 三维沙滩球的引入

以震源为球心，做一个足够小的球面，小到球内射线弯曲可以忽略不计，即假定此小球内的介质是均匀的，这个小球面就称为震源球[14]。震源向各个方向发出地震射线，震源球表面每个点都对应着一根射线的振幅，由于地震时震源周围的物理环境及应力分布并不均匀，所以振幅有正有负，将震源球上正振幅的点用黑色表示，负振幅的点用白色表示，只要点足够多，就能在震源球上形成一个均匀分布的黑白图案，称为辐射花样，将这个带有辐射花样的震源球的下半球投影到赤平面就得到了沙滩球。从这个原理可以认为，震源机制解的三维可视化结果其实就是带有辐射花样的震源球，本文称之为三维沙滩球。

2.2 辐射花样的计算

三维沙滩球的辐射花样是地震射线振幅值在震源球面的正负分布[15,16]，因此需要对振幅进行计算。Aki等[17]给出了远场P波辐射花样振幅的计算公式

A=γ·M·γ

(1)

式(1)中,A为辐射花样振幅，单位为m;M为矩张量，是NED(North East Down)坐标系(X轴朝北，Y轴朝东，Z轴垂直向下)下的3阶对称矩阵；γ为离源矢量，是离源角和方位角的函数，离源角是指从震源发射出去的射线与沿着震源向下的垂线之间的夹角。

矩张量M可以通过震源机制解的参数计算得到。震源机制解的参数有两种类型，分别是矩张量和断层参数。震源机制解的矩张量是USE(Up South East)坐标系(X轴垂直向上，Y轴朝南，Z轴朝东)下的3阶对称矩阵，包含Mrr、Mtt、Mpp、Mrt、Mrp、Mtp6个参数，那6个参数里的R、T、P对应的是USE坐标系下的三个轴，为了跟NED坐标系做区分，USE下的轴叫R、T、P，NED坐标系的就叫X、Y、Z，Mxx和Mrr都是对应轴的分量。可通过坐标转换公式得到NED坐标系下的矩张量，其转换公式为

(2)

式(2)中,左侧为NED坐标系的矩张量，右侧为USE坐标系的矩张量。

震源机制解的断层参数包括走向、倾角、滑动角3个参数，为角度制，在转成弧度制后，可通过公式转换成NED坐标系下的矩张量，其转换公式为

(3)

式(3)中,δ为倾角，λ为滑动角，φs为走向，单位均为rad；M0为常量。

离源矢量γ可通过离源角和方位角计算得到，其计算公式为

(4)

式(4)中,θ为离源角，α为方位角，单位均为rad。

2.3 三维沙滩球的绘制

三维沙滩球的绘制就是绘制带有辐射花样图案的震源球。考虑到直接在球体上绘制辐射花样图案比较困难，本文选择将辐射花样图案当做纹理贴在球体上的方式来实现，由此可将绘制流程分成两个步骤，一是制作辐射花样纹理贴图，二是绘制带有辐射花样纹理贴图的震源球。

球体的纹理其实就是球面展开后的样子，与将地球通过墨卡托投影后得到的地图类似，因此可遍历0°～360°的方位角和0°～180°的离源角，利用式(1)分别计算震源球面每一个点的振幅，然后将这些点绘制在以方位角为横轴，离源角为纵轴的图上，并分别对振幅大于0和小于0的点设置不同的颜色，即可绘制出辐射花样纹理图。关于颜色设置，目前常见的是正振幅为黑色，负振幅为白色，本文根据震源机制解的滑动角范围对断层类型进行了划分(表1)[18]，并将正断层的正振幅设为绿色，逆断层的正振幅设为红色，走滑断层的正振幅设为蓝色，负振幅均设为白色，这样通过颜色即可判断断层类型。

Cesium是开源的三维GIS引擎，通过其官方的接口可以绘制带纹理的点线面体等基础几何体，三维球体可通过向entities里添加ellipsoid类型对象进行构建，构建时需要设置球心的三维坐标、球体的半径以及纹理的路径3个参数。震源球的球心坐标就是震源的发震位置，将震源的经纬度和高程值设置给ellipsoid对象的position参数即可。震源球的半径目前没有统一的规定，比较合理的方案是根据地震震级的大小做区分，震级越大,半径越大。ellipsoid的半径参数为radii，是一个三维数组，由于是球体，三个值均设置成一致。最后将辐射花样纹理图的路径设置给ellipsoid对象的material参数，就可在震源球上绘制出辐射花样纹理。

纯色的球绕着过球心的轴无论如何转动，从各个方向看过去都一样，但是带纹理的球则不同。Cesium绘制的球体默认坐标系为ENU(X轴朝东，Y轴朝北，Z轴垂直向上)，纹理图片的左右边界(0°和360°)会在球体的正东经线处重合，而辐射花样的方位角是从正北计算的，需要将球体沿着垂直向上的轴逆时针旋转90°。ellipsoid对象的三轴转动参数为orientation，其中绕Z轴旋转的参数为heading，将其设置为-90°。这样与实际相符的带有辐射花样的震源球即三维沙滩球就绘制好了。

2.4 实现流程

综合上述分析，总结出震源机制解的三维可视化流程如图1所示，其具体步骤为：

1)输入断层参数类型或矩张量类型的震源机制解参数；

2)将参数转换成NED坐标系下的矩张量。断层参数类型使用式(3)转换，矩张量类型使用式(2)转换；

3)使用式4)计算震源球面每一个点的离源矢量；

4)使用式(1)计算震源球面每一个点的辐射振幅；

5)将球面所有点绘制在以方位角为横轴，离源角为纵轴的图上，根据断层滑动角判断断层类型，并对振幅大于0和小于0的点填充相应的颜色，得到辐射花样纹理；

6)利用Cesium绘制以震源点坐标为中心的球体，并贴上辐射花样纹理；

7)将震源球体朝北旋转90度得到最终的三维沙滩球。

3 方法验证

3.1 测试数据介绍

为了验证使用上述方法绘制出来的三维沙滩球是否正确，在欧洲地中海地震中心网站(European-Mediterranean Seismological Centre)提取了3条震源机制解数据进行测试[19]，数据如表2所示，分别为2021年8月27日发生在印度尼西亚南部海域的地震、2021年11月17日发生在中国黄海海域的地震、2022年1月9日发生在印度尼西亚的地震，数据均包含断层参数和矩张量两种形式。通过滑动角数据以及表1的对应关系，可以判断出这三个震源机制解对应的断层类型分别为正断层、逆断层和走滑断层，确保了数据的多样性。

3.2 测试数据绘制

根据图1的流程分别对表2的三条震源机制解参数进行计算，得到辐射花样纹理如图2、图3、图4所示。其中，图2(a)、图3(a)、图4(a)为断层参数对应的纹理，图2(b)、图3(b)、图4(b)为矩张量对应的纹理，横轴为0°～360°的方位角，对应球体的经度;纵轴为0°～180°的离源角，减去90°就是纬度，根据断层类型对正振幅区域填充的颜色分别为绿色、红色和蓝色。

图2 正断层辐射花样纹理Fig.2 Radiation pattern texture of normal fault parameters

图3 逆断层辐射花样纹理Fig.3 Radiation pattern texture of reverse fault parameters

图4 走滑断层辐射花样纹理Fig.4 Radiation pattern texture of strike slip fault parameters

使用表2中的经纬度和深度在Cesium框架下绘制震源球，然后分别使用上述纹理图进行贴合并向北旋转90°得到三维的沙滩球。图5给出了北半球高空视角下的三维沙滩球显示效果，考虑到实际的震源在地下，不方便查看，所以关掉了Cesium的深度检测，这样可以看到地球内部的物体，可以看出绘制出来的图案与二维的沙滩球相似，但还需要进行对比验证。

图6 断层参数沙滩球对比Fig.6 Comparison of beach balls of fault parameters

3.3 测试数据验证

由于沙滩球是震源球下半球在赤平面的投影，可以通过对比使用本文的方法绘制出来的三维沙滩球的下半球视图与使用主流的GMT绘制出来的沙滩球是否一致的方法来判断。考虑到三维沙滩球的下半球在地球内部，可以将三维沙滩球绕着其局部坐标系的Y轴和Z轴分别顺时针旋转180°使其上下反转，同时将观察视角设置为沙滩球位置的高空俯视视角，这样观察到的就是三维沙滩球的下半球。

图7 矩张量沙滩球对比Fig.7 Comparison of beach balls of moment tensor

图6给出了断层参数对应的三维沙滩球与GMT绘制的二维沙滩球的对比，为了展示三维沙滩球在不同视角下的形态，对比时没有隐藏其他两个沙滩球，为了指明图中当前视角观察的对象，添加了红色箭头来指示。通过对比可以发现两者的花纹基本一致，只是花纹的边界略有差别，其中正断层差异明显一些，相比于GMT绘制的沙滩球，本方法看到的下半球白色区域更大一些。图7给出了矩张量对应的沙滩球的对比，与图5一样，两者在形态上基本一致，但花纹的边界差异比矩张量更明显一些，尤其是逆断层，GMT绘制的沙滩球在右下角内凹很明显。考虑到地震中的沙滩球是通过乌尔夫网投影[20]得到的，是一种等角度投影，在相同角度的情况下，越靠近两侧的区域投影后面积越大，最后产生网心部分过于集中，边缘部分则较稀疏的情况，而本文使用的俯视图法不存在上述问题，所以两者对比存在的差异是合理的，使用本方法得到的三维沙滩球在一定程度上来说是正确的。

图8 汶川龙门山三维地质场景Fig.8 3D beach balls of Wenchuan earthquake

4 方法应用

三维沙滩球可作为断层面建模的参考依据，其颜色分界面可以看做断层面的走向，在实际断层面建模中可以起到约束作用，同时结合地质模型等数据可以构建三维地质场景。图8给出了由汶川龙门山的三维P波速度模型、龙门山的主断层面[21]、断层轨迹线以及对应的未上下反转的三维沙滩球[22]构建的三维地质场景，由于实际的数据在地下，不便于查看，本文将其整体抬升到了地表之上；另外由于沙滩球和断层面在地质模型内部，本文也对模型部分区域进行了切割，以便于观察。可以看出速度场和主断层面的展布与沙滩球的分布大体上一致，断层面的走向与沙滩球的分界面吻合得一般，还有待完善。

5 结 论

本文结合地震数据三维信息化发展的背景，从震源机制解沙滩球的绘制原理出发，基于目前常用的三维可视化引擎Cesium，研究了震源机制解三维可视化的方法并进行了分析讨论，最终得到如下结论：

1)震源机制解的三维可视化结果是带有辐射花样图案的震源球。

2)通过基于辐射花样计算公式得到震源球面的辐射花样纹理，并使用Cesium绘制带纹理的球体的方式来实现震源机制解的三维可视化的方法是正确可行的。

3)使用本方法绘制出来的三维沙滩球形象直观，可通过将其进行反转用于常规的震源机制解数据分析，也可用作断层面建模以及三维信息化场景的构建等，使得相关领域的研究以及向大众普及地震相关的知识更加地方便、有效。