Signature and physical mechanism of borehole dilatometers in response to super typhoon Lekima in southeastern coastal area of China
-
摘要: 依据超强台风 “利奇马” 的强度和时空演变特征，本文采用经验模态分解等方法系统地分析并揭示了该台风对我国东南沿海地区钻孔体应变影响的全貌，并在此基础上对台风扰动的机制进行了初步探讨。结果表明：① 台风演变过程中漏斗状的长周期气压波动，是造成钻孔体应变大幅张性变化的物理成因，且体应变对台风低压系统具有即时的线弹性响应特征，其变化形态与气压漏斗高度相似，弹性变形的持续时间与气压波动的历时较一致；② 在周期为103 h时，−18.2 hPa的气压变幅便可在地下62 m深处产生高达−112.1×10−9的体应变，该频点的气压影响系数为6.2×10−9/hPa；③ 在空间上，台风中心在980 km以外便能影响体应变观测，且随着台风的不断逼近或远离，其影响程度也相应地逐渐增强或减弱。
-
关键词:
- 超强台风 “利奇马” /
- 钻孔体应变 /
- 东南沿海 /
- 气压负荷 /
- 台风暴雨
Abstract: The aim of this study is to investigate how volumetric strain and super typhoon Lekima interact with each other. By combining the intensity and spatio-temporal evolution of Lekima, we systematically analyze the extracted volumetric deformation of shallow crust at five stations in the southeastern coastal area of China induced by Lekima using a decomposition analysis called the empirical mode decomposition （EMD）. Furthermore, we discuss the physical mechanism of typhoon-induced volumetric strain. The results show that: ① The super typhoon’s signature consists in a ground dilatation due to barometric pressure drop drastically, generally followed by a ground compression due to the barometric pressure recovery. The dynamic patterns for the volumetric strain and the barometric pressure are both similar to the symme-trical funnel, and the response of volumetric strain to the barometric pressure is almost instant and linear. In addition, the durations are nearly equivalent for both of them. ② The maximum magnitude in volumetric strain induced by barometric pressure with fluctuation of −18.2 hPa can reach up to −112.1×10−9 in the 62 m deep borehole, and the corresponding coefficient of barometric pressure response is about 6.2×10−9/hPa at the 103-hour period. ③ Obviously, the dilatometer records can be remotely disturbed by super typhoon at a distance of approximately 980 km away from the borehole site. When the typhoon center approaches or moves away from the borehole site, the dilatational magnitude of borehole strain will increase or decrease correspondingly. The obtained results can be used not only for identifying and determining reasonably the physical mechanism of anomalous changes induced by typhoon in low-frequency range for southeastern coastal and inland areas in China, but also for contributing the reliably observational evidences to the theoretical model research for the volumetric strain in response to barometric pressure in the low-frequency range. -
引 言
在地震波处理和成像中，快速且准确的射线追踪和走时计算是至关重要的，尤其在走时层析成像（Klibanov et al，2023）、微震源定位（Grechka et al，2015）和地震叠前偏移（Garabito，2021）等应用中起着关键作用。随着技术的不断发展，基于各向同性介质假设的技术已经不能完全满足实际生产应用的需要。各向异性介质射线追踪方法由各向同性方法发展而来，分为两大类：一是传统的两点间射线追踪方法（Julian，Gubbins，1977），该类方法可达到很高的计算精度，但对于复杂介质可能存在阴影区或走时局部最小问题（王东鹤等，2016）；二是基于网格的射线追踪方法，主要包括最短路径法和有限差分解程函方程法。最短路径法具有良好的稳健性，但为保证计算精度，需要较细的网格，致使其计算成本较高（Nakanishi，Yamaguchi，1986；Zhou，Greenhalgh，2005；张美根等，2006；赵爱华等，2006；赵后越，张美根，2014；邴琦等，2020）。相比之下，有限差分求解程函方程法可有效弥补传统射线追踪法的不足，在不增加网格密度的情况下提供高效且较为精确的数值（周洪波，张关泉，1994；赵烽帆等，2014）。
快速扫描法（fast sweeping method，缩写为FSM）和快速推进法（fast marching method，缩写为FMM）是有限差分求解程函方程最常用的两种算法。FSM利用局部解程函方程来更新格点走时信息，通过反复扫描模型以稳定走时场，通常采用上风插值方案求局部解（Zhao，2004）。因此，FSM被广泛用于求解各向异性介质的qP波走时问题（Waheed et al，2015；Hao et al，2018；Tro et al，2023）。然而，FSM在复杂模型中需要更多的迭代次数以获得高精度的走时，较为耗时。Sethian （1999）提出的FMM则采用上风差分格式来求解局部程函方程，利用窄带延拓来重建走时波前，并使用堆排序技术保存走时信息，以加速寻找最小走时。这一方法显著提高了计算效率，将计算量从波前扩展法的O（N3）减小到O（NlgN）。由于各向异性群速度确定的能量传播方向与相速度给出的波前传播方向不同，原始的FMM只适用于各向同性介质，不能可靠地用于求解各向异性方程。为解决这一问题，Sethian和Vladimirsky （2003）提出了一种适用于各向异性介质的有序迎风方法。随后，许多学者提出了改进的FMM方法，以适用于各向异性介质，并开展了进一步的研究，使其得以快速发展。
Lelièvre等（2011）将FMM方法扩展到非结构三维四面体网格中进行初至地震走时计算。Treister和Haber （2016）提出了一种适用于因式分解程函方程的FMM算法，有效地解决了由震源奇异性引起的不精确数值解问题。Waheed和Alkhalifah （2017）利用因式分解方法研究了倾斜横向各向同性（tilted transverse isotropy，缩写为TTI）方程问题，通过迭代求解一个因式分解的倾斜椭圆各向异性程函方程序列，能够更好地描述各向异性波前。Sun等（2017）提出了一种联合3D的走时计算方法，在震源附近小范围内使用计算精度较高的波前构建法计算走时，在剩余区域使用FMM算法计算走时，改善了震源附近区域计算精度不高的问题。Le Bouteiller等（2019）在间断伽辽金方法（discontinuous Galerkin method，缩写为DGM）中引入FSM算法来求解三维程函方程，降低了计算的复杂度。Waheed （2020）将TTI介质程函方程表示为椭圆各向同性方程序列后利用FMM可快速准确地计算TTI介质的走时，但计算工作量随着迭代次数的增加而增大，与各向同性FMM相比，算法实现变得更加复杂。因此，如何在求解各向异性程函方程时均衡计算效率和精准度是目前走时计算的一个挑战性问题。
针对上述问题，本文基于各向同性FMM思想提出一种基于相速度预测的各向异性FMM射线追踪方法（phase velocity prediction fast marching method，缩写为PVP-FMM）。首先，预测网格节点的相速度，再使用FMM计算网格节点波前时间；其次，计算每个网格的波前梯度，重构相速度矢量场；进一步，由相速度矢量场重构群速度矢量场；最后，由群速度矢量场插值计算射线路径和走时。该方法继承了各向异性FMM无条件稳定的特点，具有与各向同性FMM相同的计算效率和精度，能够适用于复杂构造的TTI介质射线追踪，具有良好的稳定性和较高的计算效率。
1. 各向异性程函方程建立
1.1 TTI介质程函方程
各向异性介质在高频近似条件下的程函方程可表示为（Vidale，1988）
$$ v_{\mathrm{P}}^2 （ \theta ） \nabla \tau · \nabla \tau＝1\text{，} $$ (1) 式中，$\nabla $为梯度算子，τ为走时场计算中的当前走时，vP为各向异性介质准qP波的相速度。相速度与波前法向方向有关，Thomsen （1986）在弱各向异性条件下获得相速度的近似表达式为
$$ {v_{\mathrm{P}}} （ \theta ） ＝{v_{\mathrm{n}}} （ 1 ＋ \varepsilon {\sin ^4}\theta ＋ \delta {\sin ^2}\theta {\cos ^2}\theta ） \text{，} $$ (2) 式中，θ为相角，vn为地层法向方向速度，ε和δ都为汤姆森（Thomsen）参数。
Waheed （2020）将Alkhalifah （2000）提出的TTI介质程函方程分解成一系列椭圆方程，使用椭圆波前迭代拟合各向异性波前，但随着迭代次数增加的同时，也增加了计算成本。为此，本文直接利用式（1）求解各向异性波前，将其改写为
$$ {\left( {\frac{{ {\text{∂}} \tau }}{{ {\text{∂}} x}}} \right)^2} ＋ {\left( {\frac{{ {\text{∂}} \tau }}{{ {\text{∂}} {\textit{z}} }}} \right)^2}＝\frac{1}{{v_{\mathrm{n}}^2{{ （ 1 ＋ \varepsilon {{\sin }^4}\theta ＋ \delta {{\sin }^2}\theta {{\cos }^2}\theta ） }^2}}}＝f （ \theta ） . $$ (3) 从式（3）中可以看出，如果能够知道其相速度，就可计算出相应的走时。与各向同性程函方程求解方法一样，只需在每个网格点计算二次多项式，不用在每个网格点上求解复杂的高阶多项式，与各向同性时具有相同的计算效率。因此相速度精确预测是使用式（3）求解各向异性波前的关键环节。
1.2 FMM方法离散化处理程函方程
FMM算法将整个计算域划分为三个区域：上风区、波前区和下风区。波从上风区经过波前区传播至下风区，如图1所示。上风区为已经求解出走时的已知区域，波前区是正在求解的区域，下风区是还未求解的未知区域。在波前区的窄带内，寻找最小波前时间点，将其属性改为已知点的同时继续向下风区扩展，直至所有网格点都在上风区。
![]()
图 1 快速推进法（FMM）实现原理图Figure 1. Implementation schematic diagram of fast marching method （FMM）在图1中，设节点A是波前区最小波前时间节点，节点B是其邻近的走时未知节点。由于点A的波前时间已知，可以计算波前梯度求出点A的波传播方向。当波由点A传播至点B时，如果两个节点的速度不同，传播方向会发生变化，从而使得节点B的传播方向和相角也是未知的。因此，在节点B使用式（3）求解波前时间的关键就是要预测相速度。
使用FMM方法对式（3）进行二阶上风差分格式离散。对于网格点（i， j），离散如下：
$$ \begin{gathered} { ［ \max （ D_{i \text{，} j}^{ - x}\tau \text{，} - D_{i \text{，} j}^{ ＋ x}\tau \text{，} 0 ） ］ ^2} ＋ { ［ \max （ D_{i \text{，} j}^{ - {\textit{z}} }\tau \text{，} - D_{i \text{，} j}^{ ＋ {\textit{z}} }\tau \text{，} 0 ） ］ ^2}＝f （ {\theta _{i \text{，} j}} ） \text{，} \\ f （ {\theta _{i \text{，} j}} ） ＝\frac{1}{{v_{i \text{，} j}^2{{ （ 1 ＋ {\varepsilon _{i \text{，} j}}{{\sin }^4}{\theta _{i \text{，} j}} ＋ {\delta _{i \text{，} j}}{{\sin }^2}{\theta _{i \text{，} j}}{{\cos }^2}{\theta _{i \text{，} j}} ） }^2}}} \text{，} \\ \end{gathered} $$ (4) 式中，$v_{i \text{，} j} $，${\varepsilon _{i \text{，} j}} $，${\delta _{i \text{，} j}} $，${\theta _{i \text{，} j}} $分别为地层法向速度、汤姆森参数和相角在节点（i， j）的值，$D_{i \text{，} j}^{ - x} $和$D_{i \text{，} j}^{ ＋ x} $是待求点波前时间τ在x方向上的二阶差分运算符，表示为
$$ \begin{gathered} D_{i \text{，} j}^{ - x}\tau＝\frac{{3{\tau _{i \text{，} j}} - 4{\tau _{i - 1 \text{，} j}} ＋ {\tau _{i - 2 \text{，} j}}}}{{2\Delta x}} \text{，} D_{i \text{，} j}^{ ＋ x}\tau＝- \frac{{3{\tau _{i \text{，} j}} - 4{\tau _{i ＋ 1 \text{，} j}} ＋ {\tau _{i ＋ 2 \text{，} j}}}}{{2\Delta x}} \text{，} \end{gathered} $$ (5) 式中Δx是x方向上的网格间距。z方向的前、后二阶差分计算与x方向前后二阶差分计算类似。这种迎风差分方案可确保信息从较低的波前时间值传播到较高的值。
2. 基于相速度预测的FMM波前重构
在地震波传播过程中，当地层参数变化时，其传播方向与相速度也会发生变化。当求解待求点B的波前时间时，如果依然采用已知点A的相速度进行计算，则波前时间就会存在一定误差。若这种误差累积扩散到整个计算域，则导致波前时间解的精度降低。为解决此问题，本文利用斯奈尔（Snell）方程推导出待求点与已知点间的相速度关系，进而得出相速度预测公式；再结合相速度预测公式与FMM算法求解波前时间；最后，采用群速度矢量场插值计算射线路径。
2.1 TTI介质相速度预测
如图2所示，假设地层参数在节点B处存在一个倾斜分界面。对于TTI介质，分界面两边的vn，ε和δ不同，根据斯奈尔定律，在分界面满足如下关系（邓怀群等，2000）:
$$ \begin{split} p＝\frac{{\sin {\theta _{{\rm{inc}}}}}}{{{v_{\rm{P}}} （ {\theta _{{\mathrm{inc}}}} ） }}＝\frac{{\sin {\theta _{B}}}}{{{v_{\rm{P}}} （ {\theta _{B}} ） }} \text{，} \end{split} $$ (6) 式中，p为射线参数，θinc和θB分别为入射角和出射角，${{v_{\mathrm{P}}} （ {\theta _{{\mathrm{inc}}}} ） } $和${{v_{\mathrm{P}}} （ {\theta _{{B}}} ） } $分别为节点A和B处的相速度。式中出射角θB可表示为θB＝θinc+Δθ，令x＝sinθinc，式（6）可等价表示为
$$ p＝\frac{x}{{{v_{\mathrm{n}}} （ 1 ＋ \varepsilon {x^4} － \delta {x^4}＋ \delta {x^2} ） }}. $$ (7) 为了便于计算，先对式（7）进行处理，令$f （ x ） ＝1 ＋ \varepsilon {x^4} － \delta {x^4} ＋ \delta {x^2}$，可知
$$ p{v_{\mathrm{n}}}f （ x ） ＝x\text{，} $$ (8) 对式（8）两边进行微分可得
$$ pf （ x ） {\mathrm{d}}{v_{\mathrm{n}}} ＋ p{v_{\mathrm{n}}}{x^4}{\mathrm{d}}\varepsilon ＋ （ {x^2} － {x^4} ） {\mathrm{d}}\delta \cdot p{v_{\mathrm{n}}} ＋ p{v_{\mathrm{n}}} （ 4{x^3}\varepsilon － 4{x^3}\delta ＋ 2x\delta ） {\mathrm{d}}x＝{\mathrm{d}}x. $$ (9) 根据式（8）可知，$p{v_{\mathrm{n}}}＝{x}/{{f （ x ） }} \text{，} pf （ x ） ＝{x}/{{{v_{\mathrm{n}}}}}$，将其代入式（9）有
$$ x\frac{{{\mathrm{d}}{v_{\mathrm{n}}}}}{{{v_{\mathrm{n}}}}} ＋ \frac{x}{{f （ x ） }}{x^4}{\mathrm{d}}\varepsilon ＋ \frac{x}{{f （ x ） }}{x^2} （ 1 － {x^2} ） {\mathrm{d}}\delta＝\left\{ {1 － \frac{x}{{f （ x ） }} ［ {4{x^3} （ \varepsilon － \delta ） ＋ 2x\delta } ］ } \right\}{\mathrm{d}}x\text{，} $$ (10) 进一步化简后可得
$$ \frac{\mathrm{d}x}{x}＝\frac{f （ x ） }{1＋3 （ \delta－\varepsilon ） x^4－x^2\delta}\frac{\mathrm{d}v_{\mathrm{n}}}{v_{\mathrm{n}}}＋\frac{x^4}{1＋3 （ \delta－\varepsilon ） x^4－x^2\delta}\mathrm{d}\varepsilon＋\frac{x^2 （ 1－x^2 ） }{1＋3 （ \delta－\varepsilon ） x^4－x^2\delta}\mathrm{d}\delta. $$ (11) 根据式（6）可知，待求点B的相速度为
$$ v_{\mathrm{P}} （ \theta_{\mathrm{B}} ） ＝v_{\mathrm{P}} （ \theta_{\mathrm{inc}} ） \frac{\sin （ \theta_{\mathrm{inc}}＋\Delta\theta ） }{\sin\theta_{\mathrm{inc}}}＝v_{\mathrm{P}} （ \theta_{\mathrm{inc}} ） \left(1＋\frac{\Delta x}{x}\right)\text{，} $$ (12) 在求解过程中需满足$0 {\text{≤}} \sin （ {\theta_{{\mathrm{inc}}}} ＋ \Delta \theta ） ＝x ＋ \Delta x {\text{≤}} 1$的条件。将式（11）代入式（12）有
$$ v_{\mathrm{P}} （ \theta_{\mathrm{B}} ） ＝v_{\mathrm{P}} （ \theta_{\mathrm{inc}} ） \left[\frac{f （ x ） }{1＋3 （ \delta－\varepsilon ） x^4－x^2\delta}\frac{\mathrm{d}v_{\mathrm{n}}}{v_{\mathrm{n}}}＋\frac{x^4}{1＋3 （ \delta－\varepsilon ） x^4－x^2\delta}\mathrm{d}\varepsilon＋\frac{x^2 （ 1－x^2 ） }{1＋3 （ \delta－\varepsilon ） x^4－x^2\delta}\mathrm{d}\delta\right]\text{，} $$ (13) 其中入射角θinc满足
$$ \theta_{\mathrm{inc}}＝\beta_{\mathrm{\mathit{A}}}－\alpha_{\mathrm{\mathit{B}}}\text{，} $$ (14) 式中αB为节点B处的界面倾角，βA为节点A的传播方向与z轴的夹角。由于节点A的波前时间已是最小走时，不再变化，利用x方向与z方向上的∂τ/∂x与∂τ/∂z，传播方向βA表示为
$$ \mathrm{arctan}\beta_{\mathrm{\mathit{A}}}＝\frac{\text{∂}\tau}{\text{∂}x}\mathord{\left/\vphantom{\frac{\text{∂}\tau}{\text{∂}x}\frac{\text{∂}\tau}{\text{∂}\text{z}}}\right.}\frac{\text{∂}\tau}{\text{∂}\text{z}}. $$ (15) 将式（14）和式（15）代入式（13），可得到待求点B的相速度vP（θB）。
对于待求点B的相角θB，利用式（4）以及式（13），求得待求点B的走时，再利用式（15）得到其传播方向βB，通过${\theta _{{B}}}＝{\beta _{{B}}} － {\alpha _{{B}}}$获得点B的相角。
因此，PVP-FMM方法求解波前时间最终分为三个步骤：首先，由点A的波前梯度计算点A的传播方向，再利用式（14）计算入射角；其次，根据相速度预测公式，计算出点B的相速度和相角；最后，将相速度代入式（4）中，计算出该点的波前时间，然后进行下一点预估，依此类推，直至所有网格点完成计算。
2.2 群速度矢量场和射线路径计算
各向异性射线追踪方法可分为四个步骤：首先，使用基于相速度预测的各向异性FMM方法计算网格节点波前时间；其次，计算每个网格的波前梯度，重构相速度矢量场；进一步，由相速度矢量场重构群速度矢量场；最后，由群速度矢量场插值计算射线路径和走时。具体计算步骤如下：① 利用式（4）和式（13）进行各向异性PVP-FMM波前时间计算，得到所有网格点波前时间场；② 在波前时间场中利用梯度运算可确定其整个相速度场，即
$$ \nabla\tau＝\frac{\text{∂}\tau}{\text{∂}x}\boldsymbol{x}＋\frac{\text{∂}\tau}{\text{∂}\text{z}}\boldsymbol{z}，\boldsymbol{v}_{\mathrm{P}}＝v_{\mathrm{P}}\frac{\nabla\tau}{\left|\nabla\tau\right|}\text{，} $$ (16) 式中，$\left| {\nabla \tau } \right|＝1/{v_{\mathrm{P}}}(\theta )$，${{\nabla \tau }}/{{\left| {\nabla \tau } \right|}}＝{\boldsymbol{n}}\cos \theta ＋ {\boldsymbol{t}} \sin \theta $，${\boldsymbol{ n}}$和$ {\boldsymbol{t}}$表示地层切向和法向的单位矢量；③ 通过群速度与相速度的关系可得到群速度场（Thomsen，1986），即
$$ {{\boldsymbol{v}}} _{\mathrm{g}}＝{{{\boldsymbol{v}}} _{\mathrm{P}}} （ \theta ） ＋ \frac{{{\mathrm{d}}{v_{\mathrm{P}}}}}{{{\mathrm{d}}\theta }} （ {\boldsymbol{t}} \cos \theta－ {\boldsymbol{n}} \sin \theta ） \text{，} $$ (17) 式中vg和vP（θ）分别表示群速度和相速度。对于TTI介质，相角由波前梯度方向经过坐标旋转变换后求得。${\boldsymbol{ t}}$可用相角以及${\boldsymbol{n}}$表示，即${\boldsymbol{t}}＝{1}/{{\sin \theta }}\left({{\nabla \tau }}/{{\left| {\nabla \tau } \right|}} － {\boldsymbol{n}} \cos \theta \right)$，因此式（17）中的群速度矢量也可表示为
$$ \begin{split} \boldsymbol{v}_{\mathrm{g}}＝ & \boldsymbol{v}_{\mathrm{P}} （ \theta ） ＋\frac{\mathrm{d}v_{\mathrm{P}}}{\mathrm{d}\theta} （ \boldsymbol{t}\cos\theta－\boldsymbol{n}\sin\theta ） ＝v_{\mathrm{P}}\frac{\nabla\tau}{\left|\nabla\tau\right|}＋\frac{\mathrm{d}v_{\mathrm{P}}}{\mathrm{d}\theta} （ \boldsymbol{t}\cos\theta－\boldsymbol{n}\sin\theta ） ＝v_{\mathrm{P}}^2 （ \theta ） \nabla\tau＋ \\ &\frac{\mathrm{d}v_{\mathrm{P}}}{\mathrm{d}\theta}\left(\frac{\cos\theta}{\sin\theta}\frac{\nabla\tau}{\left|\nabla\tau\right|}－\boldsymbol{n}\frac{\cos^2\theta}{\sin\theta}－\boldsymbol{n}\sin\theta\right)＝v_{\mathrm{P}}^2 （ \theta ） \nabla\tau＋\frac{\mathrm{d}v_{\mathrm{P}}}{\mathrm{d}\theta}\left(\frac{\nabla\tau}{|\nabla\tau|}\cos\theta－\boldsymbol{n}\right)\frac{1}{\sin\theta} \text{，} \end{split} $$ (18) 式中$ {\boldsymbol{n}}＝ {{\boldsymbol{z}}} \cos \alpha ＋ {\boldsymbol{x}}\sin \alpha $；④ 射线路径由接收点开始，逆着群速度矢量场方向可以一直追踪到源点，其方向由邻近节点的群速度方向插值求出。走时由每个网格的射线长度除以群速度后求和获得。
3. 数值算例
基于均匀速度TTI模型、速度倾斜渐变的TTI模型、垂直横向各向同性（vertical transverse isotropy，缩写为VTI ）模型、Marmousi模型以及近地表模型计算了五种模型的走时和射线路径，以此计算结果分析本文所提的方法在计算效率和计算精度上的性能。
3.1 均匀速度TTI模型
为了测试本文方法的准确性，设计一个均匀TTI模型进行试算。速度为2 km/s，模型横向、纵向尺寸均为0—2 km，网格间距均为10 m，倾角θ＝45°，各向异性参数ε＝0.1，δ＝0.1，震源位于模型中央（1 km，1 km）。图3是均匀模型走时及其绝对误差的对比结果，其中：图3a为PVP-FMM计算的走时与精确解的对比，可以看出使用本文方法计算得到的走时与精确解基本吻合；图3b为走时绝对误差，从中可以看出，本文方法的最大走时误差约为1.5 ms，初步验证了本文方法的正确性。
![]()
图 3 均匀模型走时及误差对比（a） 走时对比；（b） 走时绝对误差对比Figure 3. Comparison of travel time and its error for homogeneous model（a） Comparison of travel time；（b） Comparison of absolute error of travel time3.2 速度倾斜渐变的TTI模型
为了测试本文方法在TTI模型上的性能，设计了一个地层法向速度倾斜渐变的TTI模型。模型纵横向尺寸均为2 km，网格间距均为10 m，各向异性参数ε＝0.16，δ＝0.04，θ＝45°。如图4a所示，其模型速度随深度的增加而增大，初始速度为800 m/s，深度每增加1 m速度增加2 m/s。图4b为PVP-FMM与无相角预测方法计算的波前时间之差。其中，无相角预测方法直接利用波前已知的与邻近节点相速度近似的FMM，将其缩写为No-PVP-FMM。
![]()
图 4 TTI模型走时之差对比（a） 速度模型；（b） PVP-FMM与No-PVP-FMM走时之差Figure 4. Comparison of travel time difference for TTI model（a） Velocity model；（b） Travel time difference between PVP-FMM and No-PVP-FMM从图4b走时绝对误差图可以看出，在其水平位置（45°方向）的误差基本为零，这是因为相角为零时，传播误差与ε和δ无关，误差最小，但是水平位置（45°方向）到90°方向误差逐渐增大，是由相角变化所致，90°到地层法向方向的误差则是由速度变化差异较大所引起。
图5为炮点位于TTI模型中心时，利用本文方法计算所得的射线追踪过程。首先利用PVP-FMM方法计算得到所有网格点的波前时间场（图5a）；其次，利用式（16）得到其相速度矢量场（图5b）；然后根据式（18）获得群速度矢量场（图5c）；最后，射线路径由接收点开始，逆着群速度矢量场方向一直追踪到炮点，如图5d所示。由图5可知，射线、等时线的密度都与模型参数的变化相符合。
![]()
图 5 利用本文方法PVP-FMM计算所得的射线追踪过程（炮点位于模型中心）（a） 波前时间场；（b） 相速度矢量场；（c） 群速度矢量场；（d） 射线路径Figure 5. The ray tracing process calculated by PVP-FMM in this paper （the shot point is located in the center of the model）（a） Wave front time field；（b） Phase velocity vector field；（c） Group velocity vector field； （d） Ray path3.3 VTI模型
SEG Advanced modeling （SEAM）盐丘模型有多个数据集，这里我们采用其模型的一部分来验证本文方法的效果。模型参数为：模型横向、纵向坐标范围均为0—10 km，将其划分为201×201的网格，网格间距为50 m，各向异性介质参数ε＝0.2，δ＝0.12，模型的炮点位置在（3 km，0）处，模型速度如图6a所示。对该模型分别使用PVP-FMM和No-PVP-FMM方法进行走时计算和射线追踪，图6b为PVP-FMM的结果，图7为两种方法的结果与参考值的对比。图7中红色为SEAM公司所公布的参考走时解（Fehler，Keliher，2011），黑色为PVP-FMM方法，蓝色为No-PVP-FMM方法，从图中可以看出，两种算法之间差异明显，尤其在速度变化较大的地方本文方法与参考解吻合得很好，而No-PVP-FMM方法却在一些地方偏差较大。在走时误差图（图8）中也可以看出，No-PVP-FMM方法计算的解析解与参考解的最大误差约为0.14 s，本文方法计算结果与参考解的最大误差为0.018 s，可以发现本文方法比No-PVP-FMM方法的误差更小，得到的走时更加精确。
![]()
图 6 VTI模型追踪结果（a） 速度模型；（b） 波前等时线及射线追踪Figure 6. Ray tracing results of VTI model（a） Velocity model；（b） Wavefront isochron and ray path3.4 Marmousi模型
为了验证本文方法处理复杂模型的性能，对Marmousi模型进行了测试。该模型大小为2 km×2 km，网格为101×101，网格间距为20 m，ε＝0.16，δ＝0.04，速度如图9a所示。图9b是源点位于模型中央时本文方法的计算结果。
![]()
图 8 走时绝对误差图Figure 8. Comparison of absolute error of travel time（a） No-PVP-FMM；（b） PVP-FMM![]()
图 9 Marmousi模型追踪结果（a） 速度模型；（b） 波前等时线及射线追踪Figure 9. Ray tracing results of Marmousi model（a） Velocity model；（b） Wavefront isochron and ray path将震源位置移动到地面（1 km，0）后，分别使用PVP-FMM和No-PVP-FMM方法对模型进行计算。图10显示了两种方法获得的等时线与参考解之间的对比及误差图。图10a为全局图，图10b为局部放大图。从全局图中可知，两种方法的等时线结果几乎重合，但在速度变化快的地方可以看出明显的差异；在图10b的局部放大图中也可看出，No-PVP-FMM在发生速度变化大的区域内出现较大的偏差。走时绝对误差如图11所示，从图中可以看出，PVP-FMM算法明显比No-PVP-FMM的误差更小，可以说明PVP-FMM方法能明显改善速度变化剧烈区域的计算精度，针对复杂模型有更好的计算精度。
![]()
图 10 走时等值线对比图（a） 走时对比；（b） 局部放大图Figure 10. Comparison map of travel time contour（a） Comparison of travel time；（b） Local enlarged image![]()
图 11 走时绝对误差图Figure 11. Comparison of absolute error of travel time（a） No-PVP-FMM；（b） PVP-FMM3.5 近地表数值模拟
最后，为了测试本文算法在近地表各向异性模型中的实用性，根据实际地表状况，本文设计了一个地表起伏、地下结构复杂的近地表TTI模型。速度模型和地层倾角模型如图12a和12b所示，近地表各向异性参数值为ε＝0.16，δ＝0.04，随着深度的变化，参数也随之变化，如图12c和12d所示，观测系统设计了913个震源，每个震源有215个接收点。
![]()
图 12 近地表模型（a） 速度模型；（b） 倾角模型；（c） δ模型；（d） ε模型Figure 12. Near-surface models（a） Velocity；（b） Dip model；（c） δ model；（d） ε model为了与各向同性的FMM进行比较，对该模型分别使用各向同性和各向异性的波场进行有限差分法数值模拟，图13显示的是模拟各向异性介质的一个单炮记录。分别使用各向同性FMM、各向异性No-PVP-FMM和PVP-FMM方法对模型进行射线追踪，图14显示了各向异性模型任意一炮的计算初至与拾取初至的对比图，从图中可以看出，No-PVP-FMM和PVP-FMM的计算初至与从数值模拟中拾取的初至，除了在接近7 km处存在一个异常点之外，几乎全部重合。
![]()
图 14 各向异性模型的计算初至与拾取初至（单炮）Figure 14. First arrival and picked-first arrival calculated by the anisotropic model （single shot）图15给出了模拟波场初至拾取的结果。图中上部显示的是多个震源的初至时间随偏移距变化的曲线图。横坐标为测线距离，纵坐标为走时。下部是初至时间的炮点和接收点位置。从图中可观测到，与各向同性的初至时间相比，各向异性更小，二者有明显差异。本文方法计算结果与波动方程的有限差分模拟结果具有几乎相同的初至时间。个别异常的走时值源于初至拾取时产生的误差。
![]()
图 15 各向同性（a）和各向异性（b）数值模拟记录拾取的初至时间Figure 15. First arrival time picked up from isotropic （a） and anisotropic （b） numerical simulation records为了更深入地分析方法的计算误差，本文对比了所有地震道的计算结果与拾取时间的误差，结果如图16所示，可以看出，各向同性与各向异性的误差平均值和方差具有相同的量级，误差分布也相近。
![]()
图 16 射线追踪拟合的误差（913炮）（a） 各向同性FMM射线追踪拟合；（b） 各向异性No-PVP-FMM射线追踪拟合；（c） PVP-FMM射线追踪拟合Figure 16. Error of ray tracing fitting （913 shots）（a） Isotropic FMM ray tracing fitting；（b） Anisotropic No-PVP-FMM ray tracing fitting；（c） PVP-FMM ray tracing fitting为了进一步量化分析方法的性能，对913炮的各向同性FMM、各向异性No-PVP-FMM和PVP-FMM的计算误差和计算时间进行统计分析，表1列出了三种方法的计算效率和统计误差。从表中可以看出，三种方法具有相同量级的计算误差。误差均值都很小，各向同性FMM、各向异性No-PVP-FMM和PVP-FMM均方根误差分别为8.6 ms，9.3 ms和9.0 ms，No-PVP-FMM比各向同性方法的增加了0.7 ms，而PVP-FMM仅增加了0.4 ms，比No-PVP-FMM减少了43%。使用2.4 GHz CPU的计算机进行射线追踪，三种方法913炮的总耗时分别为326 s，359 s和372 s。由于No-PVP-FMM增加了群速度计算环节，耗时比各向同性FMM增加10%，而PVP-FMM又添加了相速度预测环节，耗时增加了14%。本文方法与各向同性的FMM具有相同量级计算效率，对于近地表模型同样具有良好的适应性和稳定性。
表 1 三种方法的计算精度和计算效率Table 1. Calculation accuracy and efficiency of the three methods算法 平均误差/ms 均方根误差/ms CPU计算时间/s 各向同性FMM −1.12 8.6 326 No-PVP-FMM −0.03 9.3 359 PVP-FMM −0.16 9.0 372 由以上数值算例可知，本文方法对于不同速度结构模型具有良好的适应性。从Marmousi以及VTI模型中可看出，PVP-FMM方法误差更小，具有更高的走时计算精度。在计算效率方面，虽然PVP-FMM在No-PVP-FMM计算群速度环节的基础上增加了相速度预测环节，但耗时也仅比各向同性FMM增加了14%。此外，本文所提出的方法适用于非均匀TTI介质，因为在均匀介质中，各节点的速度相同，本文所给出的相速度预测公式不起作用。
4. 结论
本文给出的TTI介质的相速度预测公式和PVP-FMM各向异性射线追踪方法，继承了各向同性FMM快速稳定的特点，具有与各向同性FMM相同量级的计算精度。通过理论模型、VTI模型和Marmousi模型测试，验证了本文方法具有小的计算误差和良好的稳定性，能够适用于复杂构造的射线追踪。通过复杂近地表多源模型测试，验证了本文方法的计算结果与有限差分模拟的波场时间能够很好地拟合，在提高FMM走时计算精度的同时，耗时也仅比各向同性FMM增加了14%。
-
![]()
图 1 超强台风 “利奇马” 的最佳路径和强度演变及东南沿海地区钻孔体应变台和气象站位置
Figure 1. The best track and intensity （colored circles） of super typhoon Lekima （from 14:00 BJT on 4 August 2019 to 11:00 BJT on 13 August 2019，equally spaced at 1 h or 3 h，respectively） marked with time as well as locations of borehole dilatometer stations （black triangles）and meteorological stations （green triangles） in southeastern coastal area of China
![]()
图 2 5个钻孔体应变台记录的 “利奇马” 低频扰动曲线
Figure 2. The low-frequency signatures of super typhoon Lekima recorded by the borehole dilatometers at five stations during 1−17 August 2019
Time series after linear curves are removed from the original one-minute-sampled volumetric strain data （black lines），and red lines show the trends of long-period changes in volumetric strain. The vertical left and right dashed lines mark the timing of Lekima's initiation and first landfall，respectively. The blue circle indicates the time when the typhoon was the nearest to the borehole dilatometer station.
![]()
图 3 东阳台（a）、溧阳台（b）和杭州气象站在 “利奇马” 过境期间的观测数据
图中第一行为钻孔体应变和气压记录及趋势，第二行为钻孔体应变和气压的变化率，第三行为日降水量
Figure 3. Records of the volumetric strain and barometric pressure at the stations Dongyang （a） and Liyang （b） as well as daily rainfall data at Hangzhou meteorological station under the passage of super typhoon Lekima during August 1−17，2019
The upper panels show the traces of the volumetric strain （black） and barometric pressure （gray），where red and green lines show the trends，respectively；the middle panels represent the variation rates of the volumetric strain （red） and barometric pressure （green）；the lower panels represent the daily rainfall
![]()
图 5 青岛台和青岛气象站在 “利奇马” 过境期间的观测数据
（a） 钻孔体应变和气压记录及趋势；（b） 钻孔体应变和气压的变化率；（c） 日降水量
Figure 5. Records of the volumetric strain and barometric pressure at Qingdao station and daily rainfall data at Qingdao meteorological station under the passage of super typhoon Lekima during August 1−17，2019
（a） Traces of the volumetric strain （black） and barometric pressure （gray），red and green lines show the trends；（b） The variation rates of the volumetric strain （red） and barometric pressure （green）；（c） The daily rainfall
![]()
图 6 非台风情况下强降水事件对东阳台观测数据的影响
（a） 钻孔体应变和气压记录及趋势；（b） 钻孔体应变和气压的变化速率；（c） 日降水量
Figure 6. Records of the volumetric strain and barometric pressure at Dongyang station and daily rainfall data at Hangzhou meteorological station under the heavy rainfall induced by non-typhoon weather during July 1−20，2019
（a） Traces of the volumetric strain （black） and barometric pressure （gray），red and green lines show the trends；（b） The variation rates of the volumetric strain （red）and barometric pressure （green）；（c） The daily rainfall
![]()
图 4 常熟台（a），南通台（b）和上海气象站在 “利奇马” 过境期间的观测数据
图中第一行为钻孔体应变和气压记录及趋势，第二行为钻孔体应变和气压的变化率，第三行为日降水量
Figure 4. Records of the volumetric strain and barometric pressure at the stations Changshu （a） and Nantong （b） as well as and daily rainfall data at Shanghai meteorological station under the passage of super typhoon Lekima during August 1−17，2019
The upper panels show the traces of the volumetric strain （black） and barometric pressure （gray），where red and green lines show the trends，respectively；the middle panels represent the variation rates of the volumetric strain （red） and barometric pressure （green）；the lower panels represent the daily rainfall
表 1 5个钻孔体应变台站及其周邻地面气象站的概况
Table 1 General information for the five borehole dilatometer sites and three neighboring meteorological stations
台站 钻孔深度/m 钻孔围岩岩性 距海岸线距离/km 邻近的气象站 距气象站距离/km 东阳 67 泥质粉砂岩 120 杭州 110 溧阳 60 安山岩 240 杭州 140 常熟 62 石英砂岩 110 上海 80 南通 94 石英砂岩 50 上海 85 青岛 56 花岗岩 0.2 青岛 6 
下载: 导出CSV
表 2 5个钻孔体应变台站对 “利奇马” 响应的特征参数
Table 2 The response patterns and magnitudes to Lekima for the five borehole dilatometer stations
台站 响应距离/km 台风
临近时
的强度历
时
/h气压
最大变幅
/hPa体应变
最大变幅
/10−9气压
变化率
/（hPa·min−1）体应变
变化率
/（10−9 min−1）累计
降水量
/mm气压影响
系数
/（10−9 hPa−1）开始响应 台风登陆 台风临近 结束响应 东阳 770 150 35 900 强热带风暴 89 −20.1 −52.8 −0.016—0.014 −0.040—0.038 167 2.6 溧阳 830 390 88 650 热带风暴 96 −17.2 −36.8 −0.012—0.009 −0.028—0.027 167 2.1 常熟 760 380 30 640 热带风暴 103 −18.2 −112.1 −0.012—0.011 −0.069—0.076 141 6.2 南通 780 410 27 - 热带风暴 120 −20.7 −76.7 −0.011—0.010 −0.037—0.034 142 3.7 青岛 980 870 26 - 热带风暴 107 −22.8 −72.1 −0.017—0.009 −0.030—0.041 54 3.2 注：“-”表示此时刻台风已停止编号，无法获取该时刻台风的具体位置。
下载: 导出CSV
表 3 各台站的台基岩石力学参数和气压波动幅值及相应的理论体应变
Table 3 The modeled volumetric strain induced by observed atmospheric loading based on corresponding rock mechanical parameters for the five borehole dilatometer stations
台站 钻孔围岩力学参数 气压波动/hPa 实测体应变/10−9 理论体应变/10−9 弹性模量/GPa 泊松比 东阳 20 0.25 −20.1 −52.8 −50.3 溧阳 40 0.25 −17.2 −36.8 −21.5 常熟 10 0.25 −18.2 −112.1 −91.0 南通 15 0.25 −20.7 −76.7 −69.0 青岛 50 0.25 −22.8 −72.1 −22.8
下载: 导出CSV
-
陈孔沫. 1981. 台风气压场和风场模式[J]. 海洋学报,<bold>3</bold>(1):44–56. Chen K M. 1981. The typhoon pressure field and wind field models[J]. <italic>Acta Oceanologica Sinica</italic>,<bold>3</bold>(1):44–56 (in Chinese).
冯志军,赵金花,荆强,陈安范,殷海涛,李杰. 2009. 青岛钻孔体应变短时畸变分析[J]. 山西地震,(4):42–47. Feng Z J,Zhao J H ,Jing Q,Chen A F,Yin H T,Li J. 2009. Study on the short time distortion anomaly of borehole volume strain in Qingdao[J]. <italic>Earthquake Research in Shanxi</italic>,(4):42–47 (in Chinese).
梁溯安. 2012. 东阳江（流域）水污染物总量控制研究[D]. 杭州: 浙江大学: 29. Liang S A. 2012. Total Amount Control of Water Pollution in Dongyangjiang River[D]. Hangzhou: Zhejiang University: 29 （in Chinese）.
彭剑文,曾飞涛,李长洪,苗胜军. 2017. 石英砂岩力学特性及各向异性试验研究[J]. 岩土力学,<bold>38</bold>(增刊1):103–112. Peng J W,Zeng F T,Li C H,Miao S J. 2017. Experimental study of anisotropy and mechanical property of quartz sandstone[J]. <italic>Rock and Soil Mechanics</italic>,<bold>38</bold>(S1):103–112 (in Chinese).
邱泽华,唐磊,张宝红,宋茉. 2012. 用小波—超限率分析提取宁陕台汶川地震体应变异常[J]. 地球物理学报,<bold>55</bold>(2):538–546. Qiu Z H,Tang L,Zhang B H,Song M. 2012. Extracting anomaly of the Wenchuan earthquake from the dilatometer recording at NSH by means of wavelet-overrun rate analysis[J]. <italic>Chinese Journal of Geophysics</italic>,<bold>55</bold>(2):538–546 (in Chinese).
邱泽华. 2017. 钻孔应变观测理论与应用[M]. 北京: 地震出版社: 1–407. Qiu Z H. 2017. The Observations of Borehole Strainmeters: Theory and Applications[M]. Beijing: Seismological Press: 1–407 （in Chinese）.
苏恺之, 李秀环, 张钧, 马相波, 李海亮. 2002. TJ-2型体应变仪的研制[G]//地壳构造与地壳应力文集. 北京: 中国地震局地壳应力研究所: 113−121. Su K Z, Li X H, Zhang J, Ma X B, Li H L. 2002. Manufacture of TJ-2 volume strain meter[G]//Bulletin of the Institute of Crustal Dynamics. Beijing: Institute of Crustal Dynamics, China Earthquake Administration: 113−121 （in Chinese）.
袁金南,林爱兰,刘春霞. 2008. 60年来西北太平洋上不同强度热带气旋的变化特征[J]. 气象学报,<bold>66</bold>(2):213–223. Yuan J N,Lin A L,Liu C X. 2008. Change characters of tropical cyclones with different intensities over the western North Pacific during the last 60 years[J]. <italic>Acta Meteorologica Sinica</italic>,<bold>66</bold>(2):213–223 (in Chinese).
袁媛,方国庆,尹京苑. 2017. 佘山台四分量钻孔应变仪记录的台风扰动特征分析[J]. 地震学报,<bold>39</bold>(5):725–737. doi: 10.11939/jass.2017.05.008 Yuan Y,Fang G Q,Yin J Y. 2017. Signatures of typhoons on strain records of four component borehole strainmeter at Sheshan station[J]. <italic>Acta Seismologica Sinica</italic>,<bold>39</bold>(5):725–737 (in Chinese).
张凌空,牛安福. 2008. 中国钻孔体应变仪同震变化观测结果[J]. 国际地震动态,(11):120. doi: 10.3969/j.issn.0253-4975.2008.11.121 Zhang L K,Niu A F. 2008. Co-seismic strain changes recorded by dilatometers in China[J]. <italic>Recent Developments in World Seismology</italic>,(11):120 (in Chinese).
张凌空,牛安福. 2019. 周期气压波对地壳岩石应变测量影响的理论解[J]. 地球物理学进展,<bold>34</bold>(4):1366–1370. doi: 10.6038/pg2019CC0244 Zhang L K,Niu A F. 2019. Theoretical solution of periodic pressure wave effect on crustal strain[J]. <italic>Progress in Geophysics</italic>,<bold>34</bold>(4):1366–1370 (in Chinese).
中央气象台台风网. 2019. 1909利奇马Lekima[EB/OL]. [2019−09−02]. http://typhoon.nmc.cn/web.html. Typhoon and Marine Weather Monitoring and Warning, National Meteorological Centre. 2019. Super typhoon Lekima （1909）[EB/OL]. [2019−09−02]. http://typhoon.nmc.cn/web.html （in Chinese）.
周龙寿,邱泽华. 2008. 地壳应变场对气压短周期变化的响应[J]. 地球物理学进展,<bold>23</bold>(6):1717–1726. Zhou L S,Qiu Z H. 2008. The response of crustal strain field to short-period atmospheric pressure variation[J]. <italic>Progress in Geophysics</italic>,<bold>23</bold>(6):1717–1726 (in Chinese).
檜皮久義,佐藤馨,二瓶信一,福留篤男,竹内新,古屋逸夫. 1983. 埋込式体積歪計の気圧補正[J]. 験震時報,<bold>47</bold>:91–111. Hikawa H,Sato K,Nihei S,Fukudome A,Takeuchi H,Furuya I. 1983. Correction due to atmospheric pressure changes of data of the borehole volume strainmeter[J]. <italic>Quarterly Journal of Seismology</italic>,<bold>47</bold>:91–111 (in Japanese).
木村一洋,露木貴裕,菅沼一成,長谷川浩,見須裕美,藤田健一. 2015. タンクモデルによる体積ひずみ計データの降水補正について[J]. 験震時報,<bold>78</bold>:93–158. Kimura K,Tsuyuki T,Suganuma I,Hasegawa H,Misu H,Fujita K. 2015. Rainfall correction of volumetric stainmeter data by tank models[J]. <italic>Quarterly Journal of Seismology</italic>,<bold>78</bold>:93–158 (in Japanese).
上垣内修. 1987. 体積歪,傾斜データに対する気圧の影響の補正に関する物理的考察[J]. 験震時報,<bold>50</bold>:41–49. Kamigaichi O. 1987. Physical considerations on the correction methods of volumetric strain and tilt data for the effects of atmospheric pressure chang[J]. <italic>Quarterly Journal of Seismology</italic>,<bold>50</bold>:41–49 (in Japanese).
Asai Y,Ishii H,Aoki H. 2009. Comparison of tidal strain changes observed at the borehole array observation system with in situ rock properties in the Tono region,central Japan[J]. <italic>J Geodyn</italic>,<bold>48</bold>(3/4/5):292–298. doi: 10.1016/j.jog.2009.09.024
Barbour A J,Crowell B W. 2017. Dynamic strains for earthquake source characterization[J]. <italic>Seism Res Lett</italic>,<bold>88</bold>(2A):354–370. doi: 10.1785/0220160155
Bonaccorso A,Linde A,Currenti G,Sacks S,Sicali A. 2016. The borehole dilatometer network of Mount Etna:A powerful tool to detect and infer volcano dynamics[J]. <italic>J Geophys Res</italic>:<italic>Solid Earth</italic>,<bold>121</bold>(6):4655–4669. doi: 10.1002/2016JB012914
Evertson D W. 1977. Borehole Strainmeters for Seismology[R]. Austin, Texas: Applied Research Lab, University of Texas: No. ARL-TR-77-62.
Hsu Y J,Chang Y S,Liu C C,Lee H M,Linde A T,Sacks S I,Kitagawa G,Chen Y G. 2015. Revisiting borehole strain,typhoons,and slow earthquakes using quantitative estimates of precipitation-induced strain changes[J]. <italic>J Geophys Res</italic>:<italic>Solid Earth</italic>,<bold>120</bold>:4556–4571. doi: 10.1002/2014JB011807
Huang N E,Shen Z,Long S R,Wu M C,Shih H H,Zheng Q,Yen N C,Tung C C,Liu H H. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. <italic>Proc R Soc Lond A</italic>,<bold>454</bold>:903–995. doi: 10.1098/rspa.1998.0193
Linde A T,Gladwin M T,Johnston M J S,Gwyther R L,Bilham R G. 1996. A slow earthquake sequence on the San Andreas fault[J]. <italic>Nature</italic>,<bold>383</bold>:65–68. doi: 10.1038/383065a0
Liu C C,Linde A T,Sacks I S. 2009. Slow earthquakes triggered by typhoons[J]. <italic>Nature</italic>,<bold>459</bold>(7248):833–836. doi: 10.1038/nature08042
Mouyen M,Canitano A,Chao B F,Hsu Y J,Steer P,Longuevergne L,Boy J P. 2017. Typhoon-induced ground deformation[J]. <italic>Geophys Res Lett</italic>,<bold>44</bold>(21):11004–11011. doi: 10.1002/2017GL075615
Matsuura T,Yumoto M,Iizuka S. 2003. A mechanism of interdecadal variability of tropical cyclone activity over the western North Pacific[J]. <italic>Clim Dyn</italic>,<bold>21</bold>:105–117. doi: 10.1007/s00382-003-0327-3
NOAA. 2019. National centers for environmental information[EB/OL]. [2019−09−02]. https://www.ncdc.noaa.gov/cdo-web/datatools/findstation.
Peduzzi P,Chatenoux B,Dao H,Bono A D,Herold C,Kossin J,Mouton F,Nordbeck O. 2012. Global trends in tropical cyclone risk[J]. <italic>Nat Clim Change</italic>,<bold>2</bold>(2):289–294.
Peng H,Kitagawa G,Takanami T,Matsumoto N. 2014. State-space modeling for seismic signal analysis[J]. <italic>Appl Math Model</italic>,<bold>38</bold>(2):738–746. doi: 10.1016/j.apm.2013.07.008
Roeloffs E A. 2006. Evidence for aseismic deformation rate changes prior to earthquakes[J]. <italic>Annu Rev Earth Planet Sci</italic>,<bold>34</bold>:591–627. doi: 10.1146/annurev.earth.34.031405.124947
Sacks I S,Suyehiro S,Evertson D W. 1971. Sacks-Evertson strainmeter,its installation in Japan and some preliminary results concerning strain steps[J]. <italic>Proc Jpn Acad</italic>,<bold>47</bold>(9):707–712. doi: 10.2183/pjab1945.47.707
Sacks I S,Suyehiro S,Linde A T,Snoke J A. 1978. Slow earthquakes and stress redistribution[J]. <italic>Nature</italic>,<bold>275</bold>:599–602. doi: 10.1038/275599a0
Takanami T,Linde A T,Sacks S I,Kitagawa G,Peng H. 2013. Modeling of the post-seismic slip of the 2003 Tokachi-Oki earthquake <italic>M</italic>8 off Hokkaido:Constraints from volumetric strain[J]. <italic>Earth Planets Space</italic>,<bold>65</bold>:731–738. doi: 10.5047/eps.2012.12.003
Wu L G,Wang B,Geng S Q. 2005. Growing typhoon influence on east Asia[J]. <italic>Geophys Res Lett</italic>,<bold>32</bold>:L18703.
Xu X D,Peng S Q,Yang X J,Xu H X,Tong D Q,Wang D X,Guo Y D,Chan J C L,Chen L S,Yu W,Li Y N,Lai Z J,Zhang S J. 2013. Does warmer China land attract more super typhoons?[J]. <italic>Sci Rep</italic>,<bold>3</bold>:1522. doi: 10.1038/srep01522
-
期刊类型引用(4)
1. 韩伟，蒋兴超，王建强，李玉宏，郭望，张云鹏，陈高潮. 汉中地区构造演化及寒武系页岩气形成地质条件研究. 地质学报. 2024(06): 1829-1839 . 
百度学术
2. 孟文，田涛，孙东生，杨跃辉，李冉，陈群策. 基于原位地应力测试及流变模型的深部泥页岩储层地应力状态研究. 地质力学学报. 2022(04): 537-549 . 百度学术
3. 程先琼，蒋科植. 基于深度降噪自编码神经网络的中国大陆地壳厚度反演. 地震学报. 2021(01): 34-47+136 . 本站查看
4. 王婷，延军平，李双双，万佳，张玉凤. 帕米尔高原Mw≥6.6级地震时间韵律特征. 高原地震. 2020(04): 6-16 . 百度学术
其他类型引用(2)
计量
- 文章访问数: 873
- HTML全文浏览量: 571
- PDF下载量: 61
- 被引次数: 6
