2018年5月28日1时50分52秒，吉林省松原市宁江区（45.27°N，124.71°E）发生M_S 5.7地震，震源深度10 km。经调查，吉林省大部分地区、黑龙江省部分地区有震感，其中，吉林松原市区、前郭县、大安市震感强烈。国内外许多研究者一直致力于地震前兆的探究，随着卫星遥感技术的日渐成熟，震前热红外异常也成为了的热点研究领域（强祖基等，1998；单新建等，2005；屈春燕等，2006）。地表温度场中存在着与活动构造带密切相关的热信息，地震前后温度场的变化与同震形变的异常间的关系得到了进一步的印证（马瑾等，2010；陈顺云等，2014）。随着该领域研究的不断深入，又出现了一些新的研究方法。张元生等(2010, 2011)采用基于离散小波变化滤波和短时傅里叶变换的相对功率谱方法研究2008年汶川M_S 8.0、2011年日本东北M_S 9.0地震前的热红外量温异常现象；郭晓等（2010）在研究强震长波辐射异常变化时提出功率谱相对变化法、历年同期亮温偏移指数法，用于研究2008年汶川M_S 8.0、2010年玉树M_S 7.1地震，并得到2次地震的亮温异常时空演化规律（温少妍，2011）；此外，解滔等（2015）使用“连续小波变换法”研究了2014年于田M_S 7.3地震的亮温异常。张元生等（2010）根据影响地表热辐射的因素所具有的不同周期，利用小波变换法去除影响因素，从而对地震热异常进行有效提取，并对2008年汶川M_S 8.0、2014年于田M_S 7.3、2008年乌恰M_S 6.8地震进行了热辐射变化分析，结果显示在3次地震发震前均存在明显的热异常，且发震区域一般位于异常区边缘及区内的活动断层上，发震时间在功率谱达到峰值之后的35天内。其它的相关研究还有很多（张璇等，2013；李青梅等，2015），这些研究在发震地点、时间上意义明确，可为短期地震预测提供一定的经验。

本文搜集了震中周围（38°—50°N，115°—135°E）2013年1月1日至2018年6月18日静止风云二号FY-2G气象卫星观测得到的亮温资料，采用小波变换和功率谱估计法作为数据处理方法，研究2018年5月28日松原M_S 5.7地震前热辐射数据，以期为东北地区中强震前热红外异常研究提供震例参考。

2018年松原M_S 5.7地震是继1975年海城7.3级地震之后东北地区震级最大的一次浅源地震。据中国地震台网测定，震中位于45.27°N、124.71°E，震源深度为10 km（http://ceic.ac.cn/）。震中位置和断层构造图如图 1所示。主震后余震次数超过823次，其中，ML≥3.0余震8次。此次地震发生在松辽断陷盆地中央沉降带内，震中位于NE向松原（扶余）—肇东断裂带和NW向第二松花江断裂带交汇部位。地震活动空间分布图像（邵博等，2016）显示，松原地区的地震活动主要沿NE向松原（扶余）—肇东断裂带分布，且具有丛集分布特征，主要分布在前郭县查干花和松原市宁江区2个区域，地震活动明显受NE向断裂控制。根据松原市活断层探测结果（万永魁等，2016），NE向松原（扶余）—肇东断裂带在松原市一带表现为2条近平行的NE向或近EW向的断层段，即扶余北断层和江北断层，前者为晚第四纪活动，后者为前第四纪活动。扶余北断层长约23 km，表现为上陡下缓的南倾特征，局部为逆断层性质。松原M_S 5.7地震等烈度线长轴为NEE向，表明此次地震受控于NE向断裂活动或至少NE向断裂参与了活动。

图 1(Fig.1)

图 1 松原MS 5.7地震震中和周围区域断层 Fig.1 The epicenter of Songyuan earthquake and the main active faults in the study area

普朗克辐射定理完整地描述了黑体辐射能量密度随频率的分布规律，其公式为

式中，B(υ, T)为黑体光谱辐射度（W·m^-2·sr·cm^-1）；υ为波数（cm^-1）；T为温度；k为玻尔兹曼常数；c为光速；h为普朗克常数。对于任意波数从υ_i1到υ_i2的波段，黑体总辐射度W_i(T)为普朗克公式对这一波段各个波长值和该波段响应函数的积分，可以表述为

式中，i为波段序号；f_i(υ)为该波段的响应函数；υ_i1、υ_i2为该波段波长上、下限对应的波数。将波段以一定的波数间隔Δυ进行离散，式（2）积分可以用以下求和表达式来近似。

其中，m_i＝(υ_i2 - υ_i1)/Δυ表示波段的离散化个数。利用式（3），在温度T_i1至T_i2范围内，以ΔT为步长，分别对每个红外波段进行积分，可以得到每个波段总辐射度W_i(T)与温度T的对应关系；再将卫星观测资料辐射定标和几何校正，得到地表像元的辐射度，然后查询辐射度与温度的对应表，即可得到对应地表的辐射温度值。

FY-2G静止气象卫星于2008年发射，运行轨道定点于105°E赤道上空，有效观测范围为60°S —60°N、45°—165°E。红外观测波段分别为10.3—11.3 μm、11.5—12.5 μm，热红外量温每小时观测1次，空间分辨率约为5 km。静止卫星相比极轨卫星温度图像及小波相对能谱计算无需轨道拼接过程。日夜温差过大是影响区分提取地震信息的主要原因，因此选取地方时00:00—04:00时间段的5次观测数据进行分析。

小波变换是一种分析非稳态信号的有效方法，在时间域和频率域中得到了很好的演化，已经被应用于地球物理学和地震勘探等领域（Kumar et al，1997）。有限时间序列的小波变换定义为

“*”表示对$\psi _{a, b}^*$取复共轭；${\psi _{a, b}}(t) = \frac{1}{{\sqrt a }}{\psi _{a, b}}\left({\frac{{t - b}}{a}} \right)$为小波母函数；a为尺度伸缩因子，用于控制变换的频率分辨率，a的离散间距越小，频率分辨率越高；b为时间平移因子（Torrence et al，1998）。这里采用Morlet小波进行分析，常用的复值Morlet小波在频率域被定义为

在ω₀≥5时，式（5）近似满足容许条件，此处ω₀将值取为6（Farge，1992）。由于Morlet小波在时间域和频率域均为复数，经小波变换后的各频段信号分量W_ψf (a, b)也是复数，因此可以提取信号的振幅和相位信息。小波能谱通常定义为振幅的平方|W_ψf (a, b)|²，由于各像元所处的纬度、海拔高度以及气候环境等因素的不同，使得不同像元的能谱信息本身存在差异，从而在空间上造成虚假的异常现象，但这并不是我们想要得到的结果。因此，对每个像元评估功率谱与其平均值的变化，这被称为相对小波能谱（RWPS）。例如，对于某一频段，6频段相对小波能谱意味着在整个分析阶段取其平均值的6倍，各频段的相对能谱变化为

其中，${{\bar W}^2}(a, b) = \frac{1}{N}\sum\limits_{l = 0}^{N - 1} {{{\left| {{W_l}f(a, b)} \right|}^2}} $，为各频段全局小波能谱，表示在整个时段内将小波能谱取平均；N为数据长度。

Torrence等（1998）讨论了连续小波变换的边界效应，结果表明在临近数据边界处小波能谱呈衰减现象，且频段越低，衰减间隔时间越长。这表明，在边界处计算的能谱值要低于数据实际的能谱值。如果在边界处出现能谱异常，那么实际的异常幅度会更大。因此，小波能谱分析方法的边界效应不会产生虚假的异常。

尺度伸缩因子a和数据长度N控制频率分辨率和整个频段数量。T为对应频段周期，计算公式为

式中，a为尺度伸缩因子；Δt为时间序列采样间隔；i为频段序号；δ_i为对a进行离散化的间隔。值得注意的是，许多频段是由在小波变换中的尺度伸缩因子和数据长度确定的，但我们仅计算了周期为8—64天的能谱（张元生等，2010）。连续小波变换是非正交变换，相邻频段信息中有重叠的成分，因此不需要太高的频率分辨率，这里取δ_i＝0.5为间隔，对a进行离散化。选取2016年1月—2018年1月的数据用于相对功率谱的分析，通过对比来识别在正常背景值下的异常变化。依据上面给定的参数，共有6个频段，对应周期依次为8.26天、11.70天、16.53天、23.38天、33.06天、46.75天。

首先，用补窗法对各像元每天5个时间点的亮温数据进行处理，以去除短时间范围内云层对数据的部分影响，再取其平均值作为每天该像元的观测值，对所有像元每天的数据作同样的处理得到各像元的日均值序列。然后，利用小波变换对每个像元进行2阶尺度函数与7阶尺度函数相减的处理，以去除连续亮温数据信息中所包含的地球基本温度场（直流部分）、年变温度场、日变温度场、雨云和寒热气流等非震因素引起的温度变化，保留由地震因素引起的温度微变化信息，得到在时间域里正负相间的亮温相对变化波形数据。考虑到短临地震异常出现的时间一般为震前10—90天，故以n=64天为窗长、m=1天为滑动窗长对小波变换后的每个像元作傅里叶变换，计算其功率谱，并对每一像元的所有频率的功率谱作相对幅值处理，得到功率谱相对变化的时空数据（张元生等，2010；郭晓等，2010）。最后，利用计算得到的时、频空间数据进行全时空和全频段扫描，提取并识别异常区域时空变化特征。

对6个频段的小波相对能谱时空演化图（图 2）分析后发现，在松原M_S 5.7地震发生前，第5频段呈现出显著的热红外异常现象。对地震前1个月的第5频段作小波能谱时空演化分析（图 3），由图 3可见，2018年5月25日，首先在吉林—丰满断裂带周围出现亮温异常，初期异常幅度及区域面积都较小，随后异常区域沿着郯庐断裂带延伸扩展。5月27日异常区域面积进一步扩大，5月28日距异常区域约200 km处发生2018年5月25日了松原M_S 5.7震。地震发生之后异常区域面积进一步延伸扩展，6月12日异常面积达到峰值，而后逐渐缩小，6月17日异常基本消失。

图 2(Fig.2)

图 2 东北地区2018年5月28日第1—6频段小波相对功率谱 Fig.2 Relative power spectrum value of 1-6 band wavelet in Northeast China on May 28, 2018

图 3(Fig.3)

图 3 松原MS 5.7地震前后第5频段小波相对功率谱的异常时空演化 Fig.3 Time-space evolution of the 5 band RWPS anomalies before Songyuan earthquake

为了分析松原M_S 5.7地震前热红外亮温功率谱的变化过程，对图 2中功率谱高值区域的平均值进行研究，其相对幅值时间变化序列如图 4所示。由图 4可见，2013—2016年研究区域亮温平均功率谱相对幅值基本在6以下，2016—2018年亮温平均功率谱相对幅值有4次较大的高幅度异常。其中，2016年3—5月研究区平均功率谱相对幅值高值异常，但是研究区该时间段的地震活动水平较低；2017年6月16日—7月27日研究区平均功率谱相对幅值高值异常，其间发生了2017年7月23日松原4.9级地震。2018年5月26日相对幅值为4.04，5月27日相对幅值达到6.13，发生地震的5月28日相对幅值为8.22，地震之后相对幅值持续增高，6月12日达到峰值，为12.36，相对幅值峰值之后持续下降，直至6月27日下降到6.01。

图 4(Fig.4)

图 4 松原MS 5.7地震前研究区域亮温平均功率谱相对幅值 Fig.4 Time series curve of mean RPSW in Songyuan seismic anomaly area

卫星观测得到的热红外亮温资料反映的是像元地表温度，小波相对能谱的计算是对每一像元单独进行的，因而每一像元的计算结果不会受到其他像元的影响。图 2中异常区域的扩大过程显示出越来越多的地区在松原M_S 5.7地震前地表温度出现了变化。在以往震例研究中发现，低频段的相对小波能谱异常要少于高频段的，并且低频段的地震热红外异常可信度也高于高频段的。然而，考虑到地震发生前出现的异常，因此仍有必要对地震前热红外异常进行分析研究。

由于岩石的导热系数很小，因此在稳定的地质构造区域由构造活动引起的额外热量很难扩散到地表。但是在活跃的构造区域，地层中的地下流体等可以通过活跃构造单元快速地扩散热量（Tronin，2000）。松原M_S 5.7地震前吉林—丰满断裂带周围开始出现异常，随着地震的临近，异常区范围加速扩大，异常幅度也快速增加，这种现象表明震源区应力已接近岩层破裂的临界状态，贯通地表的裂隙数量增多，地壳内部的温室气体溢出，由于裂隙突然增多，地下岩层与地表的热对流亦有所增强。

本文应用小波相对能谱方法对松原M_S 5.7地震震中周围（38°—50°N，116°—134°E）2013年1月1日—2018年6月18日的中国静止气象卫星FY-2G观测的亮温资料进行了分析研究，使用连续小波变换方法获得了松原M_S 5.7地震前的相对小波能谱演化。我们发现在2018年5月25日，吉林—丰满断裂带周围首次出现了明显异常。随着地震的临近，异常区域向震中区延伸。地震发生前整个异常过程持续时间仅为3天，震后异常区域面积及异常幅度持续增大，直至震后1个月才逐渐恢复正常。油气盆地地下大量天然气对震前应力的变化较为敏感，当应力积蓄到一定程度，盆地周缘的活动构造带及一些微裂隙都是天然气上涌的通道，溢出地表的甲烷、二氧化碳等温室气体辐射增温效果明显。吉林—丰满断裂带的发育程度较好，油气等通过发育较好的通道快速到达地表，这或许可以解释松原M_S 5.7地震的震中不在异常区内的原因。