基于移动窗解卷积法的2011年日本东北 MW9.0大地震场地非线性时变识别

刘启方, 陈长龙

刘启方,陈长龙. 2022. 基于移动窗解卷积法的2011年日本东北 MW9.0大地震场地非线性时变识别. 地震学报,44(1):96−110. DOI: 10.11939/jass.20210080
引用本文: 刘启方,陈长龙. 2022. 基于移动窗解卷积法的2011年日本东北 MW9.0大地震场地非线性时变识别. 地震学报,44(1):96−110. DOI: 10.11939/jass.20210080
Liu Q F,Chen C L. 2022. Identification of the temporal changes of site nonlinearity during 2011 MW9.0 Tohoku earthquake by moving time window deconvolution method. Acta Seismologica Sinica44(1):96−110. DOI: 10.11939/jass.20210080
Citation: Liu Q F,Chen C L. 2022. Identification of the temporal changes of site nonlinearity during 2011 MW9.0 Tohoku earthquake by moving time window deconvolution method. Acta Seismologica Sinica44(1):96−110. DOI: 10.11939/jass.20210080

基于移动窗解卷积法的2011年日本东北 MW9.0大地震场地非线性时变识别

基金项目: 国家自然科学基金(51978434)资助
详细信息
    通讯作者:

    刘启方,博士,教授,主要从事地震工程研究,e-mail:qifang_liu@126.com

  • 中图分类号: P315.9

Identification of the temporal changes of site nonlinearity during 2011 MW9.0 Tohoku earthquake by moving time window deconvolution method

  • 摘要: 利用模拟记录和2011年日本东北MW9.0大地震观测记录分析了基于移动窗解卷积法识别场地非线性时变特征的可行性,并与移动窗谱比法的结果进行了对比分析。研究表明:基于移动窗解卷积法可以较好地揭示场地非线性随地震动水平的变化过程,识别非线性发生的阈值、非线性变化程度及强震动后的恢复程度;与移动窗谱比法相比,移动窗解卷积法更容易获得较为稳定的土体非线性时变过程,但对于存在强阻抗比的浅表层土体,移动窗谱比法可以获得更准确的非线性程度。对2011年日本东北MW9.0大地震中8个KiK-net台站进行了非线性时变分析,结果表明;两种方法识别的非线性阈值较接近,约在40—100 cm/s2之间,且与场地vS30没有明显的相关性;在峰值加速度PGA较低的IBRH20台站,非线性引起的波速下降较小(3%)且震后几乎完全恢复;PGA 处于386—822 cm/s2之间的其余7个台站,场地等效剪切波速下降13%—37%,产生了显著的场地非线性,且震后未完全恢复;PGA大于380 cm/s2时,非线性所导致的场地波速下降、恢复与PGA无明显相关性。
    Abstract: In this paper, the feasibility of identifying nonlinear temporal changes of sites by moving time window deconvolution method is analyzed based on the simulation records and observation records from the 2011 MW9.0 Tohoku earthquake, and the results are compared with those by the moving time window spectral ratio method. The results show that the moving time window deconvolution method can reveal the nonlinear change process of the site with ground motion level. Based on this method, the threshold of nonlinearity, the degree of nonlinearity change and the recovery degree after strong ground motion process can be identified. Compared with the moving time window spectrum ratio method, the moving time window deconvolution method can obtain more stable nonlinear temporal changes process of soil, but for the shallow surface soil with strong impedance ratio, the moving time window spectrum ratio method can obtain more accurate nonlinear degree results. The nonlinear temporal changes analysis of eight KiK-net stations during the 2011 MW9.0 Tohoku earthquake shows that the nonlinear thresholds identified by the two methods are close to each other, ranging from 40 cm/s2 to 100 cm/s2, and there is no obvious correlation with site vS30. At IBRH20 station with low PGA, the decrease of wave velocity caused by nonlinearity was small (3%) and the recovery was almost complete after the earthquake. At the other seven stations with high PGA (range of 386—822 cm/s2), the site equivalent shear wave velocity decreased by 13%—37%, resulting in significant site nonlinearity and not fully recovered after the earthquake. As PGA is larger than 380 cm/s2, there is no significant correlation between PGA and the decrease/recovery of wave velocity caused by nonlinearity.
  • 通过理论分析和数值计算得到的理论地震图是地震学家们了解震源机制和地球介质性质的重要基础。自Lamb (1904)开展了自由界面对于施加在其上的力的响应的经典工作以来,经过一个多世纪的发展,特别是伴随着近50年来计算技术的迅猛发展,合成地震图技术已经取得了长足的进步。20世纪50年代以前,地震学家和数学家们研究介质响应多是采取渐近分析手段得到近似解(Lapwood,1948Pekeris,1955),而从50年代计算机引入到地震学领域中以来,各种算法不断涌现。根据解析程度的大小,大致可分成半解析-半数值方法和数值方法两类。对于半解析-半数值方法,主要有基于Cagniard-de Hoop技术的广义射线方法(Helmberger,19731974Johnson,1974)和基于矩阵运算的反射率方法(Fuchs,Müller,1971Kennett,Kerry,1979Luco,Apsel,1983Chen,1999)。对于数值方法,各种不同的计算方式层出不穷,比如有限差分法(Madariaga,1976Levander,1988Graves,1996Zhang,Chen,2006)、谱元法(Priolo,Seriani,1991Komatitsch et al,2001Capdeville et al,2003Komatitsch,Tromp,2003Tromp,Komatitsch,2008)等。因数值方法对介质的适用性更广,越来越受到地震学家们的重视。但是在很多情况下,如果介质可用一维结构近似,半解析-半数值的方法,特别是能够合成完整地震波场的反射率方法,由于其计算效率方面的优势仍然会得到广泛的应用。

    在反射率方法中,基于离散波数法的广义反透射系数法(Apsel,Luco,1983Yao,Harkri-der,1983Chen,1999)是目前应用最广泛的方法。但是,当震源和接收点的深度相近或相同时,由于被积函数收敛得非常慢,仅仅用离散波数法而不采取其它特殊技术时计算效率很低。Apsel和Luco (1983)Hisada (19941995)采用了半解析的方法,通过将积分拆成两部分的方式巧妙地解决了被积函数的慢收敛问题,但缺点是数学处理非常繁琐,且在实际的数值计算中需要特殊的处理技术。基于重复平均法的原则,张海明等(2001)提出了一种计算浅源情况下的数值积分方法,即峰谷平均法。与半解析方法相比,该方法不仅数学处理上非常简单,而且易于数值实现。

    运用峰谷平均法,需要在波数k大于阈值kc之后构造由积分值变化曲线的波峰和波谷组成的振荡序列,而这个kc在之前的研究(张海明等,2001)中是由经验公式给出的。但是,一方面,经验公式缺乏理论依据;另一方面,在有些情况下,根据经验公式确定的kc比实际的积分限阈值大得多,导致计算效率不高。因此有必要思考新的求解方法。在本文的研究中,我们发现波数积分的积分限阈值kc与反透射系数行列式的零点有关,拟通过求解频率域中瑞雷函数的零点,找到可以精确地确定积分限阈值的方案,以期通过精确确定的积分阈值能够以更高的效率计算理论地震图。

    根据广义反透射系数法,分层半空间中任意点源产生的位移在频率域中可表示为如下类型的波数积分(张海明等,2001):

    $$ Y\left( {r{\text{,}}z{\text{,}}{{{z}}_{\rm{s}}}{\text{,}}\omega } \right) {\text{=}} \mathop \int \nolimits_0^{ + \infty } F\left( {r{\text{,}}z{\text{,}}{{{z}}_{\rm{s}}}{\text{,}}\omega } \right){{\rm{J}}_n}\left( {kr} \right){\rm{d}}k{\text{,}} $$ (1)

    式中,r为震中距,zzs分别为震源和接收点的深度,k为水平波数,ω为角频率,Jnkr)为n阶贝塞尔函数,Frzzsω)是包含衰减因子${{\rm{e}}^{ {\text{-}}c(k{\text{,}} \omega) |z{\text{-}} {z_{\rm{s}}}|}}$的项,其中ckω)为与波数k和频率ω有关的实部为正的复数因子。当震源和接收点的深度相近或相同时,衰减因子作用很小,在介质响应函数和贝塞尔函数的共同作用下,整个被积函数衰减得很慢。为了得到精确的结果而不得不取很大的积分上限,这样就使得数值积分的效率很低。但是,当波数k大于某一阈值kc时,积分值随着积分上限的增大而在某一极限值附近振荡,且其包络线为单调光滑的曲线。将波峰与波谷值抽取出来组成一个慢收敛序列,利用重复平均法便可快速求出其收敛值(Zhang et al,2003)。因此,重要的一步便是求出阈值kc。以往的工作中,阈值kc由经验公式给出(张海明等,2001),

    $$ {{{k}}_{\rm{c}}} {\text{=}} 1.5\frac{{\sqrt {{\omega ^2} {\text{+}} \omega _{\rm{I}}^2} }}{{{v_{{\rm{min}}}}}}{\text{,}} $$ (2)

    式中,vmin为介质模型的最小速度,ωωI分别为角频率和虚频率。然而为了尽量照顾到所有情况,通过经验公式求出的kc可能会取得过大,从而牺牲了计算效率。

    在广义反透射系数法中(Chen,1999),为描述各层之间地震波场振幅的关系,引入了修正反透射系数。由于P-SV波和SH波的修正反透射系数和广义反透射系数具有相同的形式,且我们所研究的瑞雷波为P波和SV波耦合所形成,故下文的讨论只针对P-SV波进行。对于一个N层的分层介质,P-SV波的修正反透射系数满足[Chen (1999)中式(41)]

    $$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{T}}_{\rm{d}}^{{{(j)}}}} & {{\boldsymbol{R}}_{{\rm{ud}}}^{{{(j)}}}}\\ {{\boldsymbol{R}}_{{\rm{du}}}^{{{(j)}}}} & {{\boldsymbol{T}}_{\rm{u}}^{{{(j)}}}} \end{array}} \right]{\rm{ {\text{=}} }}{\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(j {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{12}}}^{{{(j)}}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(j {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{22}}}^{{{(j)}}}} \end{array}} \right]^{{\rm{ {\text{-}} 1}}}}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(j)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{12}}}^{{{(j {\text{+}} 1)}}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(j)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{22}}}^{{{(j {\text{+}} 1)}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_{\rm{d}}^{\left( {{j}} \right)}{\rm{(}}{{{z}}^{{{(j)}}}}{\rm{)}}} & {\rm{0}}\\ {\rm{0}} & {{\boldsymbol{K}}_{\rm{u}}^{\left( {{{j {\text{+}} 1}}} \right)}{\rm{(}}{{{z}}^{{{(j)}}}}{\rm{)}}} \end{array}} \right], $$ (3)

    式中,j表示层数,取值为1,2,···,N−1,R为修正反射系数,T为修正透射系数。当jN时,修正反透射系数满足[Chen (1999)中式(42)]

    $$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{T}}_{\rm{d}}^{{{(N)}}}}\\ {{\boldsymbol{R}}_{{\rm{du}}}^{{{(N)}}}} \end{array}} \right]{\rm{ {\text{=}} }}{\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(N {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{12}}}^{{{(N)}}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(N {\text{+}} 1)}}}} & { {\text{-}} {\boldsymbol{E}}_{{\rm{22}}}^{{{(N)}}}} \end{array}} \right]^{ {\text{-}} {\rm{1}}}}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_{{\rm{11}}}^{{{(N)}}}{\boldsymbol{K}}_{\rm{d}}^{\left( {{N}} \right)}{\rm{(}}{{{z}}^{{{(N)}}}}{\rm{)}}}\\ {{\boldsymbol{E}}_{{\rm{21}}}^{{{(N)}}}{\boldsymbol{K}}_{\rm{d}}^{\left( {{N}} \right)}{\rm{(}}{{{z}}^{{{(N)}}}}{\rm{)}}} \end{array}} \right]{\text{,}} $$ (4)

    进一步,为描述各层介质中地震波场振幅对某一层中地震波场振幅的影响,引入广义反透射系数,其表达式为:

    j=1,2,···,s-1时,满足[Chen (1999)中式(45)]

    $$ \left\{ \begin{aligned} &{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}{\rm{ {\text{=}} }} {\text{-}}{\boldsymbol{E}}{_{{\rm{21}}}^{\left( {\rm{1}} \right)^{\text{-}1}}}{\boldsymbol{E}}_{{\rm{22}}}^{\left( {\rm{1}} \right)}{\boldsymbol{K}}_{\rm{u}}^{\left( {\rm{1}} \right)}{\rm{(}}{{{z}}^{{\rm{(0)}}}}{\rm{)}}{\text{,}}\\ &{\hat{\boldsymbol T}}_{\rm{u}}^{\left( {{j}} \right)}{\rm{ {\text{=}} }}{[{\boldsymbol{I}} {\text{-}} {\boldsymbol{R}}_{{\rm{du}}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{{j}} {\text{-}} 1} \right)}]^{ {\text{-}}1}}{\boldsymbol{T}}_{\rm{u}}^{\left( {{j}} \right)}{\text{,}}\\ &{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{j}} \right)} {\text{=}} {\boldsymbol{R}}_{{\rm{ud}}}^{\left( {{j}} \right)} {\text{+}}{\boldsymbol{T}}_{\rm{d}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{{j}} {\text{-}} 1} \right)}{\hat{\boldsymbol T}}_{\rm{u}}^{\left( {{j}} \right)}{\text{;}} \end{aligned} \right. $$ (5)

    jss+1,···,N时,满足[Chen (1999)中式(46)]

    $$ \left\{ \begin{aligned} &{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{{N}} {\text{+}} 1} \right)}{\rm{ {\text{=}} }}0{\text{,}}\\ &{\hat{\boldsymbol T}}_{\rm{d}}^{\left( {{j}} \right)}{\rm{ {\text{=}} }}{[{\boldsymbol{I}} {\text{-}} {\boldsymbol{R}}_{{\rm{ud}}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{{j}} {\text{+}} 1} \right)}]^{ {\text{-}} 1}}{\boldsymbol{T}}_{\rm{d}}^{\left( {{j}} \right)}{\text{,}}\\ &{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{j}} \right)} {\text{=}} {\boldsymbol{R}}_{{\rm{du}}}^{\left( {{j}} \right)} {\text{+}} {\boldsymbol{T}}_{\rm{u}}^{\left( {{j}} \right)}{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{{j}} {\text{+}} 1} \right)}{\hat{\boldsymbol T}}_{\rm{d}}^{\left( {{j}} \right)}{\text{,}} \end{aligned} \right. $$ (6)

    式中,s为源所在层数,${\hat{\boldsymbol R}}$为广义反射系数,${\hat{\boldsymbol T}}$为广义透射系数,I为单位矩阵。在式(5)中存在矩阵求逆的情况,若相应的求逆矩阵存在行列式等于零的情况,此时该项便存在奇异性。式(5)中,${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {{j}} \right)}$${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {{j}} \right)}$由于存在矩阵求逆,因此,这三项均有可能存在某个波数使得分母为零。下文将对这三项作进一步讨论,以寻找最适于用来求解波数积分限阈值kc的项。

    以一个三层介质模型为例(Chen,1993),如表1所示。震源深度z为29.95 km,接收点深度zs为30.05 km。分别求出不同频率下${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {j} \right)}$${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {j} \right)}$各项中含有负一次幂矩阵的行列式的零点,即分别求出${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$中所含${\boldsymbol{E}}_{21}^{\left( 1 \right)}$${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {j} \right)}$中所含${\boldsymbol{I}} {\text{-}} {\boldsymbol{R}}_{{\rm{du}}}^{\left( {j} \right)}{\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( {{j} {\text{-}} 1} \right)}$${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {j} \right)}$中所含${\boldsymbol{I}}{\text{-}} {\boldsymbol{R}}_{{\rm{ud}}}^{\left( {j} \right)}{\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( {{j} {\text{+}} 1} \right)}$这三项的行列式所对应的零点。在这个三层模型中,由于s=1,N=3,则

    表  1  不含低速层的三层地壳模型参数 (Chen,1993
    Table  1.  Parameters of the three-layer crust model without low-velocity-layers(after Chen,1993
    层序 层底深度/km ρ/(g·cm−3 vS/(km·s−1 vP/(km·s−1
    1 20 2.8 3.50 6.0
    2 30 2.9 3.65 6.3
    3 45 3.1 3.90 6.7
    4 3.3 4.70 8.2
    下载: 导出CSV 
    | 显示表格
    $$ {\hat{\boldsymbol R}}_{{\rm{du}}}^{\left( 4 \right)}{\rm{ {\text{=}} }}0{\text{,}} {\hat{\boldsymbol T}}_{\rm{d}}^{\left( {\rm{3}} \right)} {\text{=}} {\boldsymbol{T}}_{\rm{d}}^{\left( {\rm{3}} \right)}{\text{.}} $$ (7)

    第三层广义透射系数${\hat{\boldsymbol T}}$对应的式子中并无矩阵求逆,故只有第一层和第二层的广义透射系数,即${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {\rm{1}} \right)}$${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {\rm{2}} \right)}$,有贡献。由${\hat{\boldsymbol T}}_{\rm{u}}^{\left( {\rm{1}} \right)}$${\hat{\boldsymbol T}}_{\rm{d}}^{\left( {\rm{2}} \right)}$这两项求得的零点有多个,只取其零点的最大值,分别记作kt1kt2,同时将通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点记为kr,则可以得到在不同频率下,各项零点与积分值的关系,如图1所示。

    图  1  不同频率下反透射系数的零点与位移的关系(位移Ur为归一化后的结果)
    图(a)−(d)分别对应f=0.1,0.4,3,5 Hz的不含低速层的三层地壳模型;图(e)和(f)分别对应f=1,5 Hz的含低速层的五层地壳模型
    Figure  1.  The relationship between zeros of R/T coefficient and the displacement in different frequencies with the normalized displacement Ur
    Figs. (a) to (d) correspond to the three-layer crust model without low-velocity-layers when f=0.1,0.4,3 and 5 Hz;Figs. (e) and (f) correspond to the five-layer crust model with low-velocity-layers when f=1 Hz and 5 Hz

    图1中可以看出,对于不同的频率取值,通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点kr均出现在积分值产生振荡行为的位置附近。再观察由透射系数求得的kt1kt2的情况:在f=0.1 Hz的情况下(图1a),kt1的值与kr相近,而kt2处于积分值产生振荡行为之前,kt2不可作为积分限阈值;在f=0.4 Hz的情况下(图1b),kt2能够满足作为积分限阈值的条件,但是kt2并不适用;在f=3 Hz的情况下(图1c),kt2的值与kr相近,而kt1处于积分值产生振荡行为之前,kt1不可作为积分限阈值;在f=5 Hz的情况下(图1d),kt2适用,kt1不适用。再考虑一个含低速层的五层地壳模型(何耀锋等,2006),如表2所示,其结果如图1ef所示。在f=1 Hz的情况下(图1e),kt2适用,而kt1不适用;在f=5 Hz的情况下(图1f),kt1可作为积分限阈值,而kt2远大于kt1kr,若取作积分限阈值,则计算效率相对较低,因此作为积分限阈值不合适。在柱坐标系下,位移可以分解为三个分量,分别为径向分量Ur,切向分量Uf和垂直分量Uz图1只显示了Ur的实部,对于Ur的虚部和另外两个方向位移分量的实部和虚部均可得到相似的结果。可以看出在所有情况中,通过反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求得的零点kr均满足其作为积分限阈值的要求,而由透射系数求得的零点只在某些情况下满足。因此我们可以认为通过${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$获得的零点kr作为积分限阈值更为合理。由此我们便找到了求解积分限阈值的另一个方法,即求解广义反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$中使逆矩阵的行列式为零的点kr,将kr作为积分限阈值。

    表  2  含低速层的五层地壳模型的参数(何耀锋等,2006
    Table  2.  Parameters of the five-layer crust model with low-velocity-layers(after He et al, 2006
    层序 层底深度/km ρ/(g·cm−3 vS/(km·s−1 vP/(km·s−1
    1 20 2.8 3.50 6.0
    2 30 2.9 3.65 6.3
    3 45 3.5 2.50 4.9
    4 70 4.0 3.80 6.5
    5 90 3.2 4.50 7.6
    6 2.4 5.60 9.3
    下载: 导出CSV 
    | 显示表格

    ${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$进一步进行研究,得到各项的含义为

    $$ {\boldsymbol{E}}_{{\rm{21}}}^{{\rm{(1)}}} {\text{=}} \frac{1}{\omega }\left[ {\begin{array}{*{20}{c}} { {\text{-}} 2{\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\gamma ^{\left( 1 \right)}}} & { {\text{-}} {\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}}\\ { {\text{-}} {\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}} & { {\text{-}} 2{\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\nu ^{\left( 1 \right)}}} \end{array}} \right] {\text{,}} $$ (8)
    $$ {\boldsymbol{E}}_{{\rm{22}}}^{{\rm{(1)}}} {\text{=}} \frac{1}{\omega }\left[ {\begin{array}{*{20}{c}} {2{\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\gamma ^{\left( 1 \right)}}} & {{\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}}\\ { {\text{-}} {\alpha ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}{\chi ^{\left( 1 \right)}}} & {{\text{-}} 2{\beta ^{\left( 1 \right)}}{\mu ^{\left( 1 \right)}}k{\nu ^{\left( 1 \right)}}} \end{array}} \right] {\text{,}} $$ (9)
    $$ {\boldsymbol{K}}_u^{\left( 1 \right)} {\text{=}} \left[ {\begin{array}{*{20}{c}} {{{\rm e}^{ {\text{-}} {\gamma ^{\left( 1 \right)}}\left( {z {\text{-}} {z^{\left( 0 \right)}}} \right)}}} & 0\\ 0 & {{{\rm e}^{{\text{-}} {\nu ^{\left( 1 \right)}}\left( {z {\text{-}} {z^{\left( 0 \right)}}} \right)}}} \end{array}} \right]{\text{,}} $$ (10)
    $$ {\nu ^{\left( 1 \right)}} {\text{=}} \sqrt {{k^2}{\text{-}} \frac{\omega }{{{{\left(\, {{\beta ^{\left( 1 \right)}}} \right)}^2}}}} {\text{,}}{\gamma ^{\left( 1 \right)}} {\text{=}} \sqrt {{k^2} {\text{-}} \frac{\omega }{{{{\left( {{\alpha ^{\left( 1 \right)}}} \right)}^2}}}}{\text{,}}(R( {{\nu ^{\left( 1 \right)}}} ){\text{>}}0 {\text{,}}R( {{\gamma ^{\left( 1 \right)}}}) {\text{>}} 0) {\text{,}} $$ (11)

    带入式(5)中第一个子式进行计算,可得瑞雷函数的表达式为

    $$ \begin{align} R\left( k \right) &{\text{=}} 4{k^2}\sqrt {\displaystyle\frac{{{{({\alpha ^{\left( 1 \right)}})}^2}{k^2} {\text{-}} {\omega ^2}}}{{{{({\alpha ^{\left( 1 \right)}})}^2}}}} {(\,{\beta ^{\left( 1 \right)}})^4}\sqrt {\displaystyle\frac{{{{(\,{\beta ^{\left( 1 \right)}})}^2}{k^2} {\text{-}} {\omega ^2}}}{{{{(\,{\beta ^{\left( 1 \right)}})}^2}}}} {\text{-}} 4{({\beta ^{\left( 1 \right)}})^4}{k^4}{\text{+}}\\ & \ 4{(\,{\beta ^{\left( 1 \right)}})^2}{k^2}{\omega ^2} {\text{-}} {\omega ^4} {\text{=}} {\left( {1 {\text{-}} 2x} \right)^2} {\text{-}} 4x\sqrt {m {\text{-}} x} \sqrt {1 {\text{-}} x}\;\;{\text{,}} \end{align} $$ (12)

    式中,$x {\text{=}} {({\beta ^{\left( 1 \right)}})^2}{\left( {k/\omega } \right)^2}$$m {\text{=}} {( {{\beta ^{(1)}}})^2}\Bigr/{( {{\alpha ^{(1)}}} )^2}$。当m的取值范围为(0,0.5)时,该瑞雷函数有且仅有一个大于1的实数根(Liu et al,2016)。因此广义反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$求逆矩阵的行列式的零点可以通过求解自由表面瑞雷函数的零点来获得。

    仍然采用表1所示的三层地壳模型。将各参数带入式(8),求出瑞雷函数的零点。同时,画出频率域下三个位移分量积分值随积分上限的变化曲线,得到瑞雷函数零点与积分值的关系,如图2所示。可以清楚地看出,积分值的震荡行为产生于瑞雷函数零点附近。因此有理由认为,通过求瑞雷函数零点来求解kc的方式是正确的。至此,本文便可以给出基于理论分析求解kc的方法。在计算含有衰减因子的积分时,若需求得积分限阈值kc,即可通过求取广义反透射系数求逆矩阵行列式的零点的方式求得,并且在实际计算中,为了保证构成慢收敛序列的最初几个波峰和波谷值的稳定,kc的取值将会稍大于零点值。我们将这种通过求瑞雷函数零点来寻求积分限阈值的方法定义为求零点法。

    图  2  f=2 Hz时波速递增的三层地壳模型下位移三分量的实部、虚部图(位移U为归一化后的结果)
    图(a),(b)分别为径向位移Ur的实部和虚部;图(c),(d)分别为切向位移Uf的实部和虚部;图(e),(f)分别为垂向位移Uz的实部和虚部
    Figure  2.  The real part and image part of three components of displacement in the three-layer crust model without low-velocity-layers with the normalized displacement U when f=2 Hz
    Figs. (a) and (b) are the real part and image part of the radial direction displacement (Ur),respectively;Figs. (c) and (d) are the real part and image part of the tangential displacement (Uf),respectively;Figs. (e) and (f) are the real part and image part of the vertical displacement (Uz),respectively

    与之前通过经验公式找积分限阈值kc的方法相比,通过求零点法来寻求kc具有一定的优势,其寻求的kc更为准确,而且针对某些特殊的情形,求零点法比经验公式法更实用,它可以更精确地找到临界点,并在保证计算精度的前提下具有更高的计算效率。

    仍然考虑一个如表1所示的三层地壳模型,震源仍置于地下深度为29.95 km处,接收点仍置于地下深度为30.05 km处,分别通过求零点法和经验公式法求出积分限阈值。为了区别两种方法求出的积分限阈值,我们将通过经验公式求得的积分限阈值记为ke,将通过求零点法求得的积分限阈值记作kc。结果如图3所示,可以清楚地看到:两种方法在求得积分限阈值后,积分值均具有震荡行为;由求零点法得到的积分限阈值kc更接近积分值开始产生震荡的位置,在运算的过程中具有更高的效率;而由经验公式法得到的积分限阈值kekc大,意味着需要更大的计算量。

    图  3  波速递增的三层地壳模型在不同频率下kekc的关系(位移Ur为归一化后的结果)
    Figure  3.  The relationship between kc and ke at different frequencies in the three-layer crust model without low-velocity-layers with the normalized displacement Ur

    若考虑如表2所示存在低速层的情况,得到的结果如图4所示。可以看出,在存在低速层的情况下,由经验公式法得到的积分限阈值ke不得不取得很大,说明在该情况下,经验公式法的适用性较差。同时,也可以看到由求零点法得到的积分限阈值kc仍然出现在产生震荡行为的积分值附近,与本文结果十分符合。

    图  4  f=5 Hz时含低速层的五层地壳模型下kekc的关系(位移Ur为归一化后的结果)
    Figure  4.  The relationship between ke and kc in the five-layer crust model with low-velocity-layers while the displace-ment has been normalized when f=5 Hz

    为了进一步检验结果的正确性,采用求零点法和经验公式法(Zhang et al,2003)分别求得理论地震图,并对结果进行对比。仍然考虑如表2所示的含低速层的5层地壳模型,震源和接收点均位于地表,震中距为200 km,结果如图5所示(图中的位移为归一化后所得)。可以看到,两种方法得到的结果相同。但是利用求零点法计算理论地震图需要12.74 s,所用计算机处理器为双核i3CPU,内存为4G;而用原有的经验公式法计算理论地震图则需要14.693 s,求零点法的加速效率为13.3%。该结果充分证明了求零点法的正确性,同时也说明该方法的确能够提升理论地震图的计算效率。

    图  5  利用求零点法和经验公式法分别获得的理论地震图径向位移Ur (a)、切向位移Uf (b)和垂向位移Uz (c)三分量归一化图
    Figure  5.  Seismogram calculated by for-zero method and empirical-formula method respectively to obtain the radial direction displacement Ur,the tangential displacement Uf,and the vertical displace-ment Uz,which are normalized

    通过比较不同深度的积分图像(图6)可以看出:当震源与接收点的相对深度越小,即震源与接收点垂直方向相对距离越小,积分收敛得越慢;当相对深度越大,即震源与接收点垂向相对距离越大,积分收敛则越快。因此,通过求零点法寻找积分限阈值的方法适用于处理震源与接收点深度相近的情况。

    图  6  波速递增的三层地壳模型下位移在不同震源与接收点相对深度h时的变化图(位移Ur为归一化后的结果)
    Figure  6.  Relative depth between source and receiver in the three-layer crust model without low-velocity-layers while the displacements Ur have been normalized

    在利用峰谷平均法计算合成理论地震图的波数积分时,找到一个合适的kc,是此方法应用的前提。本文从广义反透射系数法出发,发现由广义反射系数${\hat{\boldsymbol R}}_{{\rm{ud}}}^{\left( 0 \right)}$中求逆矩阵行列式的零点作为积分限阈值最合适,并由该行列式得到了分层介质中自由表面的瑞雷函数。同时,通过对瑞雷函数零点的求解,从图像上寻求了kc与瑞雷函数零点的关系。在不同频率、不同层数的情况下,采用本文的方法均能很好地解决求解积分限阈值的问题。同时,只要震源和接收点的相对深度不大,也即垂向相对距离使得积分值的收敛足够慢,则无论震源和接收点是否在同一层,皆不会影响求零点法的应用。与前人所用的由经验公式求解积分限阈值的方法相比,求零点法能够获取更小的积分限阈值,这样使得我们能够以更高的效率计算理论地震图;同时对于一些用经验公式求解不适用的特殊情形,例如存在低速层的情况,使用求零点法得到的阈值计算效率更高。

  • 图  1   IWTH15台站地表(上)和井下(下)的模拟记录

    A点代表非线性发生的阈值,BCD点分别代表等效剪切波速的明显下降点、最小值点和非线性恢复点,图2同此

    Figure  1.   The simulated records of surface (upper) and borehole (lower) at the station IWTH15

    The red box represents one deconvolution time window,and the black box represents the next window. The point A represents the threshold of nonlinearity, the points BC and D represent the location of obvious decrease,minimum value and nonlinearity recovery of equivalent shear wave velocity,respectively,which are the same in Fig. 2

    图  2   不同窗长(4 s,6 s,8 s,10 s和12 s)下基于移动窗解卷积的等效剪切波速与时间窗内PGA的关系

    蓝色圆点表示以该点为中心的时间窗内解卷积的等效剪切波速,黑色圆点为每个时间窗内的PGA

    Figure  2.   The relationship between the equivalent shear wave velocity by moving-window deconvolution method of different window lengths (4 s,6 s,8 s,10 s and 12 s) and PGA in each time window

    The blue dots represent the equivalent shear wave velocity by the deconvolution in the time window centered on this point,and the black dots represent the PGA in each time window

    图  3   (a) FSKH11台站地表和井下东西分量的加速度时程;(b) 基于移动窗解卷积法的等效剪切波速vS随窗内PGA的变化;(c) 基于移动窗谱比法的峰值频率随窗内PGA的变化;(d) 基于移动窗谱比法的地表和井下的谱比随时间的变化

    图(a)和(b)中,A点为非线性发生的阈值,BCD点分别代表等效剪切波速的局部最小值、小幅恢复点和最小值,E点为强震动后等效剪切波速的恢复值;图(b)中红色框为等效剪切波速vS振荡区域,红色圆圈为该区域的vS平均值;图(d)中白色圆点表示峰值频率的变化轨迹

    Figure  3.   (a) The surface and borehole accelerations of EW component at the station FSKH11;(b) Temporal changes of equivalent shear wave velocity with PGA of each time window;(c) Temporal changes of peak frequencies with PGA of each time window;(d) Color-coded surface/borehole spectral ratios plotted against time based on moving window spectral ratios method

    In Figs. (a) and (b),the point A represents the threshold of nonlinearity,the points BC and D represent the local minimum value,small recovery and minimum value of equivalent shear wave velocity vS,the point E represents the recovery vS after the strong motion. In Fig. (b) the red frame is the oscillation area of vS,the red circle is the average vS;in Fig. (d) the white dot represents the trajectory of the peak frequency

    图  4   FSKH11台站在2011年日本东北大地震前3个月内16次弱震记录(a)和震后10天内25次弱震记录(b)的东西分量解卷积波形

    Figure  4.   The deconvolved waveforms of EW component from 16 weak earthquakes recorded within three months before the 2011 Tohoku earthquake (a) and 25 weak earthquakes recorded within ten days after the earthquake (b) at the station FSKH11

    图  5   2011年日本东北大地震中台站IBRH12 (a),IWTH21 (b),IWTH27 (c)和FKSH19 (d)的场地非线性时变过程识别结果(各子图意思同图3,AF点分别表示两种方法识别到的非线性阈值的位置)

    Figure  5.   The identification results of the temporal changes in site nonlinearity at the stations IBRH12 (a),IWTH21 (b),IWTH27 (c) and FKSH19 (d) during 2011 Tohoku earthquake (The meanings of subfigures are the same as Fig. 3,and the points A and F denote the position of the nonlinear threshold by two methods)

    图  5   2011年日本东北大地震中台站IBRH20 (e),MYGH04 (f)和MYGH10 (g)的场地非线性时变过程识别结果(各子图意思同图3,AF点分别表示两种方法识别到的非线性阈值的位置)

    Figure  5.   The identification results of the temporal changes in site nonlinearity at the stations IBRH20 (e),MYGH04 (f) and MYGH10 (g) during 2011 Tohoku earthquake (The meanings of subfigures are the same as Fig. 3,and the points A and F denote the position of the nonlinear threshold by two methods)

    图  6   (a) 非线性阈值与vS30的关系;(b) 等效剪切波速vS下降比与PGA的关系;(c) vS恢复比与PGA的关系

    Figure  6.   (a) Relationship between the nonlinear threshold and vS30;(b) Relationship between the decrease ratio of equivalent shearing velocity vS and PGA;(c) Relationship between the vS recovery ratio and PGA

    表  1   IWTH15台站的场地参数

    Table  1   Site parameters of the station IWTH15

    序号层厚/m深度/m土性土类vP/(m·s−1vS/(m·s−1ρ/(g·cm−3
    144凝灰角砾岩砂土4801501.451
    2812凝灰质砂岩砂土17803602.014
    32638凝灰岩砂土17804502.014
    472110凝灰角砾岩基岩18705402.039
    512122凝灰角砾岩基岩21606802.113
    下载: 导出CSV

    表  2   土的动剪切模量比G/Gmax和阻尼比λ与剪应变γ的关系

    Table  2   Relationship between dynamic shear modulus ratio G/Gmax and damping ratio λ and shear strain γ of soil

    剪应变γ砂土基岩
    G/GmaxλG/Gmaxλ
    0.000 5%0.976 0%1.35%1.0%0.80%
    0.001 0%0.954 4%1.76%1.0%1.00%
    0.005 0%0.838 9%3.46%1.0%1.50%
    0.010 0%0.800 8%4.67%1.0%2.10%
    0.050 0%0.398 5%8.44%1.0%3.00%
    0.100 0%0.276 3%9.83%1.0%3.60%
    0.500 0%0.086 8%11.75%1.0%4.60%
    1.000 0%0.043 1%12.12%1.0%5.40%
    下载: 导出CSV

    表  3   利用DEEPSOIL进行时域非线性分析时所需拟合参数

    Table  3   Nonlinear fitting parameters in time-domain nonlinear analysis by DEEPSOIL

    土类小应变阻尼比参考压应力/MPa参考有效应变拟合参数
    βsbd
    砂土0.269 4%0.180.071.470.7200
    基岩0.778 4%0.180.350.150.7200
    注:βs为土的双曲线本构模型参数,b为参考剪应变拟合参数,d为小应变阻尼比拟合参数。
    下载: 导出CSV

    表  4   2011年日本东北MW9.0大地震中8个台站的信息及非线性识别结果

    Table  4   Information and nonlinear identification results of eight stations during the 2011 Tohoku MW9.0 earthquake

    序号台站
    编号
    北纬
    东经
    测井
    深度
    /m
    vS30
    /(m·s−1
    场地
    类别
    台站
    土层vS
    /(m·s−1
    地表
    PGA
    /(cm·s−2
    井下
    PGA
    /(cm·s−2
    非线性阈值/(cm·s−2vS
    降比
    峰值
    频率
    下降比
    vS
    复比
    移动窗
    谱比
    移动窗
    解卷积
    1 IBRH12 36.8 140.3 200 486 C 968 519.2 111.8 96 97 18% 38% 92%
    2 IWTH21 39.5 141.9 100 498 C 1085 392.4 54.3 42 40 23% 19% 94%
    3 IWTH27 39.0 141.5 100 670 C 1366 576.6 94.2 76 75 17% 23% 90%
    4 FKSH11 37.2 140.3 115 240 D 439 386.3 152.8 60 55 20% 33% 89%
    5 FKSH19 37.5 140.7 100 338 D 843 822.4 277.3 70 64 24% 39% 90%
    6 IBRH20 35.8 140.7 923 244 D 677 172.6 59.2 - 56 3% - 99%
    7 MYGH04 38.8 141.3 100 850 B 1635 393.6 92.1 58 59 37% 48% 76%
    8 MYGH10 37.9 140.9 205 348 D 585 786.9 142.9 - 97 13% - 96%
    * 场地类别引自National Earthquake Hazards Reduction Program (2015)vS代表等效剪切波速。
    下载: 导出CSV
  • 陈学良. 2006. 土体动力特性、复杂场地非线性地震反应及其方法研究[D]. 哈尔滨: 中国地震局工程力学研究所: 3−6.

    Chen X L. 2006. Study on Soil Dynamic Characteristics, Nonlinear Seismic Response of Complex Site and Its Methods[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration: 3−6 (in Chinese).

    李晓飞,孙锐,袁晓铭,李波. 2015. 现有等效线性化分析程序在实际软场地计算结果方面的比较[J]. 自然灾害学报,24(4):56–62.

    Li X F,Sun R,Yuan X M,Li B. 2015. Comparison of existing equivalent linear analysis program in calculation results of actual soft site[J]. Journal of Natural Disasters,24(4):56–62 (in Chinese).

    李小军. 1993. 非线性场地地震反应分析方法的研究[D]. 哈尔滨: 中国地震局工程力学研究所: 2−10.

    Li X J. 1993. Study on the Method for Analysing the Earthquake Response of Nonlinear Site[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration: 2−10 (in Chinese).

    苗雨,施洋,王苏阳,王海云. 2018. 基于竖向台阵地震记录的非线性场地反应研究[J]. 自然灾害学报,27(6):51–58.

    Miao Y,Shi Y,Wang S Y,Wang H Y. 2018. Assessing nonlinear soil behavior using vertical array data:A case at TCGH16 station from KiK-Net in Japan[J]. Journal of Natural Disasters,27(6):51–58 (in Chinese).

    任叶飞. 2015. 基于强震动记录的汶川地震场地效应研究[J]. 国际地震动态,(4):37–38. doi: 10.3969/j.issn.0235-4975.2015.04.008

    Ren Y F. 2015. Study on site effect in the Wenchuan earthquake using strong-motion recordings[J]. Recent Developments in World Seismology,(4):37–38 (in Chinese).

    王海云. 2014. 基于强震观测数据的土层场地反应的研究现状[J]. 地震工程与工程震动,34(4):42–47.

    Wang H Y. 2014. A review of study on soil site response estimating from strong motion data[J]. Earthquake Engineering and Engineering Dynamics,34(4):42–47 (in Chinese).

    王苏阳. 2017. 基于日本KiK-net地震动数据的场地反应研究[D]. 哈尔滨: 中国地震局工程力学研究所: 75−107.

    Wang S Y. 2017. Study on Site Effect Using Ground Motion Data From KiK-Net in Japan[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration: 75−107 (in Chinese).

    Beresnev I A,Wen K L. 1996. Nonlinear soil response: A reality?[J]. Bull Seismol Soc Am,86(6):1964–1978.

    Bonilla L F,Tsuda K,Pulido N. 2011. Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku earthquake[J]. Earth Planets Space,63:50.

    Bonilla L F,Guéguen P,Ben-Zion Y. 2019. Monitoring coseismic temporal changes of shallow material during strong ground motion with interferometry and autocorrelation[J]. Bull Seismol Soc Am,109(1):187–198. doi: 10.1785/0120180092

    Federico D, Bonilla L F, Foti S. 2016. In situ shear modulus reduction computation using seismic interferometry by deconvolution from borehole and surface data: Theory and examples[C]//Proc. 5th IASPEI/IAEE International Symposium: Effects of Surface Geology in Seismic Motion. Taiwan, China: NARLabs: 1−14.

    Hartzell S. 1998. Variability in nonlinear sediment response during the 1994 Northridge,California,earthquake[J]. Bull Seismol Soc Am,88(6):1426–1437.

    Hashash Y M A,Park D. 2001. Non-linear one-dimensional seismic ground motion propagation in the Mississippi embayment[J]. Engineering Geology,62(1/2/3):185–206.

    National Earthquake Hazards Reduction Program. 2015. Recommended Provisions for Seismic Regulations for New Buildings and Other StructuresFEMA P-1050)[S]. Washington D.C.: Building Seismic Safety Council: 14−18.

    National Research Institute for Earth Science and Disaster Resilience. 2020. Strong-motion seismograph networks (K-NET, KiK-net)[EB/OL]. [2020-10-21]. https://www.kyoshin.bosai.go.jp/kyoshin/data/index_en.html.

    Noguchi S, Sasatani T. 2008. Quantification of degree of nonlinear site response[C]//Proceedings of the 14th World Conference on Earthquake Engineering. Beijing: Seismological Press: 1−8.

    Pavlenko O,Irikura K. 2002. Nonlinearity in the response of soils in the 1995 Kobe earthquake in vertical components of records[J]. Soil Dynam Earthq Eng,22(9/12):967–975.

    Ren Y F,Wen R Z,Yamanaka H,Kashima T. 2013. Site effects by generalized inversion technique using strong motion recordings of the 2008 Wenchuan earthquake[J]. Earthq Eng Eng Vibrat,12(2):165–184. doi: 10.1007/s11803-013-0160-6

    Rubinstein J L,Beroza G C. 2004a. Evidence for widespread nonlinear strong ground motion in the MW6.9 Loma Prieta earthquake[J]. Bull Seismol Soc Am,94(5):1595–1608. doi: 10.1785/012004009

    Rubinstein J L,Beroza G C. 2004b. Nonlinear strong ground motion in the ML5.4 Chittenden earthquake:Evidence that preexisting damage increases susceptibility to further damage[J]. Geophys Res Lett,31(23):L23614.

    Sawazaki K,Sato H,Nakahara H,Nishimura T. 2006. Temporal change in site response caused by earthquake strong motion as revealed from coda spectral ratio measurement[J]. Geophys Res Lett,33(21):L21303. doi: 10.1029/2006GL027938

    Sawazaki K,Sato H,Nakahara H,Nishimura T. 2009. Time-lapse changes of seismic velocity in the shallow ground caused by strong ground motion shock of the 2000 western-Tottori earthquake,Japan,as revealed from coda deconvolution analysis[J]. Bull Seismol Soc Am,99(1):352–366. doi: 10.1785/0120080058

    Wen K L,Beresnev I A,Yeh Y T. 1994. Nonlinear soil amplification inferred from downhole strong seismic motion data[J]. Geophys Res Lett,21(24):2625–2628. doi: 10.1029/94GL02407

    Wu C Q,Peng Z,Assimaki D. 2009. Temporal changes in site response associated with the strong ground motion of the 2004 MW6.6 mid-Niigata earthquake sequences in Japan[J]. Bull Seismol Soc Am,99(6):3487–3495. doi: 10.1785/0120090108

    Wu C Q,Peng Z G,Ben-Zion Y. 2010. Refined thresholds for non-linear ground motion and temporal changes of site response associated with medium-size earthquakes[J]. Geophys J Int,182(3):1567–1576. doi: 10.1111/j.1365-246X.2010.04704.x

    Wu C Q,Peng Z G. 2011. Temporal changes of site response during the 2011 MW9.0 off the Pacific coast of Tohoku earthquake[J]. Earth Planets Space,63(7):51.

    Yamada M,Mori J,Ohmi S. 2010. Temporal changes of subsurface velocities during strong shaking as seen from seismic interferometry[J]. J Geophys Res:Solid Earth,115(B3):B03302.

图(7)  /  表(4)
计量
  • 文章访问数:  296
  • HTML全文浏览量:  140
  • PDF下载量:  46
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-05-23
  • 修回日期:  2021-10-29
  • 网络出版日期:  2022-03-10
  • 发布日期:  2022-03-17

目录

/

返回文章
返回