京棉二厂井氦、氢气体的地震前兆异常特征及与地震关系的研究.

范树全, 苏盛虎, 李霓

范树全, 苏盛虎, 李霓. 1993: 京棉二厂井氦、氢气体的地震前兆异常特征及与地震关系的研究.. 地震学报, 15(4): 490-497.
引用本文: 范树全, 苏盛虎, 李霓. 1993: 京棉二厂井氦、氢气体的地震前兆异常特征及与地震关系的研究.. 地震学报, 15(4): 490-497.

京棉二厂井氦、氢气体的地震前兆异常特征及与地震关系的研究.

  • 摘要: 系统总结了京棉二厂热水深井溶解气体中氦、氢气体异常变化与地震的关系,提出了氦、氢气体前兆异常特征与发震时间、震级关系的经验公式.用此经验公式较成功地预报了华北两次5级以上地震,具有一定的实用意义.
  • 概率地震需求模型是地震易损性分析和地震风险分析的重要组成部分,它定义了地震动强度表征参数(intensity measure,简写为IM)与结构需求参数(demand measure,简写为DM)之间的关系。IM和DM的不确定性和相关性是建立概率地震需求模型的关键之一。其中IM的选取对概率地震需求模型的建立至关重要。合理的IM不仅能够表征地震动的主要特征,而且可以使给定IM条件下结构反应需求参数的离散性变小。目前对常规地震动强度表征参数的研究较多(Kadas et al,2011李雪红等,2014Kostinakis et al,2015陈健云等,2017王亚楠等,2018),而速度脉冲型地震动与常规地震动存在显著差异,其具有周期长、幅值大、能量瞬时累积大、非平稳性强等特点,对结构有特殊的破坏作用(杨迪雄等,2005Baker,2007陈波等,2013赵晓芬,2015Yang,Zhou,2015),因此有必要针对速度脉冲型地震动强度表征参数展开进一步研究。

    目前针对速度脉冲型地震动强度表征参数的研究已取得了一定进展。Shome等(1998)及Luco和Cornell (2007)的研究表明,基阶模态周期谱加速度SaT)作为IM较其它参数更有效。但速度脉冲型地震动不同周期对应的SaT)之间的相关性随着周期的增大而逐渐减弱,特别是对于长周期结构,SaT)很难有效地解释结构反应的变异性(Baker,Cornell,20042008)。周靖等(2010)提出将SaT)和PGV作为速度脉冲型地震动强度的表征参数,但由于SaT)和PGV未包含谱形状信息(Baker,Cornell,20042008Tothong,Cornell,2008),所以将这两个参数作为IM而估算的结构地震响应的离散性随结构体系周期、延性水平的变化较大。Baker和Cornell (2008)选用SaT)和谱形参数RT1T2)作为速度脉冲型地震动强度的表征参数[其中,RT1T2)=SaT2)/SaT1),T1为一阶模态周期,T2为能够捕捉脉冲记录谱形状的谱值所对应的周期],但由于T2的确定存在很大的不确定因素,因此RT1T2)不适合作为速度脉冲型地震动强度表征参数。Tothong和Cornell (2008)采用SaT)和反应谱偏差参数ε作为速度脉冲型地震动强度表征参数,虽然ε对地震动记录的反应谱形状有重要影响(Abrahamson,1988Baker,Cornell,2006Baker,2011),但是由于ε的确定取决于新一代地震地面运动预测方程(next generation ground motion prediction equations,简写为GMPEs),而现有的经验预测模型未充分考虑速度脉冲的影响,因此,ε不适合作为速度脉冲型地震动的IM。Shahi和Baker (2013)采用速度脉冲对加速度反应谱的放大系数Af表征速度脉冲对反应谱的放大作用,Af反映了脉冲分量对地震动反应谱的放大作用,体现了速度脉冲独特的谱形特征,但以Af作为IM的有效性尚待研究。因此,探寻能够较好地表征速度脉冲谱形状的参数是确定速度脉冲型地震动IM的关键。此外,以往在研究速度脉冲型地震动强度的表征参数时,主要考察地震动本身直接得到的参数如峰值加速度、峰值速度、峰值位移、SaT)等,考察参数不够全面(Tothong,Cornell,2008周靖等,2010),因此应对综合考虑地震动三要素(振幅、频谱、持时)以及速度脉冲特性的多种地震动参数进行全面考察,从而确定速度脉冲强度表征参数及其适用范围。总之,随着速度脉冲记录的不断累积,以及对速度脉冲特性研究的深入,需要进一步研究速度脉冲型地震动IM,为基于性能的工程结构抗震设计和评估提供参考。

    鉴于此,本文拟基于NGA-West2强震动数据库中236条速度脉冲记录,通过研究这些脉冲记录的42个地震动参数,初步确定速度脉冲型地震动IM,并验证IM的有效性,及不同速度脉冲型强度表征参数的适用范围。此外,本文采用Shahi和Baker (2013)提出的速度脉冲对加速度反应谱放大作用系数Af作为IM,并进一步验证其有效性。

    本文从NGA-West2数据库中选取了52个地震事件的236条脉冲记录,基本情况列于表1

    表  1  本文所用236条速度脉冲记录的基本信息
    Table  1.  The information of 236 pulse-like velocity recordings used in this paper
    震级范围NRrup/kmVp/(cm·s−1Tp/svS30/(m·s−1
    5.0≤MW<6.0163—1723—930.2—4.4190—665
    6.0≤MW<7.01380.1—7824—1550.3—6.9139—2 017
    7.0≤MW<8.0820.3—9327—3420.8—13.2141—1 370
      注:表中N表示速度脉冲记录个数,Rrup表示断层距,Vp表示速度脉冲幅值,Tp表示脉冲周期,vS30表示场地30 m平均剪切波速。
    下载: 导出CSV 
    | 显示表格

    国内外研究人员提出的42种地震动参数列于表2,其中脉冲周期、脉冲幅值、脉冲因子、线性组合系数、最大正负速度峰值差、脉冲循环数及归一化累计平方速度差是描述速度脉冲特征的参数,其它参数是描述常规地震动特性的参数。

    表  2  42种地震动参数
    Table  2.  42 ground motion parameters
    编号名称标识符编号名称标识符编号名称标识符
    1脉冲周期Tp15速度反应谱均值Sv,avg29阿里亚斯强度IA
    2脉冲幅值Vp16位移反应谱均值Sd,avg30修正的阿里亚斯强度Imia
    3脉冲因子PI17有效峰值加速度EPA31Faifar指标If
    4线性组合系数PC18有效峰值速度EPV32累积绝对速度CAV
    5最大正负速度峰值差PPV19有效峰值位移EPD33累积绝对位移CAD
    6脉冲循环数Ncycles20豪斯纳强度SI34累积绝对动力CAI
    7归一化累计平方速度差NCSVdiff21峰值速度/峰值加速度V/A35Nau和Hall指标Vrs
    8地震动峰值加速度PGA22Mackie地震动谱强度ASI36Nau和Hall指标Drs
    9地震动峰值速度PGV23Mackie地震动谱强度VSI37地震动均方根强度指标RMSa
    10地震动峰值位移PGD24Mackie地震动谱强度DSI38地震动均方根强度指标RMSv
    11加速度反应谱峰值Sa,max25峰值位移/峰值速度D/V39地震动均方根强度指标RMSd
    12速度反应谱峰值Sv,max26加速度反应谱卓越周期Tpa40Park-Ang指标Ic
    13位移反应谱峰值Sd,max27特征周期Tg41最大增量速度MIV
    14加速度反应谱均值Sa,avg28括号持时Td42最大增量位移MID
     注:表中42种地震动参数的计算方法及意义详见杨迪雄等(2005)叶列平等(2009)Haselton等(2012)陈波(2013).
    下载: 导出CSV 
    | 显示表格

    速度脉冲对加速度反应谱的放大表现为一种窄带效应,即速度脉冲并不是对地震动所有频带一致放大,而是仅仅放大了部分频段。Shahi和Baker (2013)通过系数Af来定量地考察脉冲对加速度反应谱的放大作用,将该系数定义为

    $ \ln {{S}_{\rm{a,pulse}}}{\text{=}}\ln (\frac{{{S}_{\rm{a,pulse}}}}{{{S}_{\rm{a,res}}}}\cdot {{S}_{\rm{a,res}}}){\text{=}}\ln({{A}_{\rm{f}}}\cdot {{S}_{\rm{a,res}}}){\text{=}}\ln{{A}_{\rm{f}}}{\text{+}}\ln{{S}_{\rm{a,res}}}{\text{,}} $

    (1)

    ${A_{\rm{f}}} {\text{=}} \frac{{{S_{{\rm{a,pulse}}}}}}{{{S_{{\rm{a,res}}}}}}{\text{,}}$

    (2)

    式中,Sa,pulse为脉冲记录的加速度反应谱,Sa,res表示提取脉冲信号后残余记录的加速度反应谱。 此外,Af还体现了速度脉冲独特的谱形特征,表征了速度脉冲的主要特性。

    本研究采用具有不同延性水平的弹塑性单自由度体系(single-degree-of-freedom,简写为SDOF)分析模型,其中,分析模型的本构模型采用双线性恢复力模型,如图1所示。fouo为弹性体系的最大力和最大变形,fmum为弹塑性体系的最大力和最大变形,fyuy为弹性体系的屈服力和屈服位移,k为弹性刚度,阻尼比ξ为0.05,刚度折减系数α为0.1。选取结构自振周期0.1—6.0 s 范围内的50个周期值以及延性系数μ分别取为1,2,3,4,6,8的300个弹塑性单自由度分析模型,开展结构地震反应分析,系统最大非线性位移um作为结构地震响应需求参数。

    图  1  双线性本构模型
    Figure  1.  Bilinear constitutive model

    具有代表性是地震动强度表征参数的重要特征之一,即与其它地震动参数具有较高的相似性和相关性。相似性和相关性程度通常分别采用距离分析和秩相关性、皮尔森(Pearson)相关性进行度量。

    统计学中变量之间的相似程度通过距离分析来进行度量。距离分析是研究地震动参数指标相关关系的基础,常用的距离分析包括欧式(Euclid)距离、切比雪夫(Chebyshev)距离及布洛克(Block)距离(宋帅等,2016)。本文选取欧式距离进行度量,其表达式为

    $d(x{\text{,}}y) {\text{=}} \sqrt {\sum\limits_{i {\text{=}} 1}^n {({x_i}} {\text{-}} {y_i}{)^2}} {\text{,}}$

    (1)

    式中:xy为任意两个参数,xiyi为参数的样本值,n为样本数量。参数间的距离越小,相似性则越好。为了消除量纲不同所造成的影响,需对地震动参数指标进行如下标准化处理:

    $x_i' {\text{=}} \frac{{{x_i}{\text{-}} \overline{x} }}{\sigma }{\text{,}}$

    (2)

    This page contains the following errors:

    error on line 1 at column 1: Start tag expected, '<' not found

    Below is a rendering of the page up to the first error.

    图  2  地震动参数间距离之和的平均值
    Figure  2.  The average of the sum of distance between ground motion parameters

    秩相关性分析(宋帅等,2016)用于度量参数之间的相关性程度,其既可描述参数间的线性相关,又可描述参数间的非线性相关。设(xiyi)为二元总体的样本,秩相关系数定义为

    ${r_{\rm{s}}}(x{\text{,}}y) {\text{=}} 1 {\text{-}} \frac{6}{{n({n^2} {\text{-}} 1)}}{\sum\limits_{i {\text{=}} 1}^n {({\gamma _i} {\text{-}} {\theta _i})} ^2}{\text{,}}$

    (3)

    式中,γiθi分别表示样本值xiyi的秩观测值。通过秩相关性分析,将每个地震动参数与其它参数间的秩相关系数取绝对值后再求和,由大到小排列,结果如图3所示。可以看出,Sa,avgSv,avg,PPV,If,PGV与其它参数间的相关性最高。

    图  3  地震动参数的秩相关系数的绝对值之和
    Figure  3.  The sum of absolute values of rank correlation coefficients of ground motion parameters

    基于秩相关分析结果进一步通过皮尔森相关分析探讨地震动参数之间的线性相关性(叶列平等,2009周靖等,2010陈波,2013)。用ρxy)定义xy之间的皮尔森相关性

    $ \rho (x{\text{,}}y) {\text{=}} \frac{{{\rm {COV}}(x{\text{,}}y)}}{{\sqrt {D(x)} \sqrt {D(y)} }} {\text{=}} \frac{{\displaystyle\sum\limits_{{{i}} {\text{=}} 1}^n {({x_i} {\text{-}} {\bar x})(\,{y_i} {\text{-}} {\bar y})} }}{{\sqrt {\displaystyle\sum\limits_{i {\text{=}} 1}^n {{{({x_i} {\text{-}} {\bar x})}^2}} } \sqrt {\displaystyle\sum\limits_{i {\text{=}} 1}^n {{{({\,y_i} {\text{-}}{\bar y})}^2}} } }}{\text{,}} $

    (4)

    式中,COV(xy)为xy的协方差,Dx)和Dy)为参数的方差。当ρ=1时,xy完全正相关;当ρ=−1时,xy完全负相关。

    对于皮尔森线性相关性分析结果,将每个地震动参数与其它参数间的皮尔森相关系数的绝对值求和,由大到小排列,结果如图4所示。可以看出,Sa,avgSv,avg,PPV,If,PGV与其它参数间的皮尔森相关性最高。

    图  4  地震动参数间皮尔森相关系数绝对值和
    Figure  4.  The sum of absolute values of Pearson correlation coefficients of ground motion parameters

    综上所述,由图234可知,Sa,avgSv,avg,PPV,If,PGV与其它参数间的相似性、秩相关性和皮尔森相关性均较高。对于速度脉冲型地震动而言,传统的IM (如PGA)与其它参数的相似性、秩相关性、皮尔森相关性较低;PGV不是除SaT)外唯一可以作为脉冲型地震动IM的参数。此外,速度脉冲型地震动特有参数中的Ncycles,NCSVdiff,PC,Tp分别与其它地震动参数的相似性较低、秩相关性和皮尔森相关性也较低,这些参数虽可体现速度脉冲的主要特点,但是无法作为速度脉冲型地震动IM。因此,本研究中初步选取Sa,avgSv,avg,PPV,If,PGV作为速度脉冲型地震动IM,并对其展开进一步研究。

    为了验证Sa,avgSv,avg,PPV,If,PGV作为速度脉冲型地震动IM的条件下,结构反应需求参数的离散性程度,本文采用1.4节中的弹塑性单自由度结构模型,输入236条速度脉冲记录进行非线性时程分析。采用皮尔森相关系数(周靖等,2010宋帅等,2016)度量Sa,avgSv,avg,PPV,If,PGV及SaT)与um之间的相关程度,结果如图5所示。

    图  5  不同延性系数水平时IM与um的相关性变化趋势
    Figure  5.  The variation trend of the correlation coefficients between IM and um under different ductility coefficients levels

    图5可以看出:当延性系数恒定时,Sa,avgSv,avg,PPV,If,PGV与um的相关系数随着周期的增大呈先增大再减小的变化趋势;若延性系数发生变化,当0<T<2 s时,相关系数随着延性系数的增大而增大,在T=2 s附近达到最大,当T>2 s时,相关系数随着延性系数的增大呈减小趋势。总之,SaT)较Sa,avgSv,avg,PPV,If,PGV与um之间的相关性更好,且在0—6 s范围内的相关性变化趋势最平稳。

    上述分析表明,SaT),Sa,avgSv,avg,PPV,If,PGV这几个IM参数与上述弹塑性单自由度结构模型的最大非线性位移um均具有较好的相关性,但是以其作为IM的有效性还有待于进一步研究。IM有效性定义为:在给定地震动IM水平下,结构的响应需求离差较小或比较稳定,但并无具体标准(王建民,朱晞,2006)。通常以给定地震动强度水平下结构地震反应的离散性程度来分析地震动强度参数的有效性(Luco,Cornell,2007周靖等,2010)。离散性程度由结构响应的概率密度函数确定,即假定DM|IM是对数正态分布,结构响应的概率密度函数表达为

    $f({\rm{DM}}|{\rm{IM}}) {\text{=}} \frac{1}{{\sqrt {2\pi } {\sigma _{\ln ({\rm{DM}}|{\rm{IM}})}}}} \exp \left\{ - \frac{1}{2}\left[\frac{{\ln ({\rm{DM}}) {\text{-}} \ln({\hat{\rm {DM}}})}}{{{\sigma _{\ln ({\rm{DM}}|{\rm{IM}})}}}}\right]\right\} {\text{,}}$

    (5)

    位移需求的对数平均和标准偏差分别为

    $\ln ({\rm{DM}}) {\text{=}} \ln{a} {\text{+}} {b} \ln({\rm{IM}}){\text{,}}$

    (6)

    ${\sigma _{{\rm{DM|IM}}}} {\text{=}} \sqrt {\frac{{\displaystyle\sum\limits_{i {\text{=}} 1}^n {{{\left[\ln({\rm{DM}}) {\text{-}} \ln{a}\cdot ({\rm{I}}{{\rm{M}}^{{b}}})\right]}^2}} }}{{n - 2}}} {\text{,}}$

    (7)

    式中,ab表示回归参数,σDM表示这组回归观测值的离散程度。对不同自振周期及不同延性系数的结构的地震反应参数进行统计分析,结构地震反应的离散性如图6所示,即不同延性系数水平时,选取Sa,avgSv,avg,PPV,If,PGV作为IM的um的对数标准偏差。

    图  6  不同延性系数水平时选取Sa,avgSv,avg,PPV,If,PGV作为IM的um的对数标准偏差
    Figure  6.  The logarithmic standard deviation of um under different ductility coefficients levels when Sa,avgSv,avg,PPV,If,and PGV are taken as IM,respectively

    图6可知,在延性系数恒定的情况下,选取Sa,avgSv,avg,PPV,If,PGV作为速度脉冲型IM时,结构地震反应um的标准偏差差异不明显,均随周期的增大呈先减小再增大的趋势,且在1<T<2 s时最小。该条件下,当0<T<1 s时,与Sa,avgSv,avg,PPV,If,PGV相比,选取SaT)作为IM参数,um的标准偏差较小;当1<T<3 s时,与SaT)相比,选取Sa,avgSv,avg,PPV,If,PGV作为IM,um的标准偏差较小;当3<T<6 s时,选取Sa,avgSv,avg,PPV,If,PGV以及SaT)作为IM,um的标准偏差均较大。当延性系数变化时:在0<T<2 s范围内,选取Sa,avgSv,avg,PPV,If,PGV作为IM时um的标准偏差随着延性系数的增大而减小;当2<T<6 s时,标准偏差随着延性系数的增大而增大;在0<T<6 s范围内,选取SaT)作为IM,um的标准偏差分别随周期、延性系数的增大而增大。综上,在延性系数恒定的情形下:当0<T<1 s时,选取SaT)作为IM的有效性较好;当1<T<3 s时,选取Sa,avgSv,avg,PPV,If,PGV作为IM的有效性较好;当3<T<6 s时,选取Sa,avgSv,avg,PPV,If,PGV及SaT)作为速度脉冲型地震动IM的有效性较差。

    上述研究表明Sa,avgSv,avg,PPV,If,PGV和SaT)并不适宜在全周期段(0—6 s)内作为速度脉冲型地震动强度表征参数。不同延性水平条件下,选取Af作为IM时,结构地震反应参数um的标准差随周期及周期与脉冲周期之比的变化趋势分别如图7a图7b所示。由图7a可知,当Af作为IM时,在延性系数恒定的条件下,um的离散性呈现随着周期增大而增大的趋势;在延性系数变化的条件下,当T<2 s时,不同延性系数水平下um的离散性相差不大,当T>2 s时,um的离散性随延性系数的增大而增大。对比图7a图6可知,当0<T<6 s时,与Sa,avgSv,avg,PPV,If,PGV,SaT)相比,选取Af作为IM,um的离散性总体上较小。因此AfSa,avgSv,avg,PPV,If,PGV,SaT)更适合作为速度脉冲型地震动强度表征参数。此外由图7b可知:在延性系数变化的情形下,当T/Tp≤1时,不同延性水平条件下,um的离散性相差不大,当T/Tp>1时,um的离散性呈现随延性系数增大而增大的趋势;当延性系数恒定时,um的离散性呈现随T/Tp增大而增大的趋势,且um的离散性在T/Tp≤1时较T/Tp>1时小,即当结构的自振周期小于速度脉冲周期时,选取Af作为IM的有效性最强。

    图  7  不同延性系数时标准差随周期(a)和周期与脉冲周期之比的变化趋势(b)
    Figure  7.  The variation of the standard deviation with period (a) and the ratio of the period to the pulse period (b) for different ductility coefficient

    以236条速度脉冲型地震动记录作为输入,分析了42个地震动参数之间的相似性和相关性,初步确定了速度脉冲型地震动强度的表征参数,经分析评估采用速度脉冲对加速度反应谱的放大作用系数Af作为速度脉冲型地震动强度表征参数,并对其有效性进行了深入分析,主要结论如下:

    1) Sa,avgSv,avg,PPV,If,PGV与其它参数的相似性、相关性较好;而速度脉冲型地震动特有的参数如Ncycles,NCSVdiff,PC,Tp 与其它地震动参数的相似性较低。

    2) Sa,avgSv,avg,PPV,If和PGV与um之间的相关性受周期和延性系数的影响较大。与Sa,avgSv,avg,PPV,If,PGV相比,当周期处于0—6 s范围内时,SaT)与um之间的相关性较好,受周期和延性系数的影响较小,即相关性变化趋势最平稳。

    3) 不同参数作为速度脉冲型IM的有效性的适用范围不同。当0<T<1 s时,SaT)作为速度脉冲型地震动表征强度参数的有效性较好;当1<T<3 s时,Sa,avgSv,avg,PPV,If,PGV作为速度脉冲型地震动强度表征参数的有效性较好;当3<T<6 s时,Sa,avgSv,avg,PPV,If,PGV以及SaT)作为速度脉冲型地震动IM的情形下,随着结构自振周期以及延性系数的增大,um的离散性呈增大趋势。因此,Sa,avgSv,avg,PPV,If,PGV及SaT)不适宜在整个研究的周期范围内作为速度脉冲型地震动强度表征参数。

    4) 当0<T<6 s时,与Sa,avgSv,avg,PPV,If,PGV,SaT)相比,选取Af作为IM的情形下,um的离散性较小,因此Af是相对较优的速度脉冲型地震动IM。此外,当结构自振周期小于速度脉冲周期时,Af作为IM的情形下,um的离散性最小,其有效性最强。当T/Tp≤1时,不同延性水平条件下,um的离散性相差不大;当T/Tp>1时,um的离散性呈现随延性系数增大而增大的趋势,um的离散性受延性系数的影响较大。

    综上所述,与以往的研究相比,本研究通过采用较充分的脉冲记录,从振幅、频谱、持时及速度脉冲特性等多方面考察,从而确定了适用于较大周期范围的速度脉冲型地震动强度表征参数。对于速度脉冲型地震动记录,Af可在较大周期范围内作为速度脉冲型强度表征参数,但是,对于Af如何在速度脉冲型地震动记录输入中应用,尚需进一步研究。

    太平洋地震工程研究中心(Pacific Earthquake Engineering Research Center,简写为PEER)为本文提供了速度脉冲记录,斯坦福大学Jack Wesley Baker教授和Shrey Kumar Shahi博士为本文提供了多分量速度脉冲识别方法的程序,作者在此一并表示感谢。

  • [1] 北京市地质局水文一大队,1973.北京地区的热矿水.水文地质I:程地质.1:1——3.

    [2] 国家地震局.1985.地震水文地球化学观测技术规范.16——17.地震版社.北京,

    [3] 郑炳华、敏顺民、徐好民,1981.燕山地区北西向和北西西向断裂构.造获本特征初步探讨.地震地质.3.2.31——36

    [4] QrKOnur, B. A.,1971. Геохмия Прпрогиых Газов,,290——291. Нздателъетво Недра,Моекна.

    [1] 北京市地质局水文一大队,1973.北京地区的热矿水.水文地质I:程地质.1:1——3.

    [2] 国家地震局.1985.地震水文地球化学观测技术规范.16——17.地震版社.北京,

    [3] 郑炳华、敏顺民、徐好民,1981.燕山地区北西向和北西西向断裂构.造获本特征初步探讨.地震地质.3.2.31——36

    [4] QrKOnur, B. A.,1971. Геохмия Прпрогиых Газов,,290——291. Нздателъетво Недра,Моекна.

  • 期刊类型引用(1)

    1. 刘长生,刘学谦,周晨,李娇,高峰. 绥化地震台视电阻率数据跟踪分析. 地震地磁观测与研究. 2024(06): 72-82 . 百度学术

    其他类型引用(0)

计量
  • 文章访问数:  1217
  • HTML全文浏览量:  10
  • PDF下载量:  124
  • 被引次数: 1
出版历程
  • 发布日期:  2011-09-02

目录

/

返回文章
返回