| 辐射花样的计算与震源球的绘制（转） | 您所在的位置：网站首页 › 沙滩如何画 › 辐射花样的计算与震源球的绘制（转） |
辐射花样的计算与震源球的绘制（转）
|
辐射花样的计算与震源球的绘制（转）
P 波辐射花样计算公式坐标系的选取震源机制解辐射花样计算代码震源球的绘制投影方式的选取绘图脚本图像格式转换
P 波辐射花样计算公式
Quantitative Seismology (Aki and Richards,1980) 中式（4.29）给出了零迹地震矩 M 所产生的 n 分量位移公式： 式（4.91）中给出了远场 P 波位移的矢量形式，看上去更加直观一些： 与时间相关的部分可以分离成震源时间函数），离源矢量是与方位角和离源角有关的矢量，所以求辐射花样的本质就是矢量和张量的乘法。如何选定坐标系是一个关键问题。 按照 Aki&Richards(1980) 中图 4.20 的方式定义坐标系，如下图，定义 X 轴为北向，Y 轴为东向，Z 轴为垂直向下，即 NED 坐标系。 震源机制解一般有两种表达方式，一种是矩张量形式，另一种是断层参数形式。 矩张量形式是震源机制的通用表示方式，需要六个分量。对于地震震源而言，多限制矩张量为零迹张量，即去除爆炸源的成分，只保留 double couple 和 CLVD 部分。断层参数形式需要三个分量 (strike,dip,rake)，只能表示 double couple 位错源。Global CMT 给出了零迹矩张量解和断层参数解。 若使用 GCMT 给出的断层参数 (strike,dip,rake) 解，则可根据 Aki&Richards(1980) P117 Box4.4 中式 1 将其转换成 NED 坐标系下的矩张量。 若使用 GCMT 给出的矩张量解，由于 GCMT 给出的是 (Mrr, Mtt, Mff, Mrt, Mrf, Mtf) 解，即 USE 坐标系下的矩张量，需要转换成 NED 坐标系的矩张量，方可使用。
不同的文献给出的坐标系可能不同，比如这里提到的 NED 坐标系和 USE 坐标系。即便相同的坐标系所使用的符号也可能不同，比如 GCMT 的 （r,t,f）坐标系和 Aki&Richards(1980) 中给出的(r,) 坐标系其实都是 USE 坐标系。 辐射花样计算代码获得矩张量以及离源矢量的表达式之后，即可求出震源球上任一点的辐射振幅。代码如下： #include #include #include #define PI 3.14159265358979323846 #define DEG2RAD PI/180.0 int main (int argc, char *argv[]) { int i, j; float m[3][3]; if (argc != 7 && argc != 4) { fprintf(stderr,"Usage: %s mrr mtt mff mrt mrf mtf\n", argv[0]); fprintf(stderr," Or: %s strike dip rake\n", argv[0]); exit(1); } if (argc == 7) { // moment tensor sscanf(argv[1], "%f", &m[2][2]); // mrr -> mzz sscanf(argv[2], "%f", &m[0][0]); // mtt -> mxx sscanf(argv[3], "%f", &m[1][1]); // mff -> myy sscanf(argv[4], "%f", &m[2][0]); // mrt -> mzx m[0][2] = m[2][0]; sscanf(argv[5], "%f", &m[2][1]); // mrf -> -Mzy m[2][1] = -m[2][1]; m[1][2] = m[2][1]; sscanf(argv[6], "%f", &m[0][1]); // mtf -> -Mxy m[0][1] = -m[0][1]; m[1][0] = m[0][1]; } else if (argc == 4) { // strike, dip, rake float strike, dip, rake; sscanf(argv[1], "%f", &strike); sscanf(argv[2], "%f", &dip ); sscanf(argv[3], "%f", &rake ); strike *= DEG2RAD; rake *= DEG2RAD; dip *= DEG2RAD; m[0][0] = - sin(dip)*cos(rake)*sin(2*strike) - sin(2*dip)*sin(rake)*sin(strike)*sin(strike); m[0][1] = sin(dip)*cos(rake)*cos(2*strike) + 0.5*sin(2*dip)*sin(rake)*sin(2*strike); m[0][2] = -cos(dip)*cos(rake)*cos(strike) - cos(2*dip)*sin(rake)*sin(strike); m[1][1] = sin(dip)*cos(rake)*sin(2*strike) - sin(2*dip)*sin(rake)*cos(strike)*cos(strike); m[1][2] = -cos(dip)*cos(rake)*sin(strike) + cos(2*dip)*sin(rake)*cos(strike); m[2][2] = sin(2*dip)*sin(rake); m[1][0] = m[0][1]; m[2][0] = m[0][2]; m[2][1] = m[1][2]; } fprintf(stdout," / %6.3f %6.3f %6.3f \\ \n", m[0][0], m[0][1], m[0][2]); fprintf(stdout,"M = | %6.3f %6.3f %6.3f | \n", m[1][0], m[1][1], m[1][2]); fprintf(stdout," \\ %6.3f %6.3f %6.3f / \n", m[2][0], m[2][1], m[2][2]); FILE *fop; fop = fopen("pattern.dat", "wb"); double az, theta; float p[3]; // 离源矢量 for (i=0; i amp += p[k]*m[k][l]*p[l]; } fwrite(, sizeof(float), 1, fop); } fclose(fop); return 0; }此代码可以正确处理断层参数和矩张量两种形式的震源机制解，二者均可被正确转换为 NED 坐标系下的矩张量解。对 360 度的方位角以及 90 度的离源角进行遍历，计算每一点的振幅值，并保存到 pattern.dat 中待用。 关于离源角，需要注意两点： 离源角的取值范围为 [0,90]，即只计算震源球的下半球，这是因为多数情况下绘制震源球都使用下半球投影（上半球辐射的能量无法传播到大震中距处）。 离源角与纬度的对应关系为：纬度 = 离源角 - 90。 震源球的绘制 投影方式的选取目前已经拥有了震源球的下半球上任意一点的振幅（未归一化），还需要选择合适的投影方式将数据投影到 “赤道” 面上。 绘制震源球有两种投影方式，分别是 Schmidt 投影和 Wulff 投影。前者是等面积投影，后者是等角度投影。在 GMT 中分别对应 JA 和 JS 。这里以 Wulff 投影为例，想要使用 Schmidt 投影只需要把 JS 改成 JA 即可。 绘图脚本 #!/bin/bash R=0/360/-90/0 J=S0/-90/15c B=a30f10 name=pattern PS=${name}.ps gmt set MAP_FRAME_TYPE=plain gmt set FORMAT_GEO_MAP=+D gmt xyz2grd ${name}.dat -G${name}.nc -I6m/6m -R$R -ZLBxf gmt grd2cpt ${name}.nc -Cpolar -E100 > ${name}.cpt gmt psxy -R$R -J$J -T -K -P > $PS gmt grdimage ${name}.nc -R$R -J$J -C${name}.cpt -B$B -BN -K -O >> $PS gmt grdcontour ${name}.nc -R$R -J$J -L-0.001/0.001 -C1 -K -O -W2p >> $PS gmt psxy -R$R -J$J -T -O >> $PS rm gmt.conf gmt.history ${name}.cpt ${name}.nc绘图脚本的一些说明： 设置 FORMAT_GEO_MAP 使得方位角范围是 0 到 360，而不是 -180 到 180。其中 0 度指向正北方向。这里 R 的横向范围是 0 到 360，实际上 360 度处与 0 度处是同一个经度，所以网格中没有计算 360 度处的振幅。同时在 -Z 选项中使用了 x 以表明 X 轴的周期性。在振幅为 0 处绘制了等值线。 图像格式转换利用 psconvert 命令可以将 PS 文件转换为其它格式的图像，最好选择透明的 PNG 格式 gmt psconvert -A -TG beachball.ps。转至Seisman |
其中等式右边共计 5 项，第一项为近场项，第二、三项分别是 P、S 波的中间场项，第四、五项分别为 P、S 波的远场项。一般研究辐射花样大多关注于 P 波远场辐射花样，即第四项。
其中，除去震源时间函数以及绝对振幅，只留下辐射花样相关的因子：
R
a
d
=
r
.
M
.
r
Rad = r.M.r
Rad=r.M.r 其中，r 为离源矢量，是离源角和方位角的函数，表征了地震射线从震源的出射方向，M为矩张量。
可以得到，此坐标下，离源矢量 r 的具体形式:【本文地址】