单纯形定位方法在新疆数字地震台网的测定精度分析

张志斌; 金花; 王晓飞

doi:10.11899/zzfy20190119

单纯形定位方法在新疆数字地震台网的测定精度分析

doi: 10.11899/zzfy20190119

新疆维吾尔自治区地震局, 乌鲁木齐 830011

基金项目:

测震台网青年骨干专项 CEA-JC/QNCZ-18313

新疆地震科学基金 201707

详细信息

作者简介:
张志斌, 男, 生于1988年。工程师。主要从事地震定位和震源机制方面研究。E-mail:306920839@qq.com

计量
- 文章访问数: 111
- HTML全文浏览量: 38
- PDF下载量: 9
- 被引次数: 0
出版历程
- 收稿日期: 2018-06-28
- 刊出日期: 2019-03-01

Analysis of Precision of the Simplex Location Method in the Xinjiang Digital Seismic Network

Earthquake Agency of Xinjiang Uygur Autonomous Region, Urumqi 830011, china

摘要

摘要: 地震定位是地震学中最基本的问题之一，现阶段大部分采用基于走时的方法。单纯形定位方法作为一种直接搜索类的方法，在新疆测震台网日常地震目录产出中发挥了重要作用。本文采用数值模拟的方法，利用新疆测震台网现有布局，结合“3400”走时表，分析单纯形定位方法在新疆测震台网的测定精度。研究表明，对于新疆测震台网网内浅源地震，单纯形定位方法能够得出较为精确的震中位置，而得出的网外地震震中位置存在一定偏差；初至折射波对震源深度有一定的控制能力，速度模型的准确性对震源深度和发震时刻的测定影响较大。
- 新疆地区 /
- 地震定位 /
- 震源深度 /
- 单纯形定位法
Abstract: Seismic location is one of the most basic issues in seismology. Currently, most of earthquake location used a time-based approach. However, the simplex location method, as a direct search method, plays an important role in the daily earthquake catalog of the Xinjiang Digital Seismic Network. This paper applies numerical simulation methods, using the Xinjiang Digital Seismic Network and combing the "3400" travel timetable, and analyzes the measurement accuracy of the simplex location method in the Xinjiang Digital Seismic Network. The results show that for the shallow earthquakes in the Xinjiang Digital Seismic Network, the simplex location method can get a high-accurate epicenter location, However, there exists a certain deviation in outside the seismic network. The first arrival refracted wave has a certain ability to control the depth of the epicenter. The accuracy of the velocity model has a great influence on the depth of the epicenter and the determination of the time of earthquake occurrence.
- Xinjiang region /
- Seismic location /
- Focal depth /
- Simplex location method

HTML全文

引言

近年来，随着世界经济的发展，中国近海结构的研究得到了快速发展，如跨海桥梁、海上风电、人工岛和石油平台等。然而，中国福建、广东沿海和中国台湾等地区的近海结构面临着地震作用的威胁。在地震作用下，桥梁与周围水体的相互作用会对桥墩产生动水压力。地震作用下水体与结构的动力相互作用会对结构产生附加动水压力，其不仅会改变结构的动力特性，还会对水中结构的动力反应产生较大影响(Liaw等，1974；Han等，1996；黄信等，2011a；江辉等，2014)。因此，充分认识和研究地震激励下近海结构的动水压力，对于近海结构的抗震设计具有重要意义。

在实际近海工程中存在各种截面形式的墩柱结构，如圆形和椭圆形等。目前，国内外学者对圆柱结构地震动水压力的研究已取得大量的研究成果。Liaw等(1974)基于辐射波浪理论推导了圆柱结构的动水压力公式，结果表明，对于细长结构可以忽略水体压缩性的影响；表面波仅在低频时对动水压力的影响比较明显。忽略水体压缩性和表面波时，该动水压力可视为部分水体质量与结构加速度的乘积，这部分水体称为“附加质量”。Williams(1986)采用边界积分方法研究了地震作用下水中圆柱结构的动力反应。Tanaka等(1988)分析了水平地震作用下弹性圆柱体动水压力的附加质量系数和阻尼系数。Han等(1996)提出了计算水中圆柱结构自振频率的简化公式。赖伟等(2004)提出了一种圆形桥墩上地震动水压力的半解析半数值方法。黄信等(2011b, 2012)讨论了水体压缩性、表面波和水底吸收边界对圆形桥墩地震动水压力的影响。杜修力等(2012)、Du等(2014)提出了可压缩水体条件下圆柱结构地震动水压力的时域算法和时域简化公式。Goyal等(1989)、Li等(2013)和Jiang等(2017)提出了圆柱形桥墩上地震动水压力的附加质量简化计算公式。另外，Liao(1985)研究了水中多个圆柱结构的动力相互作用；Sun等(1991)研究了轴对称结构的地震动水压力；Avilés等(2001)讨论了海底吸收条件对轴对称结构地震动水压力的影响；Wei等(2015)提出了水中轴对称结构的地震设计和分析的简化方法。

上述分析表明，地震作用下圆柱结构与水体的动力相互作用问题已取得了可供实际应用的研究成果，但对椭圆柱体结构动水压力方面仍鲜有研究。本文基于椭圆坐标系，根据线性辐射波浪理论，利用分离变量法推导了椭圆柱体结构动水压力的解析解，并采用有限元方法建立了地震作用下水体与椭圆柱体结构动力相互作用的简化分析方法。

1. 椭圆形柱体动水压力求解

地震作用下椭圆柱体结构与水体动力相互作用的分析模型如图 1所示，a和b分别表示椭圆外径的半长轴和半短轴，a₁和b₁分别表示椭圆内径的半长轴和半短轴，h为水深，H为柱体高度。直角坐标系下，z轴沿柱体轴线向上，坐标原点位于柱体底部；地基为刚性，地面运动加速度为$\ddot{u}$_g。水体假定为无旋、无粘、可压缩的小扰动流体，并忽略表面重力波的影响。

图 1 椭圆柱体与水体相互作用分析模型

Figure 1. Analytical model of the interactiin of elliptic cylinder and water

下载: 全尺寸图片幻灯片

1.1 控制方程和边界条件

由于柱体截面呈椭圆形，故在椭圆坐标系下求解椭圆柱体的地震动水压力。椭圆坐标系如图 2所示。直角坐标系与椭圆坐标系直接的变化关系为：

$$ x=\mu \cos h\xi \cos \eta $$

(1)

$$ y=\mu \text{sin}h\xi \text{sin}\eta $$

(2)

图 2 椭圆坐标系

Figure 2. Elliptical cylindrical coordinate

下载: 全尺寸图片幻灯片

式中，$\mu =\sqrt{{{a}^{2}}-{{b}^{2}}}$，$\xi $和η分别为椭圆坐标系的径向和环向坐标，取值范围分别为$0\le \xi <\infty $和$ 0\le \eta < 2\text{ }\pi\text{ }$。其中，椭圆柱体表面的环向坐标为：

$$ {{\xi }_{0}}=\text{ta}{{\text{n}}^{-1}}(b/a) $$

(3)

在椭圆柱坐标系下，以动水压力p表示的流体控制方程为(Bhatta，2005)：

$$ \frac{2}{{{\mu }^{2}}(\text{cos}h2\xi-\text{cos}2\eta)}\left(\frac{{{\partial }^{2}}p}{\partial {{\xi }^{2}}}+\frac{{{\partial }^{2}}p}{\partial {{\eta }^{2}}} \right)+\frac{{{\partial }^{2}}p}{\partial {{z}^{2}}}=0 $$

(4)

水体底部、静水表面和无穷远边界条件为：

$$ \frac{\partial p}{\partial z}\left| _{z=0} \right.=0 $$

(5)

$$ p\left| _{z=h} \right.=0 $$

(6)

$$ p\left| _{\xi \to \infty } \right.=0 $$

(7)

地震动沿长轴方向时，η=0和η=0.5π面的对称边界条件为：

$$ \frac{\partial p}{\partial \eta }\left| _{\eta =0} \right.=0 $$

(8a)

$$ p\left| _{\eta =0.5\pi } \right.=0 $$

(8b)

地震动沿短轴方向时，η=0和η=0.5π面的对称边界条件为：

$$ \frac{\partial p}{\partial \eta }\left| _{\eta =0.5\pi } \right.=0 $$

(9a)

$$ p\left| _{\eta =0} \right.=0 $$

(9b)

地震动沿长轴和短轴方向时，水体与结构交界面边界条件分别为：

$$ \frac{\partial p}{\partial \xi }\left| _{\xi ={{\xi }_{0}}} \right.=-\rho \ddot{u}b\text{cos}\eta $$

(10)

$$ \frac{\partial p}{\partial \xi }\left| _{\xi ={{\xi }_{0}}} \right.=-\rho \ddot{u}a\text{sin}\eta $$

(11)

式中，ρ表示水体密度，$\ddot{u}$表示结构的加速度。

1.2 分离变量求解

在椭圆柱坐标系下，$p(\xi, \eta, z)$可分离变量为：

$$ p=R(\xi)G(\eta)Z(z) $$

(12)

将式(12)代入式(4)整理得到3个解耦的方程：

$$ {Z}''+\lambda _{j}^{2}Z=0 $$

(13)

$$ {G}''+({{a}_{0}}+2q\text{cos}2\eta)G=0 $$

(14)

$$ {R}''-({{a}_{0}}+2q\text{cos}h2\xi)R=0 $$

(15)

式中，${{\lambda }_{j}}$和${{a}_{0}}$是分离变量常数，$q={{\mu }^{2}}\lambda _{j}^{2}/4$为一无量纲参数。

由方程(13)和边界条件式(5)、(6)可得：

$$ Z={{d}_{j}}\text{cos}{{\lambda }_{j}}z $$

(16)

式中，d_j为待定系数，${{\lambda }_{j}}=(2j-1)\text{ }\pi\text{ }/2h$，j=1，2，…。通过正交归一化可得${{d}_{j}}=\sqrt{2/h}$。

方程(14)为修正的角向马蒂厄方程(熊天信，2014)，其解为第一类角向马蒂厄函数$c{{e}_{n}}(\eta, -q)$和$s{{e}_{n}}(\eta, -q)$，即：

$$ c{{e}_{2n}}(\eta, -q)={{(-1)}^{n}}\sum\limits_{k=0}^{\infty }{{{(-1)}^{k}}A_{2k}^{(2n)}}\text{cos}2k\eta $$

(17a)

$$ c{{e}_{2n+1}}(\eta, -q)={{(-1)}^{n}}\sum\limits_{k=0}^{\infty }{{{(-1)}^{k}}B_{2k+1}^{(2n+1)}}\text{cos}(2k+1)\eta $$

(17b)

$$ s{{e}_{2n+2}}(\eta, -q)={{(-1)}^{n}}\sum\limits_{k=0}^{\infty }{{{(-1)}^{k}}B_{2k+2}^{(2n+2)}}\text{sin}(2k+2)\eta $$

(17c)

$$ s{{e}_{2n+1}}(\eta, -q)={{(-1)}^{n}}\sum\limits_{k=0}^{\infty }{{{(-1)}^{k}}A_{2k+1}^{(2n+1)}}\text{sin}(2k+1)\eta $$

(17d)

函数$c{{e}_{n}}(\eta, q)$和$s{{e}_{n}}(\eta, q)$的归一化正交关系为：

$$ \begin{array}{l} \int_0^{2\pi } {c{e_m}(\eta ,q)c{e_n}(\eta ,q){\rm{d}}\eta } = \left\{ \begin{array}{l} \pi ,m = n\\ 0,m \ne n \end{array} \right.\\ \int_0^{2\pi } {s{e_m}(\eta ,q)s{e_n}(\eta ,q){\rm{d}}\eta } = \left\{ \begin{array}{l} \pi ,m = n\\ 0,m \ne n \end{array} \right. \end{array} $$

(18)

$$ \begin{array}{l} \int_0^{2{\rm{\pi }}} {c{e_n}(\eta ,q){\rm{cos}}k\eta {\rm{d}}\eta = {\rm{\pi }}A_k^{(n)}} \\ \int_0^{2{\rm{\pi }}} {s{e_n}(\eta ,q){\rm{sin}}k\eta {\rm{d}}\eta = {\rm{\pi }}B_k^{(n)}} \end{array} $$

(19)

方程(15)为修正的径向马蒂厄方程(熊天信，2014)。整数阶径向马蒂厄方程的完全解为：

$$ R=\left\{ \begin{align} & \sum\limits_{n=0}^{\infty }{C_{n}^{1}I{{e}_{n}}}(\xi, -q)+C_{n}^{2}K{{e}_{n}}(\xi, -q) \\ & \sum\limits_{n=0}^{\infty }{C_{n+1}^{3}I{{o}_{n+1}}}(\xi, -q)+C_{n+1}^{4}K{{o}_{n+1}}(\xi, -q) \\ \end{align} \right. $$

(20)

式中，$C_{n}^{1}$、$C_{n}^{2}$、$C_{n}^{3}$和$C_{n}^{4}$是任意常数；函数$I{{e}_{n}}(\xi, -q)$和$I{{o}_{n}}(\xi, -q)$称为第一类变形贝塞尔型径向马蒂厄函数，是单调递增函数；函数$K{{e}_{n}}(\xi, -q)$和$K{{o}_{n}}(\xi, -q)$称为第二类变形贝塞尔型径向马蒂厄函数，是单调递减函数。

1.3 柱体的动水压力

当地震动沿长轴方向时，根据边界条件式(5)—(8)和(10)，并利用式(18)的正交性可得式(12)的解为：

$$ p=-\rho b\sum\limits_{j=1}^{\infty }{{{u}_{j}}\frac{B_{1}^{(1)}K{{e}_{1}}(\xi, -q)}{K{{e}_{1}}^{\prime }({{\xi }_{0}}, -q)}c{{e}_{1}}(\eta, -q)Z} $$

(21)

式中，${{u}_{j}}=\int_{0}^{h}{{\ddot{u}}}Z\text{d}z$，$K{{e}_{1}}^{\prime }(\xi, -q)$为函数$K{{e}_{1}}(\xi, -q)$的一阶导数。

当地震动沿短轴方向时，根据边界条件式(5)—(7)、(9)和(11)，并利用式(18)的正交性可得式(12)的解为：

$$ p=-\rho a\sum\limits_{j=1}^{\infty }{{{u}_{j}}\frac{B_{1}^{(1)}K{{o}_{1}}(\xi, -q)}{K{{o}_{1}}^{\prime }({{\xi }_{0}}, -q)}s{{e}_{1}}(\eta, -q)Z} $$

(22)

式中，$K{{o}_{1}}^{\prime }(\xi, -q)$为函数$K{{o}_{1}}(\xi, -q)$的一阶导数。

地震作用沿长轴方向时，椭圆截面柱体表面单位高度上的动水力为：

$$ {{f}_{x}}(z)=-\int_{0}^{2\pi }{p({{\xi }_{0}}, \eta, z)}b\text{cos}\eta \text{d}\eta $$

(23)

将式(21)代入式(23)整理得：

$$ {{f}_{x}}(z)=-\rho \text{ }\pi\text{ }{{b}^{2}}\sum\limits_{j=1}^{\infty }{{{u}_{j}}{{S}_{xj}}}Z $$

(24a)

$$ {{S}_{xj}}=-\frac{{{[B_{1}^{(1)}]}^{2}}K{{e}_{1}}({{\xi }_{0}}, -q)}{K{{e}_{1}}^{\prime }({{\xi }_{0}}, -q)} $$

(24b)

地震作用沿短轴方向时，椭圆截面柱体表面单位高度上的动水力为：

$$ {{f}_{y}}(z)=-\int_{0}^{2\pi }{p({{\xi }_{0}}, \eta, z)}a\text{sin}\eta \text{d}\eta $$

(25)

将式(22)代入式(25)整理得：

$$ {{f}_{y}}(z)=-\rho \text{ }\pi\text{ }{{a}^{2}}\sum\limits_{j=1}^{\infty }{{{u}_{j}}{{S}_{yj}}}Z $$

(26a)

$$ {{S}_{yj}}=-\frac{{{[A_{1}^{(1)}]}^{2}}K{{o}_{1}}({{\xi }_{0}}, -q)}{K{{o}_{1}}^{\prime }({{\xi }_{0}}, -q)} $$

(26b)

1.4 方法验证

首先，通过圆形桥墩动水力的解析解(Li等，2013)验证本文提出的椭圆形截面柱体动水压力的解析解。假定结构为刚性，将单位高度动水力沿高度积分，可得到椭圆形柱体的均布附加质量系数为：

$$ {{C}_{M}}=-\frac{1}{h{{m}_{0}}{{{\ddot{u}}}_{g}}}\int_{0}^{h}{f(z)\text{dz}} $$

(27)

式中，${{m}_{0}}$为单位高度水体的附加质量，沿长轴方向时${{m}_{0}}=\rho \text{ }\pi\text{ }{{b}^{2}}$，沿短轴方向时${{m}_{0}}=\rho \text{ }\pi\text{ }{{a}^{2}}$。图 3为本文圆柱体附加质量系数的解与圆柱解析解的对比。由图 3可以看出，本文解与圆柱解析解很好地吻合。

图 3 圆柱解析解与本文解的对比

Figure 3. The proposed solution compared with the analytical solution for circular cylinders and this paper

下载: 全尺寸图片幻灯片

进一步通过Wang等(2019)提出的数值方法验证本文提出的椭圆形截面柱体动水压力的解析解，该方法实质上是用垂向特征函数展开，将三维辐射问题简化为二维问题；然后用外域特征函数展开解与结构截面附近内域有限元联合求解。图 4为本文椭圆柱体附加质量系数的解与数值解的对比，椭圆柱的长轴和短轴尺寸分别为a=20m、b=10m。由图 4可以看出，本文的解与椭圆柱数值解很好地吻合。

图 4 椭圆柱数值解与本文解的对比

Figure 4. The proposed solution compared with the numerical solution for elliptical cylinders and this paper

下载: 全尺寸图片幻灯片

2. 水中椭圆柱体结构的地震反应求解

采用有限元方法将柔性柱体结构离散为梁单元，则地震作用下水中椭圆柱体的动力方程为(王勖成，2003)：

$$ {\bf{M}}{{{\bf{\ddot{u}}}}_{s}}(t)\text{+}{\bf{C}}{{{\bf{\dot{u}}}}_{s}}(t)+{\bf{K}}{{{\bf{u}}}_{s}}(t)=-{\bf{M}}{{{\bf{\ddot{u}}}}_{g}}(t)+{\bf{F}} $$

(28)

式中，M、C和K分别为结构质量、阻尼和刚度矩阵；${{\bf{u}}_{\rm{s}}}$为结构柔性位移列向量；F为地震动水力列向量。将地震动水力式(24)或(26)进行有限元离散，则F可以表示为

$$ {\bf{F}}\text{=}-{{{\bf{M}}}_{p}}[{{{\bf{\ddot{u}}}}_{g}}(t)+{{{\bf{\ddot{u}}}}_{s}}(t)] $$

(29)

$$ {{{\bf{M}}}_{p}}\text{=}{\bf{W}}\left[\sum\limits_{j=1}^{\infty }{{{m}_{0}}S_{j}^{{}}{\bf{Z}}_{{}}^{\text{T}}{\bf{Z}}} \right]\bf{W} $$

(30)

$$ {\bf{W}}\text{=}\int_{0}^{h}{{{{\bf{N}}}^{\text{T}}}{\bf{N}}\text{d}}z $$

(31)

式中，M_p为附加质量矩阵，N为形函数列向量；S_j如式(24b)或(26b)。

将式(29)代入式(28)整理得：

$$ ({\bf{M}}+{{{\bf{M}}}_{\text{p}}}){{{\bf{\ddot{u}}}}_{\text{s}}}+{\bf{C}}{{{\bf{\dot{u}}}}_{\text{s}}}+{\bf{K}}{{{\bf{u}}}_{\text{s}}}=-({\bf{M}}+{{{\bf{M}}}_{\text{p}}}){{{{\bf{\ddot{u}}}}}_{\text{g}}} $$

(32)

式(32)可通过数值积分方法Newmark-β方法求解(刘晶波等，2005)。需要指出的是，附加质量矩阵M_p是满阵的，难以在商业有限元中实现。本文将附加质量矩阵M_p的每一行元素进行集中化形成1个集中的附加质量矩阵M_g，即：

$$ {{\bf{M}}_{\text{p}}}=\left[\begin{array}{*{35}{l}} {{m}_{1, 1}} & {{m}_{1, 2}} & \cdots & {{m}_{1, L-1}} & {{m}_{1, L}} \\ {{m}_{2, 1}} & {{m}_{2, 2}} & \cdots & {{m}_{2, L-1}} & {{m}_{2, L}} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ {{m}_{L-1, 1}} & {{m}_{L-1, 2}} & \cdots & {{m}_{L-1, L-1}} & {{m}_{L-1, L}} \\ {{m}_{L, 1}} & {{m}_{L, 2}} & \cdots & {{m}_{L, L-1}} & {{m}_{L, L}} \\ \end{array} \right] $$

(33)

$$ {{\bf{M}}_{\text{g}}}=\left[\begin{array}{*{35}{l}} {{m}_{1}} & 0 & 0 & 0 & 0 \\ 0 & {{m}_{2}} & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & {{m}_{L-1}} & 0 \\ 0 & 0 & \cdots & 0 & {{m}_{L}} \\ \end{array} \right] $$

(34)

式中，${{m}_{i, j}}\text{=}{{m}_{j, i}}$，${{m}_{i}}={{m}_{i, 1}}+{{m}_{i, 2}}+\cdots \cdots +{{m}_{i, L-1}}+{{m}_{i, L}}$，L表示柱体结构水下的节点数目。

通过柔性椭圆柱体验证提出的集中附加质量矩阵的精度。图 5为地面运动位移，脉冲持时0.2s，频谱覆盖了地震作用频段。椭圆柱体尺寸为H=80m、a=40m、b=20m、a₁=30m、b₁=15m；密度和弹性模量分别为2500kg/m³和30GPa；水深h=80m；梁单元长8m；不考虑阻尼作用。图 6为地面运动沿长轴方向时，采用精确附加质量模型和集中附加质量模型计算得到的结构顶部位移时程。由图可见，采用集中附加质量方法计算的位移反应的周期偏大，即结构柔性运动引起的水体附加质量偏大。

图 5 地面位移时程和傅里叶谱

Figure 5. The time history of displacement of the ground motion and its Fourier spectrum

下载: 全尺寸图片幻灯片

图 6 参考解和集中附加质量模型位移时程的对比

Figure 6. The comparison of time history of displacement of the cylinder obtained by reference solution and lumped added mass method

下载: 全尺寸图片幻灯片

因此，本文进一步提出修正的集中附加质量方法，即将式(32)修正为如下形式：

$$ ({\bf{M}}+\alpha {{{\bf{M}}}_{\text{g}}}){{{\bf{\ddot{u}}}}_{\text{s}}}+{\bf{C}}{{{\bf{\dot{u}}}}_{\text{s}}}+{\bf{K}}{{{\bf{u}}}_{\text{s}}}=-({\bf{M}}+{{{\bf{M}}}_{\text{g}}}){{{\bf{\ddot{u}}}}_{\text{g}}} $$

(35)

式中，$\alpha \le 1$为集中附加质量矩阵修正系数。图 7为集中附加质量矩阵修正系数随无量纲参数宽深比(l)和长短轴比(δ)的变化；无量纲参数宽深比(l)和长短轴比(δ)的定义如公式(36)和(37)所示：

$$ l=\frac{D}{h} $$

(36)

$$ \delta =\frac{a}{b} $$

(37)

图 7 集中附加质量修正系数

Figure 7. Correction factors for the lumped added mass matrix

下载: 全尺寸图片幻灯片

式中，地震方向沿x轴时$D=2a$，地震方向沿y轴时$D=2b$。

图 8 参考解和修正集中附加质量模型的位移时程比较

Figure 8. The comparison of time history of displacement of the cylinder obtained by reference solution and modified-lumped added mass method

下载: 全尺寸图片幻灯片

图 8为地面运动沿长轴方向时，采用满阵附加质量模型和修正集中附加质量模型计算得到的结构顶部位移时程，满阵附加质量模型解为参考解。由图可见，采用修正集中附加质量方法计算的位移反应与参考解很好地吻合。

3. 椭圆形柱体均布附加质量简化公式

式(34)所示的集中附加质量矩阵与通过刚性结构动水力解析公式计算得到的附加质量是一致的。需要指出的是，动水力附加质量解析公式的数学表达复杂，难以在工程中进行应用。因此，提出附加质量的简化计算公式。假定附加质量沿高度均匀分布，相应的附加质量系数定义如式(27)所示。将式(24a)或(26a)代入式(27)整理得：

$$ {{C}_{M}}=\sum\limits_{j=1}^{\infty }{\frac{8{{S}_{j}}}{{{(2j-1)}^{2}}{{\text{ }\pi\text{ }}^{2}}}} $$

(38)

式中，$S_{j}^{{}}$的定义见式(24b)或(26b)。需要指出的是，简化公式是在$0.2\le l\le 2$和$0.2\le \delta \le 5$的范围内拟合得到的。

通过曲线拟合，圆形柱体均布附加质量系数的简化公式为：

$$ {{C}_{M1}}=0.6{{\text{e}}^{-0.93l}}+0.403{{\text{e}}^{-0.156l}} $$

(39)

当地震动沿长轴方向时，通过曲线拟合得到椭圆形柱体沿长轴方向的均布附加质量系数的简化公式为：

$$ {{C}_{Mx}}/{{C}_{M1}}={{p}_{11}}{{\delta }^{2}}+{{p}_{12}}\delta +{{p}_{13}} $$

(40)

$$ {{p}_{11}}=0.00367{{l}^{1.554}}+0.0221 $$

(41a)

$$ {{p}_{12}}=-0.185{{l}^{0.507}}-0.041 $$

(41b)

$$ {{p}_{13}}=0.157{{l}^{0.505}}+1.037 $$

(41c)

当地震动沿短轴方向时，通过曲线拟合，得到椭圆形柱体沿短轴方向的均布附加质量系数的简化公式为：

$$ {{C}_{My}}/{{C}_{M1}}={{p}_{21}}{{\delta }^{{{p}_{22}}}}+{{p}_{23}} $$

(42)

$$ {{p}_{21}}=-0.277{{\text{e}}^{-0.0186l}}+0.293{{\text{e}}^{-1.102l}} $$

(43a)

$$ {{p}_{22}}=-0.008{{l}^{2}}+0.186l-1.056 $$

(43b)

$$ {{p}_{23}}=1.295{{{\rm{e}}}^{-0.0106l}}-0.31{{{\rm{e}}}^{-1.052l}} $$

(43c)

图 9为椭圆柱体附加质量系数的解析解与简化公式的对比，图 10则为简化公式的误差，可以看出简化公式与解析解吻合较好。

图 9 附加质量系数解析解和简化公式的对比

Figure 9. The comparison of coefficient of the added mass obtained by analytical solution and simplified formula

下载: 全尺寸图片幻灯片

图 10 附加质量系数简化公式的误差

Figure 10. The error of the simplified formula for the added mass coefficient

下载: 全尺寸图片幻灯片

4. 结语

(1) 基于椭圆坐标系，采用分离变量法将拉普拉斯方程转换为马蒂厄方程。通过求解马蒂厄方程，提出了椭圆柱体结构地震动水压力的解析解。

(2) 建立了地震作用下椭圆柱体结构与水体相互作用的动力有限元方程，结构的动水力通过附加质量矩阵施加，该矩阵是满阵的。

(3) 为便于椭圆柱体结构附加质量矩阵在商业有限元中实现，提出了集中附加质量矩阵的方法，该方法中柔性运动引起的附加质量为集中附加质量矩阵和修正系数的乘积。

(4) 基于刚性柱体结构动水力的解析解，通过曲线拟合的方法建立了椭圆柱体结构动水力的均布附加质量简化公式，公式中的系数仅与无量纲参数宽深比和长短轴比相关。

图 1 地震目录震源深度直方图

Figure 1. Histogram of focal depth in earthquake catalog

下载: 全尺寸图片幻灯片

图 2 “3400”走时表速度模型

Figure 2. Velocity model from"3400"ravel time table

下载: 全尺寸图片幻灯片

图 3 虚拟台网及虚拟震中位置

Figure 3. Virtual network and suppositional earthquake location

下载: 全尺寸图片幻灯片

图 4 新疆测震台网分布及虚拟震中位置

Figure 4. Distribution of seismic station in Xinjiang and suppositional earthquake location

下载: 全尺寸图片幻灯片

表 1 虚拟台网不同初始深度的定位结果

Table 1. Location results at different initial depths of the virtual network

虚拟震中	虚拟发震时刻	初始深度 /km	定位结果
虚拟震中	虚拟发震时刻	初始深度 /km	北纬/°	东经/°	深度/km	发震时刻	残差
40.0719°N 116.0939°E H：13km	2010-05-08 10:30:30.0	10	40.0839	116.0902	11	10:30:30.3	0.2
		20	40.0672	116.0916	13	10:30:29.8	0.06
		30	40.0713	116.0934	12	10:30:30.0	0.006
		40	40.0713	116.0939	13	10:30:29.9	0.01

下载: 导出CSV

表 2 各研究区域不同初始深度的定位结果

Table 2. Results of location in different regions at different initial depths of Xinjiang Digital Seismic Network

区域	虚拟震中	虚拟发震时刻	初始深度/km	定位结果
区域	虚拟震中	虚拟发震时刻	初始深度/km	北纬/°	东经/°	深度/km	发震时刻	残差
南天山西段	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30 40	40.54 40.54 40.54 40.54	79.65 79.65 79.65 79.65	9 9 10 9	23:22:3.0 23:22:3.0 23:22:3.1 23:22:3.0	0.21 0.22 0.25 0.21
天山中段	43.83°N 86.35°E H：10km	2016-12-08 13:15:03	10 20 30 40	43.82 43.82 43.82 43.82	86.35 86.35 86.35 86.35	10 10 10 9	13:15:3.1 13:15:3.1 13:15:3.1 13:15:3.1	0.46 0.46 0.46 0.46
西昆仑	36.12°N 82.49°E H：10km	2014-02-12 17:19:48	10 20 30 40	36.13 36.12 36.14 36.13	82.49 82.49 82.51 82.49	8 9 10 9	17:19:48.3 17:19:48.4 17:19:48.5 17:19:48.4	0.37 0.38 0.47 0.38

下载: 导出CSV

表 3 加入随机误差后不同初始深度的定位结果

Table 3. Positioning results at different initial depths after adding random errrors

区域	虚拟震中	虚拟发震时刻	初始深度/km	定位结果
区域	虚拟震中	虚拟发震时刻	初始深度/km	北纬/°	东经/°	深度/km	发震时刻	残差
南天山西段	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30 40	40.54 40.54 40.54 40.54	79.64 79.65 79.66 79.68	11 11 10 9	23:22:3.1 23:22:3.1 23:22:3.1 23:22:2.9	0.41 0.41 0.42 0.45
天山中段	43.83°N 86.35°E H：10km	2016-12-08 13:15:03	10 20 30 40	43.83 43.83 43.82 43.82	86.36 86.36 86.36 86.35	12 12 10 11	13:15:3.1 13:15:3.1 13:15:3.1 13:15:3.0	0.58 0.56 0.52 0.55
西昆仑	36.12°N 82.49°E H：10km	2014-02-12 17:19:48	10 20 30 40	36.12 36.12 36.12 36.13	82.49 82.49 82.51 82.51	7 8 9 8	17:19:48.3 17:19:48.3 17:19:48.2 17:19:48.4	0.49 0.42 0.46 0.46

下载: 导出CSV

表 4 不同数据源的定位结果

Table 4. Positioning results from different data source

处理方法	虚拟震中	虚拟发震时刻	初始深度/km	定位结果
处理方法	虚拟震中	虚拟发震时刻	初始深度/km	北纬/°	东经/°	深度/km	发震时刻	残差
1	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30	40.54 40.54 40.53	79.65 79.64 79.66	8 11 12	23:22:2.9 23:22:3.2 23:22:2.0	0.4 0.43 0.41
2	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30	40.54 40.53 40.54	79.65 79.66 79.63	8 8 9	23:22:3.0 23:22:2.8 23:22:3.0	0.4 0.38 0.25
3	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30	40.56 40.56 40.55	79.64 79.63 79.64	23 19 21	23:22:2.0 23:22:1.5 23:22:1.7	1.3 1.3 1.2
4	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30	40.55 40.53 40.51	79.64 79.70 79.63	10 10 10	23:22:2.9 23:22:2.0 23:22:3.4	0.35 0.65 0.53
5	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30	40.54 40.55 40.55	79.60 79.64 79.63	5 4 9	23:22:2.8 23:22:2.9 23:22:3.6	0.27 0.28 0.19
6	40.54°N 79.65°E H：10km	2010-01-02 23:22:03	10 20 30	40.50 40.50 40.49	79.39 79.38 79.37	5 6 6	23:22:4.5 23:22:5.8 23:22:5.7	0.69 0.6 0.54

下载: 导出CSV

表 5 虚拟震源位于下地壳时不同数据源的定位结果

Table 5. Positioning results from different data source in which the suppositional earthquake is located in the lower crust

处理方法	虚拟震中	虚拟发震时刻	初始深度 /km	定位结果
处理方法	虚拟震中	虚拟发震时刻	初始深度 /km	北纬/°	东经/°	深度/km	发震时刻	残差
01	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.54 40.54 40.54	79.65 79.65 79.65	30 29 29	23:22:3.0 23:22:2.9 23:22:2.9	0.06 0.2 0.12
02	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.54 40.53 40.54	79.65 79.65 79.65	29 31 29	23:22:2.9 23:22:3.0 23:22:2.9	0.31 0.29 0.28
1	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.54 40.53 40.54	79.67 79.66 79.65	18 25 27	23:22:2.3 23:22:2.8 23:22:2.8	0.9 0.4 0.5
2	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.53 40.54 40.54	79.66 79.65 79.64	10 31 30	23:22:3.4 23:22:2.9 23:22:3.0	0.6 0.3 0.2
3	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.55 40.56 40.54	79.64 79.63 79.64	39 30 41	23:22:1.6 23:22:1.0 23:22:1.8	0.9 1.4 1.0
4	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.54 40.51 40.56	79.61 79.59 79.64	31 31 31	23:22:3.4 23:22:3.9 23:22:2.9	0.3 0.3 0.3
5	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.59 40.54 40.55	79.70 79.62 79.63	27 30 33	23:22:1.8 23:22:3.3 23:22:3.2	0.5 0.2 0.3
6	40.54°N 79.65°E H：30km	2010-01-02 23:22:03	10 20 30	40.53 40.49 40.47	79.46 79.45 79.34	37 36 36	23:22:4.0 23:22:4.3 23:22:4.9	0.4 0.5 0.7
注：处理方法01未进行任何处理，处理方法02加入±0.5s随机误差，下同。

下载: 导出CSV

表 6 不同虚拟震源深度的定位结果（固定初始震源深度10km）

Table 6. Positioning results in different virtual source locations at a fixed initial depth of 10km

震源位置及发震时刻	处理方法	定位结果
震源位置及发震时刻	处理方法	北纬/°	东经/°	深度/km	发震时刻	残差
40.54°N 79.65°E H：5km 2010-01-02 23:22:03	01 02 1 2 3 4 5 6	40.54 40.54 40.53 40.53 40.54 40.57 40.53 40.50	79.65 79.65 79.65 79.66 79.64 79.65 79.67 79.43	5 5 6 7 20 6 7 19	23:22:3.1 23:22:3.1 23:22:3.0 23:22:3.0 23:22:2.9 23:22:2.8 23:22:2.7 23:22:4.9	0.2 0.3 0.5 0.6 1.6 0.5 0.1 0.6
40.54°N 79.65°E H：15km 2010-01-02 23:22:03	01 02 1 2 3 4 5 6	40.54 40.53 40.54 40.54 40.55 40.55 40.53 40.52	79.65 79.65 79.65 79.66 79.63 79.64 79.65 79.51	15 15 16 12 8 15 10 15	23:22:3.0 23:22:3.0 23:22:3.1 23:22:3.0 23:22:1.7 23:22:3.0 23:22:3.2 23:22:3.9	0.2 0.3 0.2 0.4 1.9 0.3 0.2 0.6
40.54°N 79.65°E H：20km 2010-01-02 23:22:03	01 02 1 2 3 4 5 6	40.54 40.53 40.54 40.53 40.56 40.55 40.53 40.52	79.65 79.65 79.65 79.66 79.64 79.62 79.63 79.48	20 19 19 17 16 20 14 24	23:22:2.9 23:22:3.0 23:22:3.0 23:22:3.1 23:22:2.0 23:22:3.2 23:22:3.2 23:22:4.2	0.3 0.47 0.4 0.3 1.7 0.4 0.2 0.7

下载: 导出CSV

参考文献(19)

蔡明军, 山秀明, 徐彦等, 2004.从误差观点综述分析地震定位方法.地震研究, 27(4):314-317. doi: 10.3969/j.issn.1000-0666.2004.04.005

陈向军, 上官文明, 宋秀青等, 2014.新疆全区和分区地壳速度模型的分析.中国地震, 30(2):178-187. doi: 10.3969/j.issn.1001-4683.2014.02.005

李艳永, 热依木江, 唐明帅等, 2016.利用震相方位角改善地震台网稀疏地区地震定位精度.地震地磁观测与研究, 37(2):57-62. http://d.old.wanfangdata.com.cn/Periodical/dzdcgcyyj201602010

田玥, 陈晓非, 2002.地震定位研究综述.地球物理学进展, 17(1):147-155. doi: 10.3969/j.issn.1004-2903.2002.01.022

田玥, 陈晓非, 2005.水平层状介质中的快速两点间射线追踪方法.地震学报, 27(2):147-154. doi: 10.3321/j.issn:0253-3782.2005.02.004

王桂丹, 2016.基于JOPENS/MSDP地震定位方法的定位结果分析.内陆地震, 30(1):74-80. http://www.cnki.com.cn/Article/CJFDTotal-LLDZ201601009.htm

夏仕安, 张佑龙, 隆爱军, 2011. MSDP在测震台网的应用.四川地震, (3):13-16. doi: 10.3969/j.issn.1001-8115.2011.03.003

谢辉, 金春华, 蔡新华等, 2011. MSDP单纯型和HYP2000定位方法对比分析.地震地磁观测与研究, 32(5):15-19. doi: 10.3969/j.issn.1003-3246.2011.05.003

新疆地震局分析预报室, 1982. "3000公里走时表"及其应用、新疆走时表工作的设想.地震地磁观测与研究, 3(4):25-30. http://www.cnki.com.cn/Article/CJFDTotal-DZGJ198204005.htm

曾融生, 滕吉文, 阚荣举等, 1965.我国西北地区地壳中的高速夹层.地球物理学报, 14(2):94-106. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX196502002.htm

张炳, 汪贵章, 戚浩等, 2012. HYPO2000与单纯形定位方法对比及作为速报算法的可行性分析.地震地磁观测与研究, 33(5-6):13-17. http://d.old.wanfangdata.com.cn/Periodical/dzdcgcyyj201205003

张诚, 张伶, 邓齐赞等, 1979.甘肃及邻近地区的地壳厚度.西北地震学报, 1(2):22-26. http://www.cnki.com.cn/Article/CJFDTOTAL-ZBDZ197902003.htm

张志斌, 朱浩清, 李艳永, 2015.基于地震波反演研究南天山中西段的震源深度.震灾防御技术, 10(S1):701-711. doi: 10.11899/zzfy2015s101

赵珠, 丁志峰, 易桂喜等, 1994.西藏地基定位-一种使用单纯形优化的非线性方法.地震学报, 16(2):212-219. http://www.cnki.com.cn/Article/CJFD1994-DZXB402.010.htm

赵仲和, 牟磊育, 2010.地震定位结果中误差描述的标准化.见: 中国地球物理学会第二十六届年会暨中国地震学会第十三次学术大会论文集.宁波: 中国地球物理学会, 中国地震学会, 299. http://cpfd.cnki.com.cn/Article/CPFDTOTAL-ZGEM201010001267.htm

中国地震局监测预报司, 2017.测震学原理与方法.北京:地震出版社.

Geiger L., 1912. Probability method for the determination of earthquake epicenters from arrival time only. Bulletin of Saint Louis University, (8):60-71.

Nelder J. A., Mead R., 1965. A simplex method for function minimization. The Computer Journal, 7(4):308-313. doi: 10.1093/comjnl/7.4.308

Prugger A. F., Gendzwill D. J., 1988. Microearthquake location:a nonlinear approach that makes use of a simplex stepping procedure. Bulletin of the Seismological Society of America, 78(2):799-815.

施引文献

资源附件(0)

访问统计