不同地区人工爆炸与天然地震记录特征及识别研究

王婷婷, 边银菊, 杨千里, 任梦依

王婷婷,边银菊,杨千里,任梦依. 2021. 不同地区人工爆炸与天然地震记录特征及识别研究. 地震学报,43(4):427−440. DOI: 10.11939/jass.20210169
引用本文: 王婷婷,边银菊,杨千里,任梦依. 2021. 不同地区人工爆炸与天然地震记录特征及识别研究. 地震学报,43(4):427−440. DOI: 10.11939/jass.20210169
Wang T T,Bian Y J,Yang Q L,Ren M Y. 2021. Research on seismic characteristics and identification of artificial explosion in different areas and natural earthquake. Acta Seismologica Sinica43(4):427−440. DOI: 10.11939/jass.20210169
Citation: Wang T T,Bian Y J,Yang Q L,Ren M Y. 2021. Research on seismic characteristics and identification of artificial explosion in different areas and natural earthquake. Acta Seismologica Sinica43(4):427−440. DOI: 10.11939/jass.20210169

不同地区人工爆炸与天然地震记录特征及识别研究

基金项目: 中央级公益性科研院所基本科研专项(DQJB19B17)和核查项目共同资助
详细信息
    通讯作者:

    边银菊: e-mail:bianyinju@cea-igp.ac.cn

  • 中图分类号: P315.3+1

Research on seismic characteristics and identification of artificial explosion in different areas and natural earthquake

  • 摘要: 本文分析了河北怀来多次爆炸、河北三河采石场多次爆炸和低震级天然地震事件的记录特征和时频差异。结果显示:河北怀来爆炸的P波能量强、衰减快、S波发育弱;河北三河采石场爆炸的P波、S波主频均低于怀来爆炸,S波与面波混淆,不同震中距的台站记录低频发育明显;而天然地震的有效频带更宽,频率成分更为复杂。将Pg/Sg谱比判据应用于小震级地震与爆炸的识别中,探索交叉频带谱比对不同地区爆炸的识别。结果表明:高频(>5 Hz)Pg/Sg谱比判据可将研究数据中的爆炸与小震级地震完全区分;与Sg低频(0—2 Hz)有关的交叉频带谱比可对两个不同地区的爆炸进行识别,交叉频带的谱比判据较传统的单一频带谱比判据能够更好地反映出不同类型事件的特征差异。
    Abstract: The differences of the seismic characteristics and frequency of the Huailai explosions, the Sanhe quarry explosions and the natural earthquakes with low magnitude are discussed. The results show that the two different area explosions have obviously different seismic characteristics and frequency distribution, Huailai explosion has stronger P wave energy than S wave and fast attenuation; The main frequencies of P wave and S wave in Sanhe quarry explosion are lower than that in Huailai explosion, S wave and surface wave are confused, and the low frequency developed obviously at different distances; While for natural earthquakes, the effective frequency band is wider and the frequency components are more complex than explosions. Pg/Sg spectral ratios in small-magnitude earthquakes and explosions were studied and cross-band spectral ratios were explored. Results obtained show that the high frequency (>5 Hz)Pg/Sg spectral ratio discriminants can completely distinguish explosions from low magnitude earthquakes; The spectral ratios of the cross-band related to low frequency (0−2 Hz) of Sg can effectively identify explosions in these two areas, Pg/Sg discriminants of the crossed frequency band can better reflect the difference characteristics of different types of events than that of the traditional single frequency band.
  • 爆炸一般涉及化学爆炸、物理爆炸和核爆炸,相较于天然地震,爆炸的震源更为简单。然而,爆炸产生地震效应的能量比例和释放机制仍是一个复杂问题。在爆炸近区,介质变形和爆炸作用之间的物理过程很难从理论上给予描述。当地质条件一定时,炸药的性质、爆炸位置、装药结构、爆炸形式(集中和分散爆炸)以及起爆方式等条件的改变都会引起爆炸所产生的地震波特征的变化(林大超,白春华,2007)。当起爆方式不同时,单点爆炸(地下核爆炸)较多孔延迟爆炸(工业中广泛采用)的记录具有较强的高频特征;除起爆方式外,路径对地震记录影响也很大,天然地震具有一定的埋深,其高频能量衰减速度慢,而工业爆炸发生在地表,其高频能量随距离衰减较快。不同的爆炸条件,不同的传播路径都会引发爆炸记录特征的不同,在地震学中,主要通过地震波记录的时频特征来对不同类型震动事件的性质进行识别(Pomeroy et al,1982Taylor et al,1989Douglas,2013Koper et al,2016)。

    在地震和地下核爆炸的众多识别判据中,对P/S震相幅值比的研究最为深入、应用最为广泛。爆炸理论上属于膨胀源,纵波发育;天然地震属于剪切破裂源,横波发育。针对地下核爆炸和天然地震理论上的纵横波发育差异,Willis等(1963)率先提出用振幅比判据识别地下核爆炸与天然地震,分析了70次内华达核爆炸和附近200多次天然地震,得出78%的天然地震Lgmax/Pmax振幅比大于核爆炸。Murphy和Bennett (1982)指出P波振幅相当的情况下天然地震波其Lg高频含量较爆炸高,高频振幅比的识别效果更好。Fisk等(1996)将0.5—10 Hz分成7个频带,对比了不同核试验场区核爆炸和天然地震的Pn/Lg,Pg/Lg,表明核爆炸与天然地震事件的分类能力随着振幅比频率的增大而增高,当频率大于3 Hz时,可以很好地进行区分。研究显示P/S谱比在高频处可以有效地区分地震与爆炸事件,但对于它的频率依赖性缺乏相应的理论依据。后续不断的研究与探索表明地下核爆炸与地震P/S谱比的频率依赖性主要与P、S波的拐角频率差异和爆炸的波谱超调有关(Xie,Patton,1999Fisk,20062007)。2006年以来朝鲜共进行了六次地下核试验,大部分研究均使用P/S谱比判据对这几次试验的性质进行判定(Kim,Richards,2007Chun et al,2011Zhao et al,2017Kim et al,2018)。

    长期以来,国内外研究人员对爆破地震波的传播理论、爆炸机制、识别依据等方面进行了大量的、富有成效的研究(Fisk et al,2008Bonner et al,2013Pasyanos,2018),但大多数集中在大当量核爆和小当量地震波近场勘探试验中,相关成果也较难推广应用到与采矿炸山等息息相关的工业爆炸中。P/S谱比可以有效地识别天然地震与地下核爆炸,但震源类型识别所面临的挑战不仅仅是识别天然地震和地下核爆炸,还包括对采矿爆炸、突发爆炸灾害、矿震塌陷、滑坡等事件以及其它类型震动源的识别。由于地壳浅层介质的复杂性及非天然震动类型的多样性,对低震级天然地震和工业爆炸的识别更为困难,研究不同工业爆炸记录特征、准确识别天然地震和工业爆炸事件,对研究活断层划分、爆炸安全评估、矿震灾害评估、爆炸震源机制等有着重要的意义。

    北京遥测地震台网初建于1966年邢台地震后,1998年开始对其进行数字化改建,2008年完成了“十五”台网扩建工程,包括了北京、天津和河北等地的107个数字地震台站,可以更好地监测首都圈周边的地震活动。本文从北京数字遥测台网获得的天然地震和非天然地震事件目录显示,首都圈地区每月发生的天然地震事件约数百次,其中1.5≤ML≤4.0的天然地震事件和非天然地震事件(包含爆炸和矿塌)的发震比例约为3∶1 (郑秀芬等,2006)。河北怀来县和三河市是距离北京市中心两处最近的爆炸区,其中怀来爆炸区距北京市中心约80 km,发生过大当量人工爆炸“明灯一号”事件;三河采石场爆炸区距北京市中心约70 km,长期进行采石活动。本文拟分析河北怀来、河北三河不同地区的爆炸与天然地震的地震波记录特征、频率特征的差异,并应用P/S谱比判据对不同事件进行性质识别。

    从北京遥测数字地震台网、国家数字测震台网数据备份中心(Zheng et al,2010)获取了首都圈地震和爆炸事件目录及波形数据,本文选取的怀来爆炸、三河采石场爆炸、周边的天然地震事件及观测台站的位置分布如图1所示,可以看出:人工爆炸成团分布较为集中,其中怀来爆炸位于北京西部山区,三河采石场爆炸位于华北裂陷盆地(徐锡伟等,2002孙海霞等,2018);天然地震则分布较广;研究区内观测台站分布较为均匀,因此本文的研究结果有一定的适用性。具体数据筛选如下:

    图  1  研究所用爆炸、地震及台站位置分布
    Figure  1.  Distribution of  explosions ,seismic events and stations in this article

    1) 中国地震局地球物理研究所于2007年12月12日凌晨3时在河北怀来实施了一次名为“明灯一号”的人工爆破。此次爆破的单次爆炸药量达到50 t,远大于一般工程爆破,北京地震台网测定其近震震级为ML2.8。本文选取此次事件周边地区2002—2015年发生的爆炸事件共47次,震级范围为ML1.4—2.5。

    2) 河北三河地区富产冶金用白云岩矿(赵爱华等,2012),有些矿床达到了大型甚至超大型规模,此地区长期有采石活动,爆炸属性为采石场爆炸。由于城市的快速扩建,2008年以前河北三河采石场每年发生采石场爆炸次数约500次,2008年以后每年发生约50次。本文选取2010年的53次爆炸事件,震级范围为ML1.7—2.5。

    3) 北京及周边天然地震次数较多,本文选取了2010年的54次天然地震,震级ML1.5—3.0,震源深度为4—10 km。通过北京遥测数字地震台网及国家测震台网数据备份中心(Zheng et al,2010)获得相应的北京及周边地区区域台站波形记录,最终整理出有效垂向记录为:怀来爆炸681条,震中距为10—300 km,平均震中距为98 km;三河采石场爆炸1 213条,震中距为11—260 km,平均震中距为112 km;天然地震1 242条,震中距为11—260 km,平均震中距为114 km。

    理论上天然地震是剪切源,有一定的埋深;爆炸是膨胀源,发生在近地表。传播距离相同时,爆炸的高频损失更为严重,面波也较为发育。目前工业爆破经常采用延时微差爆破,虽不能严格将其认为点源,但并不会像天然地震一样具有近乎完全的双力偶震源机制和四象限振幅分布特征。工业活动的延时爆炸和非延时爆炸发生在近地表,埋深浅,故其记录的高频成分随传播路径衰减快,除近场外,爆炸频谱相较于天然地震更缺乏高频能量(Smith,1993唐兰兰,2018)。

    对于怀来爆炸、三河采石场爆炸和天然地震事件,北京、天津、河北、山西等部分区域台站均有信噪比高的宽频带记录。本文选取三次不同事件:怀来爆炸“明灯一号”事件(2007-12-11,ML2.8)、三河采石场爆炸(2010-05-12,ML2.2)、天然地震(2010-05-13,ML2.1)的地震记录(图2),对不同事件尽量选取震中距相近的台站记录进行对比,最小震中距为17.79 km,最大为214 km,波形滤波段为0.5—25 Hz。由于近震中距P,S波干扰严重,本文对震中距10 km以内的数据进行了剔除,地震与爆炸记录震中距均分布在10—300 km以内,其中P波主要震相包含Pg,Pn,Pb和PmP等。全文的波形特征和谱比研究对不同震中距的地震与爆炸均基于理论速度窗口提取P、S波段,未详细区分震相,另外近震震相中以直达波震相Pg和Sg为主,因此全文中以Pg和Sg震相指代P,S波段(图2)。

    图  2  不同类型事件的地震记录波形
    (a)“明灯一号”大当量爆炸;(b)三河采石场爆炸;(c)天然地震。波形上方为事件发生时间及震级,波形右上角为记录台网、台站和震中距,下同
    Figure  2.  Seismic waveforms for different types of events
    (a) “Bright Lamp No.1” large yield explosion;(b) Quarry blast;(c) Earthquake. The time and the magnitude of the events are net on the top of traces,and the station and epicentral distance are on the top right corner of the waveforms,the same below

    分析三次事件按震中距排序的波形,可以看出:① 天然地震事件整段波形频率较高,S波发育,S波振幅大于P波,大约2倍以上;② “明灯一号”爆炸事件P波能量强,衰减快,震中距较近的台站有明显的瑞雷面波,整体S波成分较弱;③ 三河采石场爆炸事件记录不论P波还是S波,与前两类事件相比低频较为发育,S波与面波混淆,不同震中距均低频发育;④ 在不同的震中距范围内,同类事件都表现出类似的特征,如天然地震P波较弱,频率高,两类爆炸P波发育且衰减快等。

    河北三河2002—2016年发生过约3 000次爆炸,其中2008年之前每年约发生500次,2008年之后每年约50次。作为长期的采石矿场,场地介质被大量爆炸事件充分松动,介质损伤严重可能是引起三河爆炸低频发育的主要因素。本文分析了2003年(2003-09-02,ML2.5)、2004年(2004-09-18,ML2.4)、2005年(2005-06-05,ML2.8)三次典型爆炸的地震记录,如图3所示,与2010年爆炸(图2b)对比可知,不同年份不同震中距的三河采石场爆炸记录,都显示出低频发育的特征。大量爆炸事件造成场地介质疏松以致地震记录低频发育,但不是造成三河采石场爆炸记录低频产生的主要因素。

    图  3  2003 (a),2004 (b)和2005 (b)年3次爆炸事件的波形对比
    Figure  3.  Waveform comparison of three explosion events occured in 2003 (a),2004 (b) and 2005 (c)

    频率特征提取是地震数据分析处理中的重要工作之一,为了增强谱估计的稳定性,在经典的谱估计基础上Thomson (1982)提出了多窗口功率谱估计方法,使用的窗函数是一组相互正交的离散椭球序列,又称Slepian窗。利用多个正交窗口获得各自独立的近似功率谱,综合多个窗口结果得到最终功率谱估计,这种功率谱估计具有更大的自由度,并且估计精度和估计波动稳定性方面的效果较好。普通的功率谱估计利用单一窗口,在序列始端和末端均会丢失相关信息,多窗口谱估计法通过增加窗口尽量减少丢失的信息,更好地解决频率泄露等问题,可以获得更为平滑稳定的谱信息(McCoy et al,1998Zanandrea et al,2004王云专等,2009)。

    本文采用Thomson (1982)提出的多窗口谱分析方法计算不同事件的振幅谱信息。首先,对所有波形去均值、去趋势和尖灭等初始化波形处理;其次,分别计算每条记录的P波、S波多窗口谱信息,并进行归一化处理;最后,将同一类的所有记录的归一化频谱叠加在一起,得到不同类型事件具有代表性的P波、S波频率分布图。频谱图不仅能反映出不同类型事件的频率分布差异,而且随着更多样本的加入,还能对地域性的谱比判据进行动态调整。

    不同类型事件的P波、S波归一化振幅谱随频率的分布如图4所示,以最高幅值的0.707倍作为主频率频带分布区间,对不同事件的主频成分进行统计。天然地震P波能量主要集中在2.5—13.2 Hz,衰减较为缓慢;怀来爆炸的P波能量主要分布在2.5—7.2 Hz;三河采石场爆炸的P波主频较低,集中在1.5—5.8 Hz,在6—9 Hz的频带内迅速衰减(图4a)。天然地震S波能量主要频率集中在2.5—9 Hz,衰减缓慢;怀来爆炸的S波能量主要分布在2—3.8 Hz,在4—6 Hz的频带内迅速衰减;三河采石场爆炸的S波主频较低,集中于1—2 Hz,在2—3 Hz的频带内迅速衰减(图4b)。综上可得:① 天然地震频率成分复杂,P、S波主要能量频带均较宽,相较于人工爆炸记录存在高频成分;② 爆炸整体频带较天然地震窄,其中怀来爆炸P波主频在2.5—7.2 Hz,S波主频在2—3.8 Hz;三河采石场爆炸P波主频在1.5—5.8 Hz、S波主频在1—2 Hz;怀来爆炸的P、S波频率均高于三河采石场爆炸(表1)。

    图  4  不同类型事件P波(a)和S波(b)主频分布
    Figure  4.  Main frequency distribution of P (a) and S (b) waves for different types of events
    表  1  天然地震、怀来爆炸、三河爆炸主频分布
    Table  1.  The main frequency distribution of earthquake,Huailai explosion and Sanhe explosion
    事件P波主频分布/HzS波主频分布/Hz
    天然地震 2.5—13.22.5—9.0
    怀来爆炸2.5—7.22.0—3.8
    三河采石场爆炸1.5—5.81.0—2.0
    下载: 导出CSV 
    | 显示表格

    相同条件下,集中瞬时爆炸相较于延时微差爆炸,地震波记录有更多的高频成分(Kim et al,1994)。“明灯一号”大当量科学探测爆炸试验为集中瞬时爆炸,其它怀来爆炸与“明灯一号”地震记录波形及振幅谱特征相似,高频P波能量大于S波,因此本文推测怀来地区爆炸事件大部分为单次爆炸或者多孔瞬时爆炸。延时微差爆炸适用于开挖岩石地基等爆破面积较大的开采工程,是工业活动中的常用方式,三河采石场可能采用该方式进行爆炸作业,故其地震记录与主频分布均含有更多的低频成分,但是不同的场地介质、传播路径也会造成不同地区爆炸地震记录的特征差异。

    理论上爆炸是膨胀源,产生强烈P波。构造地震会产生巨大的剪切能量,即辐射出强烈的S波,基于该思想,研究人员较早应用P/S谱比对天然地震与地下核爆炸进行识别研究,并在不同核试验场取得了有效的应用。工业爆炸与地下核爆炸在爆炸方式、爆炸环境、炸药类型、目的等方面有很大的不同,在P/S谱比识别应用中虽然可以参考核爆的识别方式,但是震源模型的建模上却有很大不同。

    在小震级事件与工业爆炸的识别中,Kim等(19931997)对美国东部和俄罗斯南部的地震和矿震的研究表明,在5—25 Hz频率范围内,高频P/Lg谱值比是一种可靠的分类方法。潘常周等(2007a)检验了P/S谱比判据对新疆地区低震级地震事件的适用性,得到低震级工业爆炸也具有与核爆大致相当的识别效果。在P/S谱比分析中,Pn/Lg分类效果要优于Pg/Lg (Bennett,Murphy,1986)。但低震级事件由于其震级小,有效地震记录的震中距较近,无法产生信噪比高的Pn震相。本研究中的小震级事件,纵波主要震相为Pg,横波主要震相为Sg,在P/S谱比研究中只考虑了Pg与Sg震相,具体计算步骤如下:

    利用震相速度窗口对地震记录提取Pg,Sg有效波段,进行傅里叶变换后,以1 Hz为间隔计算1—20 Hz频带内对应的Pg/Sg幅值比。由于P波S波随震中距的衰减速率不同,需要对P/S谱比进行震中距校正。本文参考Fisk等(2000)对P/S幅值比随距离变化关系进行建模,即

    $$ \mathrm{lg}A{\text{=}}a{\text{+}}b\mathrm{l}\mathrm{g}{\varDelta }{\text{+}}c{\varDelta }{\text{,}} $$ (1)

    式中,$ a $$ b $$ \mathrm{c} $为拟合参数,由实际数据特定频带的震相幅值比A确定,如Pg/Sg (5—8 Hz),$ \varDelta $是震中距。利用天然地震的Pg/Sg谱比随震中距的关系拟合得到式中系数值。取参考距离R0为100 km,将不同距离R上的幅值比均校正到距离R0处(Hartse et al,1997王婷婷,边银菊,2015),即:$A_{{{R}}}^{\rm{rev}} $AR$A_{100}^{\rm{o}} $$A_{R}^{\rm{o}} $(其中,AR$A_{{{R}}}^{\rm{rev}} $分别为校正前、后幅值,$A_{100}^{\rm{o}} $$A_{{R}}^{\rm{o}} $分别为100米和当前距离的理论幅值)。天然地震的P、S振幅随着方位的不同存在一定的差异,为了减少方位的影响,本文对每次事件计算了多台平均的幅值比结果。衰减校正后,1—20 Hz多台平均的Pg/Sg谱比如图5所示。

    图  5  天然地震与爆炸Pg/Sg谱比值分布
    Figure  5.  Pg/Sg spectral ratio value of earthquakes and explosions

    分析天然地震与人工爆炸Pg/Sg谱比值随频率分布可知:频率大于5 Hz时Pg/Sg谱比可将此批天然地震与人工爆炸完全区分,这主要因为爆炸P、S波在5 Hz以上能量差异较大(图4),P波能量约5 Hz时达到最大并开始衰减,而S波能量在大于5 Hz时能量最弱,天然地震P、S波能量随频带的分布差异较小,因此当人工爆炸记录频率大于5 Hz时Pg/Sg谱比值大于天然地震。其中怀来爆炸与三河采石场爆炸表现出相近的谱比分布特征,无法通过同频带谱比对这两个区域爆炸进行识别。

    图5中2005年3月12日发生在怀来的一次爆炸记录其各频带谱比值均落在地震群里,地震目录不排除有误识的情况,因此对此次事件性质进行进一步核实。此事件震中附近既发生过爆炸事件也发生过天然地震,选取两类不同事件分别与2005-03-12事件进行波形对比,三次事件参数列于表2。北京西拨子台(XBZ)记录的三次事件波形如图6所示。

    表  2  三次对比事件参数
    Table  2.  Parameters of three comparison events
    事件性质发震时间 (UTC)震中位置ML深度/km
    年-月-日 时:分:秒北纬/°  东经/°
    怀来爆炸2004-03-02 08:03:35.640.50  115.471.8--
    可疑事件2005-03-12 13:46:27.440.51  115.461.6--
    天然地震2010-07-10 09:53:43.040.47  115.551.87
    下载: 导出CSV 
    | 显示表格
    图  6  XBZ台站的不同事件地震波形记录(滤波频带为:0.5—25 Hz)
    (a) 2004-03-02怀来爆炸;(b) 2005-03-12可疑事件;(c) 2010-07-10天然地震事件。左上角为事件震级和震中距
    Figure  6.  Seismic waveforms of different events recorded by XBZ (filtered:0.5—25 Hz)
    (a) 2004-03-02 Huailai explosion;(b) 2005-03-12 Suspicious event ; (c) 2010-07-10 natural earthquakeThe magnitude and epicentral distance are on the top left corner of the wavforms

    图6可以看出,2004-03-02怀来爆炸三分向P波均发育,垂直向S波能量弱,整个波段频率成分简单;2010-07-10天然地震三分向S波能量远大于P波,频率成分复杂;而2005-03-12可疑事件整体P、S波包络形态、震相细节均与2010-07-10天然地震极为相似。

    可疑事件震中位置既不是天然地震频发区,也不是人工爆炸高发区,无法通过震中位置来判定其属性。通过分析对比震中附近天然地震、人工爆炸多台记录波形,发现可疑事件初动清晰的几个台站初动均向下,不同震中距地震记录P波发育较S波弱,具备天然地震记录特征。因此通过XBZ台站记录的震相细节波形对比、初动方向、多台地震记录分析,以及P/S谱比值分布,本文推测2005-03-12可疑事件属性应为天然地震。

    在同一频带内计算P/S谱比,可以先验地消除许多潜在的困难,包括P波与S波之间与频率相关的影响(如地震矩)等,但将其限制在一个频带内,同时也丢失了很多信息。Hartse等(1997)注意到高频P与低频S谱比也能有效区分地震与爆炸,尤其对于低震级事件,比高频同频带P/S的分类效果更好。Fisk (20062007)为全球主要核试验场的P/S谱比的频率依赖性进行了一致的理论模型解释,得到P/S谱比的识别性能主要与爆炸P波与S波拐角频率差异以及爆炸的震源谱形状有关。根据地下核爆炸P与S谱随频率变化显著的特征,Fisk等(20082009)建立了地下核爆炸与天然地震的二维网格P/S谱比理论模型,得到核爆炸与天然地震的二维P/S谱比存在差异,表明交叉频带谱比较传统的单一频带谱比有更好的识别能力。

    上节研究表明同频带Pg/Sg谱比能将北京及周边的地震与爆炸完全区分,但是无法对两类爆炸进行识别。对比怀来地区与三河采石场爆炸时域波形记录和主频分析,可以看出这两个区域爆炸有不同的P,S频谱特征。为进一步对两个区域的工业爆炸进行识别,本文开展交叉P/S谱比计算。在0—20 Hz内以2 Hz为间隔分别计算所有频率下的Pg/Sg谱比值,共100个,所有谱比值计算中均进行了距离校正。分析交叉频带谱比值的识别能力,得到与低频Sg (0—2 Hz)有关的比值可有效地对两类爆炸进行分类。Pg (0—2 Hz)/Sg (0—2 Hz),Pg (6—8 Hz)/Sg (0—2 Hz)谱比值分布如图7所示,阈值线所对应的正确识别率分别为92%和94%。

    图  7  交叉频带谱比分布
    Figure  7.  Cross-band spectral ratio values
    (a) Pg (0—2 Hz)/Sg (0—2 Hz);(b) Pg (6—8 Hz)/Sg (0—2 Hz);

    两类爆炸识别中交叉频带谱比正确识别率大于85%的频带列于表3,大部分判据集中在与Sg低频频带(0—2 Hz)有关的比值,其中0—10 Hz的Pg与0—2 Hz的Sg的谱比值识别率可达90%以上,图4b显示三河采石场爆炸S波能量主要集中在1—2 Hz,因此2 Hz内的Pg/Sg谱比值小于怀来爆炸。交叉频带的P/S谱比更能反映出怀来爆炸、三河采石场爆炸不同的P、S波频率分布差异,在性质识别中也显示出较好的分类能力。理论上集中爆炸、延时爆炸等不同爆炸机制产生的P、S谱模型也会存在差异,但实际地震记录还会受到场地条件、埋藏深度、传播路径等影响,因此在未来工作中还需对不同爆炸地震记录的各种影响因素进行校正,才能进一步获得物理机制的明确解释。

    表  3  怀来爆炸与三河采石场爆炸正确识别率大于85%的频带分布
    Table  3.  Frequency bands with the correct discrimination rate greater than 85% between Huailai and Sanhe quarry explosions
    交叉频谱比三河采石场爆
    炸正确识别数
    怀来爆炸正
    确识别数
    正确识
    别率
    Pg (0—2 Hz)/Sg (0—2 Hz) 51 44 92%
    Pg (0—2 Hz)/Sg (4—6 Hz) 46 46 89%
    Pg (0—2 Hz)/Sg (6—8 Hz) 51 40 88%
    Pg (2—4 Hz)/Sg (0—2 Hz) 47 45 89%
    Pg (2—4 Hz)/Sg (4—6 Hz) 53 36 86%
    Pg (4—6 Hz)/Sg (0—2 Hz) 47 47 91%
    Pg (6—8 Hz)/Sg (0—2 Hz) 52 45 94%
    Pg (8—10 Hz)/Sg (0—2 Hz) 50 44 91%
    Pg (10—12 Hz)/Sg (0—2 Hz) 50 39 86%
    Pg (12—14 Hz)/Sg (0—2 Hz) 55 33 85%
    下载: 导出CSV 
    | 显示表格

    基于天然地震、怀来爆炸、三河采石场爆炸的地震记录特征差异,本文利用同频带Pg/Sg谱比判据对天然地震与人工爆炸进行了有效地区分;对于两类不同地区的人工爆炸,通过交叉频带谱比得到与Sg低频带(0—2 Hz)有关的谱比判据可以有效地区分怀来爆炸与三河爆炸;分析同频带谱比、交叉谱比分类能力较好的频带,表明主频分布不仅能反映出不同事件频谱分布的差异,还能对谱比判据的频带选择提供依据。随着更多样本的加入,不同事件的频谱差异界限将更加清晰,进行识别工作时可根据频谱差异直接调整更为适合的P/S谱比频带,有助于提高不同类型事件的性质分类能力。

    河北怀来和三河地区是距离北京市最近的两处爆炸区域,河北怀来山区2002—2015年发生过48次人工爆炸,其中包含2007年用于科学探测研究的大当量爆炸“明灯一号”。河北三河地区发生的主要是采石场爆炸,爆炸数目较多。通过对比怀来爆炸、三河采石场爆炸的地震记录和主频分析,得出怀来爆炸记录频率高于三河采石场爆炸,三河采石场爆炸低频发育。Kim等(1994)指出集中瞬时爆炸相较于延时微差爆炸,其地震记录有更多的高频成分。“明灯一号”爆炸试验方式为集中瞬时爆炸,其它怀来爆炸与“明灯一号”地震记录波形及频谱特征相似,可推测怀来山区爆炸事件大部分为单次集中或者多孔瞬时爆炸的炸山爆炸。

    延时微差爆炸作业方式可能是三河采石场地震记录低频发育的一个因素;怀来地处山区,炸山爆炸一般会有一定埋深,而三河采石场位于平原区,采石场爆炸埋深相对更浅一些,也可能是三河爆炸低频产生的一个因素;首都圈地区地质构造复杂,西北部为山区,东南为凹陷盆地,地震波衰减系数显示西北部山区衰减Q值大,盆地Q值小(朱新运,2011)。怀来地区爆炸事件传播路径多为山区,山区硬岩成分较多,高频衰减慢,而处于平原区沉积层发育的三河地区则高频衰减快,也会造成三河采石场爆炸的低频发育。三河采石场长期有采石活动,大量爆炸会造成场地介质疏松,通过对比历史爆炸地震记录,发现爆炸造成的场地疏松并不是三河采石场低频发育的主要因素,但是场地条件对地震波传播的影响至关重要,因此爆炸方式、爆炸深度、传播路径、场地条件等都有可能是导致这两个区域爆炸地震记录时频域特征差异的因素。

    P/S幅值比是地震与爆炸识别的有效依据。本文对两个不同地区爆炸和天然地震的谱比研究表明,同频带高频Pg/Sg谱比值可以有效地区分地震与爆炸,但是无法对两类爆炸进行识别。针对两类爆炸P,S震相的记录差异,本文计算了交叉频带P/S谱比值,获得了较好的识别结果。P/S幅值比的有效利用,需要处理由于距离、路径和台站效应而引起的显著变化,在P/S幅值比的推广应用中,探索了不同的路径效应校正方法。最常用的是基于一维几何扩散的震中距校正(Rodgers et al,1999),如本文衰减校正所示。但对于介质结构复杂的区域,除了震中距的影响,地壳厚度、沉积层厚度、地形起伏等都会影响区域震相的发育,进而影响震相幅值比的识别效果。Taylor等(2002)提出了震级与距离修正方法(magnitude and distance amplitude correction,缩写为MDAC),即对震源谱进行震源机制、传播路径、几何扩散等校正,使P/S谱更能反映出震源的差异,Walter和Taylor (2001)对MDAC进行了修正,使得搜索参数过程更加灵活,Fisk等(2009)利用MDAC方法校正提高了全球不同核试验场地下核爆炸与地震谱比的识别效能。潘常周等(2007b)采用贝叶斯克里金技术建立了P/S震相幅值比校正曲面,得到在传播路径差异较大的情况下,经过曲面修正后可以降低P/S震相幅值比的离散度。对于小震级事件,尤其是随工业需求设定的爆炸方式,震源函数的模拟以及介质高频响应的区域衰减模型计算都更为复杂,因此还需要不断探索小震级P/S谱比的路径校正方法。

    本文主要对怀来爆炸、三河采石场爆炸以及天然地震三类事件开展了时域地震记录分析、频域多窗口谱分析及Pg/Sg谱比在不同事件性质识别中的应用研究,获得以下结论:

    1) 通过对比相近震中距的地震记录特征得到天然地震整段波形频率较高,S波振幅大于P波振幅;“明灯一号”爆炸P波能量强,衰减快,近台有明显的瑞雷面波,整体S波发育弱;三河采石场爆炸记录无论P波还是S波,低频成分强,S波与面波混淆,持续时间长。

    2) 利用多窗口谱分析技术计算主频成分,得到天然地震频率成分复杂,P,S波主频频带均较宽,P波主要能量集中在2.5—13.2 Hz,S波集中在2.5—9 Hz;爆炸主频频带较窄,怀来爆炸P波主频在2.5—7.2 Hz,S波在2—3.8 Hz,三河采石场爆炸P波主频在1.5—5.8 Hz,S波在1—2 Hz,怀来爆炸P波、S波主频都高于三河采石场爆炸。

    3) 计算天然地震、人工爆炸的Pg/Sg谱比,去除了传播路径中几何扩散的影响,得到高频处(>5 Hz)可将本文地震与爆炸数据完全区分,两个地区爆炸同频带谱比值分布特征相近,无法区分。

    4) 鉴于两类爆炸地震记录和主频分布均存在差异,本文探索了交叉频带谱比对不同区域爆炸的识别,得到与Sg低频(0—2 Hz)有关的谱比判据可进一步区分怀来爆炸及三河采石场爆炸,交叉频带的谱比判据较传统的单一频带能更好地反映出不同类型事件的特征差异。

    本文研究结果显示,不同地区、不同类型爆炸地震记录存在差异,但是影响因素众多,因此对于工业爆炸以及突发爆炸灾害(如天津爆炸、响水爆炸),无论是爆炸机制、还是传播路径等的影响都值得深入研究,下一步工作将继续探索消除传播路径、场地介质等因素的影响,提高谱比判据对事件的识别效能,加深对不同爆炸方式震源机制的认知。

    中国地震局地球物理研究所北京遥测数字地震台网和国家数字测震台网数据备份中心为本研究提供了地震目录和波形数据,作者在此一并表示感谢。

  • 图  7   交叉频带谱比分布

    Figure  7.   Cross-band spectral ratio values

    (a) Pg (0—2 Hz)/Sg (0—2 Hz);(b) Pg (6—8 Hz)/Sg (0—2 Hz);

    图  1   研究所用爆炸、地震及台站位置分布

    Figure  1.   Distribution of  explosions ,seismic events and stations in this article

    图  2   不同类型事件的地震记录波形

    (a)“明灯一号”大当量爆炸;(b)三河采石场爆炸;(c)天然地震。波形上方为事件发生时间及震级,波形右上角为记录台网、台站和震中距,下同

    Figure  2.   Seismic waveforms for different types of events

    (a) “Bright Lamp No.1” large yield explosion;(b) Quarry blast;(c) Earthquake. The time and the magnitude of the events are net on the top of traces,and the station and epicentral distance are on the top right corner of the waveforms,the same below

    图  3   2003 (a),2004 (b)和2005 (b)年3次爆炸事件的波形对比

    Figure  3.   Waveform comparison of three explosion events occured in 2003 (a),2004 (b) and 2005 (c)

    图  4   不同类型事件P波(a)和S波(b)主频分布

    Figure  4.   Main frequency distribution of P (a) and S (b) waves for different types of events

    图  5   天然地震与爆炸Pg/Sg谱比值分布

    Figure  5.   Pg/Sg spectral ratio value of earthquakes and explosions

    图  6   XBZ台站的不同事件地震波形记录(滤波频带为:0.5—25 Hz)

    (a) 2004-03-02怀来爆炸;(b) 2005-03-12可疑事件;(c) 2010-07-10天然地震事件。左上角为事件震级和震中距

    Figure  6.   Seismic waveforms of different events recorded by XBZ (filtered:0.5—25 Hz)

    (a) 2004-03-02 Huailai explosion;(b) 2005-03-12 Suspicious event ; (c) 2010-07-10 natural earthquakeThe magnitude and epicentral distance are on the top left corner of the wavforms

    表  1   天然地震、怀来爆炸、三河爆炸主频分布

    Table  1   The main frequency distribution of earthquake,Huailai explosion and Sanhe explosion

    事件P波主频分布/HzS波主频分布/Hz
    天然地震 2.5—13.22.5—9.0
    怀来爆炸2.5—7.22.0—3.8
    三河采石场爆炸1.5—5.81.0—2.0
    下载: 导出CSV

    表  2   三次对比事件参数

    Table  2   Parameters of three comparison events

    事件性质发震时间 (UTC)震中位置ML深度/km
    年-月-日 时:分:秒北纬/°  东经/°
    怀来爆炸2004-03-02 08:03:35.640.50  115.471.8--
    可疑事件2005-03-12 13:46:27.440.51  115.461.6--
    天然地震2010-07-10 09:53:43.040.47  115.551.87
    下载: 导出CSV

    表  3   怀来爆炸与三河采石场爆炸正确识别率大于85%的频带分布

    Table  3   Frequency bands with the correct discrimination rate greater than 85% between Huailai and Sanhe quarry explosions

    交叉频谱比三河采石场爆
    炸正确识别数
    怀来爆炸正
    确识别数
    正确识
    别率
    Pg (0—2 Hz)/Sg (0—2 Hz) 51 44 92%
    Pg (0—2 Hz)/Sg (4—6 Hz) 46 46 89%
    Pg (0—2 Hz)/Sg (6—8 Hz) 51 40 88%
    Pg (2—4 Hz)/Sg (0—2 Hz) 47 45 89%
    Pg (2—4 Hz)/Sg (4—6 Hz) 53 36 86%
    Pg (4—6 Hz)/Sg (0—2 Hz) 47 47 91%
    Pg (6—8 Hz)/Sg (0—2 Hz) 52 45 94%
    Pg (8—10 Hz)/Sg (0—2 Hz) 50 44 91%
    Pg (10—12 Hz)/Sg (0—2 Hz) 50 39 86%
    Pg (12—14 Hz)/Sg (0—2 Hz) 55 33 85%
    下载: 导出CSV
  • 林大超, 白春华. 2007. 爆炸地震效应[M]. 北京: 地质出版社: 70–78.

    Lin D C, Bai C H. 2007. Explosion Seismic Effect [M]. Beijing: Geological Publishing House: 70–78 (in Chinese).

    潘常周,靳平,王红春. 2007a. P/S震相幅值比判据对低震级地震事件的适用性检验[J]. 地震学报,29(5):521–528.

    Pan C Z,Jin P,Wang H C. 2007a. Applicability of P/S amplitude ratios for the discrimination of low magnitude seismic events[J]. Acta Seismologica Sinica,29(5):521–528 (in Chinese).

    潘常周,靳平,肖卫国. 2007b. 利用克里金技术标定新疆及附近地区P/S震相幅值比及其在地震事件识别中的应用[J]. 地震学报,29(6):625–634.

    Pan C Z,Jin P,Xiao W G. 2007b. Calibration of P/S amplitude ratios for seismic events in Xinjiang and adjacent areas based on a Bayesian Kriging method[J]. Acta Seismologica Sinica,29(6):625–634 (in Chinese).

    孙海霞,杨选,林向东,侯丽娟,鲁跃. 2018. 北京及邻区Q值及台站场地响应分析[J]. 华北地震科学,36(3):1–11. doi: 10.3969/j.issn.1003-1375.2018.03.001

    Sun H X,Yang X,Lin X D,Hou L J,Lu Y. 2018. Analysis on the Q value and station site response in Beijing and its adjacent areas[J]. North China Earthquake Sciences,36(3):1–11 (in Chinese).

    唐兰兰. 2018. 诱发地震检测及地震震源分类[D]. 合肥: 中国科学技术大学: 35–44.

    Tang L L. 2018. Induced Seismicity Detection and Discrimination of Seismic Source[D]. Hefei: University of Science and Technology of China: 35–44 (in Chinese).

    王婷婷,边银菊. 2015. 振幅衰减特性在地震与爆破识别中的应用[J]. 地震学报,37(1):169–179. doi: 10.11939/j.issn:0253-3782.2015.01.015

    Wang T T,Bian Y J. 2015. Amplitude attenuation and its application to earthquake and explosion discrimination[J]. Acta Seismologica Sinica,37(1):169–179 (in Chinese).

    王云专,王珊,董相杰,于承业. 2009. 多窗谱分析在Q值估算中的应用[J]. 地球物理学进展,24(6):2156–2162. doi: 10.3969/j.issn.1004-2903.2009.06.031

    Wang Y Z,Wang S,Dong X J,Yu C Y. 2009. The application of the multiple tapers spectral analysis in the Q estimation[J]. Progress in Geophysics,24(6):2156–2162 (in Chinese).

    徐锡伟, 吴卫民, 张先康, 马胜利, 马文涛, 于贵华, 顾梦林, 江娃利. 2002. 首都圈地区地壳最新构造变动与地震[M]. 北京: 科学出版社: 5–10.

    Xu X W, Wu W M, Zhang X K, Ma S L, Ma W T, Yu G H, Gu M L, Jiang W L.2002. Recent Tectonic Changes and Earthquakes in the Capital Circle Area[M]. Beijing: Science Press: 5–10 (in Chinese).

    赵爱华,郭永霞,孙为国,杨立强,刘哲. 2012. 华北地区矿山爆破活动的时空特征[J]. 地球物理学进展,27(3):917–923. doi: 10.6038/j.issn.1004-2903.2012.03.012

    Zhao A H,Guo Y X,Sun W G,Yang L Q,Liu Z. 2012. Spatio-temporal characteristics of mine blast activity in North China[J]. Progress in Geophysics,27(3):917–923 (in Chinese).

    郑秀芬,傅瑀,许绍燮. 2006. 地震记录中小爆破的识别与判据研究[J]. 地震地磁观测与研究,27(5):29–33. doi: 10.3969/j.issn.1003-3246.2006.05.006

    Zheng X F,Fu Y,Xu S X. 2006. The discrimination and criteria study between blasts and small earthquakes[J]. Seismological and Geomagnetic Observation and Research,27(5):29–33 (in Chinese).

    朱新运. 2011. 衰减、场地响应等地震波传播相关信息综合研究[D]. 北京: 中国地震局地球物理研究所: 17–24.

    Zhu X Y.2011. Comprehensive Study of Seismic Wave Propagation Related Information Such as Attenuation and Site Response[D]. Beijing: Institute of Geophysics, China Earthquake Administration: 17–24 (in Chinese).

    Bennett T J,Murphy J R. 1986. Analysis of seismic discrimination capabilities using regional data from western United States events[J]. Bull Seismol Soc Am,76(4):1069–1086. doi: 10.1785/BSSA0760041069

    Bonner J L,Russell D R,Reinke R E. 2013. Modeling surface waves from aboveground and underground explosions in Alluvium and Limestone[J]. Bull Seismol Soc Am,103(6):2953–2970.

    Chun K Y,Wu Y,Henderson G A. 2011. Magnitude estimation and source discrimination:A close look at the 2006 and 2009 North Korean underground nuclear explosions[J]. Bull Seismol Soc Am,101(3):1315–1329. doi: 10.1785/0120100202

    Douglas A. 2013. Forensic Seismology and Nuclear Test Bans[M]. Cambridge: Cambridge University Press: 340–368.

    Fisk M D,Gray H L,McCartor G D. 1996. Regional event discrimination without transporting thresholds[J]. Bull Seismol Soc Am,86(5):1545–1558.

    Fisk M D, Bottone S, McCartor G D. 2000. Regional seismic event characterization using a Bayesian Kriging approach[C]//Proceedings of the 22nd Annual DOD/DOE Seismic Research Symposium: Planning for Verification of and Compliance with the Comprehensive Nuclear-Test-Ban TreatyCTBT). Louisiana: U.S. Government or Federal Rights: 1–11.

    Fisk M D. 2006. Source spectral modeling of regional P/S discriminants at nuclear test sites in China and the former soviet union[J]. Bull Seismol Soc Am,96(6):2348–2367. doi: 10.1785/0120060023

    Fisk M D. 2007. Corner frequency scaling of regional seismic phases for underground nuclear explosions at the Nevada test site[J]. Bull Seismol Soc Am,97(3):977–988. doi: 10.1785/0120060186

    Fisk M D, Taylor S R, Patton H J, Walter W R. 2008. Applications of a next-generation MDAC discrimination procedure using two-dimensional grids of regional P/S spectral ratios[C]//Proceedings of 2008 Monitoring Research Review: Ground-based Nuclear Explosion Monitoring Technologies: 583–592.

    Fisk M D, Taylor S R, Walter W R, Randall G E. 2009. Seismic event discrimination using two-dimensional grids of regional P/S spectral ratios: Applications to Novaya Zemlya and the Korean Peninsula[C]//Proceedings of the 2009 Monitoring Research Review: Ground-Based Nuclear Explosion Monitoring Technologies: 465–474.

    Hartse H E,Taylor S R,Phillips W S,Randall G E. 1997. A preliminary study of regional seismic discrimination in central Asia with emphasis on western China[J]. Bull Seismol Soc Am,87(3):551–568.

    Kim W Y,Simpson D W,Richards P G. 1993. Discrimination of earthquakes and explosions in the eastern United States using regional high-frequency data[J]. Geophys Res Lett,20(14):1507–1510. doi: 10.1029/93GL01267

    Kim W Y,Simpson D W,Richards P G. 1994. High-frequency spectra of regional phases from earthquakes and chemical explosions[J]. Bull Seismol Soc Am,84(5):1365–1386.

    Kim W Y,Aharonian V,Lerner-Lam A L,Richards P G. 1997. Discrimination of earthquakes and explosions in southern Russia using regional high-frequency three-component data from the IRIS/JSP Caucasus network[J]. Bull Seismol Soc Am,87(3):569–588.

    Kim W Y,Richards P G. 2007. North Korean nuclear test:Seismic discrimination low yield[J]. Eos,Trans Amer Geophys Union,88(14):158–161.

    Kim W Y,Richards P G,Schaff D,Jo E,Ryoo Y. 2018. Identification of seismic events on and near the North Korean test site after the underground nuclear test explosion of 3 September 2017[J]. Seismol Res Lett,89(6):2120–2130.

    Koper K D,Pechmann J C,Burlacu R,Pankow K L,Stein J,Hale J M,Roberson P,McCarter M K. 2016. Magnitude-based discrimination of man-made seismic events from naturally occurring earthquakes in Utah,USA[J]. Geophys Res Lett,43(20):10638–10645. doi: 10.1002/2016GL070742

    McCoy E J,Walden A T,Percival D B. 1998. Multitaper spectral estimation of power law processes[J]. IEEE Trans Signal Process,46(3):655–668. doi: 10.1109/78.661333

    Murphy J R,Bennett T J. 1982. A discrimination analysis of short-period regional seismic data recorded at Tonto Forest Observatory[J]. Bull Seismol Soc Am,72(4):1351–1366. doi: 10.1785/BSSA0720041351

    Pomeroy P W,Best W J,McEvilly T V. 1982. Test ban treaty verification with regional data-a review[J]. Bull Seismol Soc Am,72(6):S89–S129.

    Pasyanos M E. 2018. The coupled location/depth/yield problem for North Korea’s declared nuclear tests[J]. Seismol Res Lett,89(6):2059–2067.

    Smith A T. 1993. Discrimination of explosions from simultaneous mining blasts[J]. Bull Seismol Soc Am,83(1):160–179.

    Rodgers A J, Walter W R, Schultz C A, Myers S C,Lay T. 1999. A comparison of methodologies for representing path effects on regional P/S discriminants [J]. Bull Seismol Soc Am,89(2):394–408.

    Taylor S R,Denny M D,Vergino E S,Glaser R E. 1989. Regional discrimination between NTS explosions and western U.S. earthquakes[J]. Bull Seismol Soc Am,79(4):1142–1176. doi: 10.1785/BSSA0790041142

    Taylor S R,Velasco A A,Hartse H E,Phillips W S,Walter W R,Rodgers A J. 2002. Amplitude corrections for regional seismic discriminants[J]. Pure Appl Geophys,159(4):623–650. doi: 10.1007/s00024-002-8652-8

    Thomson D J. 1982. Spectrum estimation and harmonic analysis[J]. Proc IEEE,70(9):1055–1096. doi: 10.1109/PROC.1982.12433

    Walter W R, Taylor S R. 2001. A Revised Magnitude and Distance Amplitude Correction MDAC2) Procedure for Regional Seismic Discriminants: Theory and Testing at NTS[R]. Livermore, CA: Lawrence Livermore National Lab: 1–13.

    Willis D E,DeNoyer J,Pollack H,Wilson J T. 1963. Differentiation of earthquakes and underground nuclear explosions on the basis of amplitude characteristics[J]. Bull Seismol Soc Am,53(5):979–987. doi: 10.1785/BSSA0530050979

    Xie J K,Patton H J. 1999. Regional phase excitation and propagation in the Lop Nor region of central Asia and implications for P/Lg discriminants[J]. J Geophys Res:Solid Earth,104(B1):941–954. doi: 10.1029/1998JB900045

    Zanandrea A,Da Costa J M,Dutra S L G,Rosa R R,Saotome O. 2004. Spectral and polarization analysis of geomagnetic pulsations data using a multitaper method[J]. Comput Geosci,30(8):797–808. doi: 10.1016/j.cageo.2004.03.016

    Zhao L F,Xie X B,Wang W M,Fan N,Zhao X,Yao Z X. 2017. The 9 September 2016 North Korean underground nuclear test[J]. Bull Seismol Soc Am,107(6):3044–3051. doi: 10.1785/0120160355

    Zheng X F,Yao Z X,Liang J H,Zheng J. 2010. The role played and opportunities provided by IGP DMC of China national seismic network in Wenchuan earthquake disaster relief and researches[J]. Bull Seismol Soc Am,100(5B):2866–2872. doi: 10.1785/0120090257

  • 期刊类型引用(1)

    1. 刘守财. 区域断裂带上的土心墙石渣坝抗震安全分析研究. 湖南水利水电. 2024(03): 90-93+106 . 百度学术

    其他类型引用(0)

图(7)  /  表(3)
计量
  • 文章访问数:  809
  • HTML全文浏览量:  357
  • PDF下载量:  144
  • 被引次数: 1
出版历程
  • 收稿日期:  2020-11-03
  • 修回日期:  2021-04-02
  • 网络出版日期:  2021-07-29
  • 发布日期:  2021-07-14

目录

/

返回文章
返回