2014年1月3日，北美局部地区发生显著的大地震颤现象，并伴随有巨大轰爆声响。该现象主要发生在加拿大多伦多及其附近地区，自瑞克斯代尔(Rexdale)、贝尔维尔(Belleville)、朔姆贝格(Schomberg)、列治文山(Richmond Hill)、新市(Newmarket)，至克雷茵堡(Kleinburg)、斯多夫(Stouffville)、阿贾克斯(Ajax)、奥沙瓦(Oshawa)、宾顿(Brampton)和奥兰治(Orangeville)。当地多数居民被轰爆声响惊醒，并在描述这些声响时将其比作“两辆火车相撞”或“屋内大梁崩断”。室外的行人也反映听到类似“巨大的枪炮声”的轰爆声响，并感受到大地震颤。这些类似破坏性大地震前兆的现象造成了部分居民的恐慌，经新闻报道也引起了政府管理人员及相关科研人员的担忧。虽然这一伴随巨大轰爆声响的大地震颤现象与地震现象相似，但却并不是传统定义的地震，而是一种新型的自然灾害，称为冻震(frostquake)或称为霜冻地震(cryoseism)^[1]。

冻震是由地表结构霜冻作用激发的大地震颤现象，其中霜冻作用即指地表结构中的土壤或岩石中的水分冻结成冰后体积膨胀挤压土体或岩层的作用^[1-6]。冻震与冰震(icequake)不同，冰震指的是冰川和冰层在冰流滑动的激励下或在地震波穿过冰层的影响下产生的震动，主要发生在极地地区^[7-9]。因此，2014年1月多伦多大地震颤不是冰震，而是冻震。冻震与传统定义的构造型地震的区别在于成因不同，冻震的成因是地表结构中土-水材料或岩-水材料的剧烈霜冻作用，而非地壳构造变动，因此冻震也可被划分归属于非构造型地震动一类，其他如塌陷地震动(collapse-earthquake)、水库地震动(reservoir-induced earthquake)等。但冻震现象却与传统地震现象相似，如前所述，尤其类似于破坏性大地震前兆的现象。根据修正Mercalli烈度(modified Mercalli intensity，MMI)表的定义，2014年1月多伦多冻震依照其新闻报道及居民震感描述可以被评定为IMM达Ⅴ度，有较强震感、无明显震害。

此次多伦多冻震并非史无前例，根据目前的文献记载，人类纪录的首次冻震现象可以追溯到1819年美国Massachusetts Derfield River Region冻震^[1]。虽然冻震不属于传统定义的地震范畴，但历史上多次出现误将冻震记载为地震的情况^[2-5]，这一现象说明剧烈霜冻作用诱发的强烈冻震具有较强地震的特征，如震感、甚至震害，因此也引起了相关研究人员的关注与担忧。Lacroix^[1]在1980年依据当时已有的文献报道，首次对冻震现象进行了较全面系统的定义与归纳总结，提出了形成冻震的地质条件、气象条件以及冻震的震颤特征。Barosh^[5]在2000年通过对New England地区自1958年至1999年间发生的冻震事件的整理汇总，首次给出了冻震作用导致房屋结构开裂的实证照片，以及冻震引起的地表结构断裂的纵剖照片。目前文献中对于冻震的研究多为现象描述、震害报道以及归纳统计，尚未针对其形成机理提出相关的理论模型，无法为定量估测冻震作用可能引起的地面震颤的峰值加速度水平提供理论基础，也难以进一步预测其可能达到的最强烈度等级与可能造成的震害大小。

本文旨在估测冻震作用可能引起的地面峰值加速度(peak ground acceleration，PGA)的最大值，预测冻震可能达到的最大烈度，并分析评估冻震造成的最高震害水平。

1 冻震的形成

图 1所示为2014年1月北美地区冻震城市郊区情况，图 2所示为其中多伦多地区路面、地面、山体岩石断裂产生裂缝的情况。由图 1可以看出冻震发生时的两点气象特征：1) 大规模降水与降温引起霜冻作用；2) 地表少量或无积雪覆盖。图 2证明土壤、岩石缝隙中水分的霜冻作用会导致地表结构断裂，冻震过程中巨大的声响则是在断裂过程中发生的，彻夜持续的轰爆声响则说明同一地区多处发生地表结构断裂情况，且各处断裂在一个较长的时间区间内先后发生。

图 3所示则为文[5]纪录的1979年2月美国Massachusetts Pepperell冻震造成的民用房屋结构室内楼板开裂情况，表明冻震作用可能引起房屋结构产生震害。表 1给出了历史上被误录为地震的部分较强冻震事件及其对应的MMI值^[2-5]，其中有2次冻震的MMI达到Ⅵ度，分别为1855年2月Hudson River Valley冻震与1952年1月Burlington冻震，表 1中其他冻震的MMI值在Ⅳ~Ⅴ度之间，表明冻震的MMI可达Ⅳ~Ⅵ度。表 2所示为经本文整理的Lacroix^[1]在1980年给出的冻震三特征，包括气象条件、地质条件及震颤特征，可为下一步研究提供支持与验证。

2 冻震的模型

图 4所示为一种可能的冻震形成机制, 主要可分为4步：1) 大规模降水，原初地质条件为表层结构有裂纹或裂缝，对应表 2中的2a项，气象条件则为短时间内的大规模降水，对应表 2中的1c项；2) 快速渗透、大量蓄水，其所需的地质条件为地表结构材料易渗透，对应表 2中的2b项，气象条件为大规模降水导致地表结构材料迅速饱和，仍对应表 2中的1c项；3) 大规模降温、发生霜冻，所需条件主要为气象条件，包括霜冻气候、微少或无积雪、快速大降温，分别对应表 2中的1a、1b、1d项；4) 形成冻震，其过程可以概括为挤压变形引起的能量蓄积与结构断裂伴随的能量释放，包括地表结构中的土壤及岩层，表现出的现象主要为多处产生断裂、断裂伴随轰爆声响、断裂分时发生、总持续时间较长，整体形成大地震颤。另外，由图 4可以看出，冻震的成因与形成位于地表结构，地表结构的断裂与震颤也位于该处，因此冻震的震源与震场基本位于同一空间位置，没有明显的传播现象。这与传统的构造型地震显著不同。

2.1 最大峰值速度

依据材料应变能完全转换为大地运动动能的假定得：

$ {{G}_{\text{m-fault}}}=\frac{1}{2}M{{(V_{\text{p}}^{*})}^{2}}. $

(1)

其中：G_m-fault为应变能，M为质量，V_p^*为名义峰值速度。再将弹性变形假定代入式(1) 得：

$ \begin{align} &\frac{\frac{1}{2}\sigma \varepsilon {{V}_{\text{total}}}}{M}=\frac{1}{2}{{\left( V_{\text{p}}^{*} \right)}^{2}} \\ &\ \ \ \ \ \Rightarrow \frac{\sigma \varepsilon }{\rho }={{\left( V_{\text{p}}^{*} \right)}^{2}} \\ &\ \ \ \ \Rightarrow \frac{E{{\varepsilon }^{2}}}{\rho }={{\left( V_{\text{p}}^{*} \right)}^{2}}. \\ \end{align} $

(2)

其中：σ为霜冻作用应力，ε为霜冻作用应变，V_total为土-水、岩-水或岩-土-水材料总体积，ρ为材料平均密度，E为材料弹性模量。由式(2) 解得V_p^*为

$ V_{\text{p}}^{*}=\varepsilon \sqrt{\frac{E}{\rho }}. $

(3)

2.2 霜冻作用应变

霜冻作用应变主要由材料中的水分凝固结冰体积膨胀受到约束产生。首先假设膨胀不受约束，根据体积叠加原理，可得：

$ \begin{matrix} {{V}_{\text{total}}}={{V}_{\text{soil}}}+{{V}_{\text{water}}} \\ \Rightarrow \text{d}{{V}_{\text{total}}}=\text{d}{{V}_{\text{water}}}. \\ \end{matrix} $

(4)

其中：V_soil为土、岩或土-岩的体积，V_water为水的体积，dV_total为总体积变化量，dV_water为水的体积变化量即水凝固结冰的体积膨胀量。由于受到约束作用，真实的材料总体积不变，而体积膨胀趋势引起材料压应变，即霜冻作用应变ε，由于水凝固结冰体积膨胀25%，ε可按下式计算：

$ \varepsilon = \frac{1}{3}\frac{{{\rm{d}}{V_{{\rm{total}}}}}}{{{V_{{\rm{total}}}}}} = \frac{1}{3}\frac{{{\rm{d}}{V_{{\rm{water}}}}}}{{{V_{{\rm{total}}}}}} = \frac{1}{3}\frac{{0.25{V_{{\rm{water}}}}}}{{{V_{{\rm{total}}}}}} \Rightarrow \varepsilon = \frac{1}{{12}}n $

(5)

其中：n为材料初始体积含水率，n=V_water/V_total。式(3) 表明V_p^*与ε成正比例关系，式(5) 表明ε与n成正比例关系，所以V_p^*与n也成正比例关系，即含水率越高，V_p^*越大，冻震颤振现象越强。这与表 2“冻震的形成条件与震颤特征”中1c项及2b项^[1]条件对高渗水率及高含水量的需求一致。

2.3 震颤函数

用3种典型函数描述震颤在一个周期内的地面运动，用一系列周期性重复的此种函数近似冻震的地面运动时程。这3种微断层震颤函数分别为信号函数(式(6)，图 5)、简谐函数(式(7)，图 6)与复合函数(式(8)，图 7)。通过式(6)—(8)，可以将地面峰值加速度PGA与有效峰值速度V_p联系起来，从而计算PGA。

信号函数：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;v = \left\{ \begin{array}{l} {V_{\rm{p}}}-\frac{{4{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}t, \;\;\;\;\;\;\;\;\;\;0 \le t < \frac{{{T_{\rm{p}}}}}{2};\\ -3{V_{\rm{p}}} + \frac{{4{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}t, \;\;\;\;\;\;\;\;\;\;\frac{{{T_{\rm{p}}}}}{2} \le t < {T_{\rm{p}}}\; \end{array} \right.\\ \Rightarrow u = \left\{ \begin{array}{l} \frac{1}{8}{V_{\rm{p}}}{T_{\rm{p}}}-\frac{{2{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}{\left( {t - \frac{{{T_{\rm{p}}}}}{4}} \right)^2}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le t < \frac{{{T_{\rm{p}}}}}{2};\\ - \frac{1}{8}{V_{\rm{p}}}{T_{\rm{p}}} + \frac{{2{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}{\left( {t - \frac{{3{T_{\rm{p}}}}}{4}} \right)^2}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{{T_{\rm{p}}}}}{2} \le t < {T_{\rm{p}}} \end{array} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \Rightarrow a = \left\{ \begin{array}{l} - \frac{{4{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}, \;\;0 \le t < \frac{{{T_{\rm{p}}}}}{2};\\ \frac{{4{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}, \;\;\;\;\;\frac{{{T_{\rm{p}}}}}{2} \le t < {T_{\rm{p}}}. \end{array} \right. \end{array} $

(6)

其中：u为震颤位移，v为震颤速度，a为震颤加速度，T_p是震颤周期。因为冻震的震场也是其震源所在，所以近似取场地特征周期为T_p的取值，此处取四类场地值T_p=0.9 s。

简谐函数：

$ \begin{array}{*{20}{c}} {v = {V_{\rm{p}}}\cos \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}t} \right), 0 \le t < {T_{\rm{p}}}}\\ { \Rightarrow u = \frac{{{V_{\rm{p}}}{T_{\rm{p}}}}}{{2{\rm{\pi }}}}\sin \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}t} \right), 0 \le t < {T_{\rm{p}}}}\\ { \Rightarrow a = \frac{{2{\rm{\pi }}{V_{\rm{p}}}}}{{{T_{\rm{p}}}}}\sin \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}t} \right), 0 \le t < {T_{\rm{p}}}.} \end{array} $

(7)

复合函数：

$ \begin{array}{*{20}{c}} {v = \left\{ \begin{array}{l} {V_{\rm{p}}}{{\rm{e}}^{- {n_1}({n_3}{T_{\rm{p}}}- t + {t_0})}}\sin \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}(t- {t_0})} \right), 0 \le t < {n_3}{T_{\rm{p}}};\\ {V_{\rm{p}}}{{\rm{e}}^{ - {n_2}(t - {t_0} - {n_3}{T_{\rm{p}}})}}\sin \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}(t - {t_0})} \right), {n_3}{T_{\rm{p}}} \le t < 2{T_{\rm{p}}} \end{array} \right.}\\ { \Rightarrow a = \left\{ \begin{array}{l} {V_{\rm{p}}}{{\rm{e}}^{ - {n_1}({n_3}{T_{\rm{p}}} - t + {t_0})}}\left[{{n_1}\sin \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}(t-{t_0})} \right) + \frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}\cos \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}(t-{t_0})} \right)} \right], 0 \le t < {n_3}{T_{\rm{p}}};\\ {V_{\rm{p}}}{{\rm{e}}^{ - {n_2}(t - {t_0} - {n_3}{T_{\rm{p}}})}}\left[{-{n_2}\sin \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}(t-{t_0})} \right) + \frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}\cos \left( {\frac{{2{\rm{\pi }}}}{{{T_{\rm{p}}}}}(t-{t_0})} \right)} \right], {n_3}{T_{\rm{p}}} \le t < 2{T_{\rm{p}}}. \end{array} \right.} \end{array} $

(8)

其中：n₁、n₂、n₃为曲线参数，t₀为初始相位参数。本文参考Menun与Fu^[10]在研究近场地震动给出的建议取值(n₁=0.5，n₂=0.5，n₃=0.75，t₀=0.25T_p)。

3 冻震-地震震害比较研究

本节选取了一系列参数取值来估算冻震作用下最大可达到的V_p范围(n=0.05~0.2, E=496~1 428 MPa, ρ=2 000 kg/m³)，由式(3) 和(5) 计算得到冻震作用下最大可达到的V_p范围是0.19~1.41 cm/s。由式(6)—(8) 可以计算得到冻震作用下最大可达的PGA范围是0.84~6.26 gal(信号函数，1 gal≈4.5 L)、1.33~9.83 gal(简谐函数)、2.33~17.26 gal(复合函数)。根据已有研究给出的PGA-MMI相关关系公式(表 3)，将3种函数得到的最大PGA分别绘于图 8所示的MMI-PGA相关关系图中，可以看出，预测冻震最大PGA范围0.84~6.26 gal(信号函数)、1.33~9.83 gal(简谐函数)、2.33~17.26 gal(复合函数)对应的MMI均不超过Ⅵ度，虽然不同的函数选择会影响预测的结果，但其影响不超过1度，即采用信号函数、简谐函数的预测结果为Ⅴ度，较采用复合函数的预测结果Ⅵ度相差不超过1度。图 8表明冻震可能产生的最重烈度不会超过MMI的Ⅵ度，这与表 2“冻震的形成条件与震颤特征”中3c项^[1]对冻震烈度的归纳总结结论一致。

图 9对比了冻震地面加速度模拟时程曲线与典型构造型破坏性大地震1952年的Taft-339地面加速度实际记录时程曲线，对比结果表明冻震可能引起的最大PGA不会超过该典型构造型破坏性大地震的PGA的10%。因此，冻震作用可能会造成城镇建筑开裂，但不会导致更严重的、或破坏性的结构损伤。但是，当冻震由地震诱发时，冻震作用叠加于地震作用将会在一定程度上增大地震的破坏性。

4 结论

冻震是一种新型的自然灾害。本文针对冻震可能引起的最大PGA、冻震可能达到的最大烈度IMM进行定量估测研究，提出了基于能量假定的冻震模型。该模型能较好地解释冻震的各类特征，并预测冻震可能达到的最大烈度MMI不会超过Ⅵ度，与已有的冻震相关文献记载及新闻报道的结果一致，即自1819年至今尚未发现烈度MMI高于Ⅵ度的冻震现象，从理论研究的角度给出冻震最大烈度水平，为专业人员进行相关结构研究提供设计参考。不同的震颤函数选择可能会降低预测结果的具体数值，但差别不大，如采用信号函数、简谐函数预测的最大烈度较采用复合函数预测的结果降低不超过1度，总体较为一致。同时，本文补充了一条冻震震颤特征，即“产生多处地表结构断裂，各处断裂分时发生，总持续时间相对较长”，并应用所提出的模型较好地解释了包括该条在内的冻震的地质条件、气象条件以及震颤特征。总体而言，单独的冻震作用不会造成城镇建筑的明显破损，但有可能叠加于地震作用上从而增大后者的破坏性。