祝杰 李婉秋 陶志刚 刘洋洋 马未宇
1)中国地震台网中心,北京 100045
2)山东建筑大学测绘地理信息学院,济南 250101
地震的孕育和发生是地应力长期累积、失稳,并快速释放的结果,目前地应力无法通过直接观测快速监测。Gorny等(1988)发现加兹利地震前出现大范围红外长波辐射异常增强现象,开拓了利用卫星热红外遥感技术研究地震的新途径。刘德富等(1997)首次利用地面射出长波辐射(OLR)资料分析我国1976—1985年龙陵、唐山、松潘、乌恰、共和5个7级以上强震震前一个月的OLR月均距平辐射结果,发现在强震区月均距平辐射显著增强并达到历年同期最高值。为探索震前热红外异常机理,各国学者通过分析多物理参量、多卫星数据资料源,探测地表面热异常分布与地震的关系(马瑾等,2020;Ouzounov et al,2007; 康春丽等,2007; 马未宇等,2020),相关的震例研究显示震前热红外异常与地震存在关联性。
以往研究缺乏震前红外预测与地震形变监测对比的实证性观测结果。震前红外长波辐射异常是电磁辐射能量的变化,其热力来源是下垫面的热量交换,包含了构造和非构造因素(如地形地貌、气象、季节等环境因素),从构造和非构造杂糅的热红外信号中找到与地震活动真正相关的信息,获得震前预测的红外辐射异常区域是否与地震构造形变区域对应是关键性问题(马瑾等,2005)。InSAR技术因其高精度、高时空分辨率、全天时全天候广域监测等特性,在地震同震、震间和震后形变监测、地面沉降监测、地质灾害隐患识别等领域取得诸多成功应用(王琪等,2020; 王超等,2002; 张景发等,2008; 许强等,2019; 徐小波等,2020)。联合InSAR技术与OLR辐射两种空间对地观测技术是研究地表形变多参量异常空间耦合关系的新思路(申旭辉等,2018)。然而,InSAR提取的震前地表微弱形变信号可能会湮没在大气延迟相位或地形误差相位中(Zebker et al,1997; 单新建等,2009),试图通过震前的形变参量异常分析震前红外异常与地震构造应力变化的关系,往往存在较大困难。InSAR技术为研究物质迁移的地壳应力应变场响应提供了新手段,通过发震断层破裂后的InSAR同震形变反推发震断层发生的构造应力积累及变化过程,学者在此方面开展了一系列研究(许才军等,2014; 单新建等,2017; 余鹏飞等,2022)。因此,本文尝试通过分析门源MS6.9 地震震前红外长波辐射演化过程和InSAR同震形变场,探讨震前红外预测异常和同震地表形变(构造应力)的对应关系。
北京时间2022年1月8日1时45分,青海省海北州门源县发生MS6.9 地震,震中位于37.77°N、101.26°E,震源深度10km,距离震中最近的乡镇为门源县皇城蒙古族乡和峨堡镇。地震造成地表破裂约22km,极震区最大烈度达Ⅸ度。此次地震是该地区近年来发生的最大地震,宁夏银川市、中卫市、甘肃张掖市、武威市、金昌市、兰州市等多个城市震感强烈。震后,国内外多家机构发布了门源地震的震源机制,不同机构测定震源机制使用的数据不同,导致震源机制存在一定差异(表1),但震源机制解均表明门源地震是一次走滑为主的地震事件。
表1 不同机构提供的门源地震震源机制
此次地震震中发生在祁连山内部冷龙岭断裂带北侧(图1),大地构造上位于青藏高原东北缘与阿拉善地块南缘之间隆起区的褶皱带内(郭鹏等,2017; 韩帅等,2022),区域内部由西向东发育有肃南—祁连断裂、托莱山断裂、冷龙岭断裂、金强河断裂、古浪断裂、毛毛山断裂、老虎山断裂、海原断裂等,组成以左旋走滑为主要特征,总长超1000km的祁连—海原断裂带在青藏高原板块东北缘应力分配和地壳形变中起着极其重要的转换调节作用(徐锡伟等,2007; 孙赫等,2017; 李智敏等,2022),此褶皱带在NE向构造应力作用下,发生了NE向的挤压缩短、顺时针旋转和SEE向挤出等构造变形(袁道阳等,2004),构造环境极为复杂。发震断裂冷龙岭断裂位于祁连—海原断裂带中西段,是祁连—海原断裂带“天祝地震空区”(Gaudemer et al,1995)的重要组成部分。冷龙岭断裂东起甘肃省天祝县双龙煤矿以东,西至门源—祁连公路的八道班,全长约120km,总体走向 60°N~70°W,倾向NE,倾角50°(何文贵等,2010),由一组近乎平行的NWW向分布的断裂组成,全新世以来活动强烈,主要表现为左旋走滑运动并伴有正倾滑性质(李振洪等,2022)。
射出长波辐射(Outgoing Long-wave Radiation,OLR)是指“地-气”系统通过红外窗口通道(4~120μm)向外部空间发射的电磁波能量密度,单位为W/m2。OLR对海洋表层和近地表温度的变化响应敏感,常用于监测与“热”成因现象相关的地学灾害征兆(马未宇等,2014)。本文利用美国NOAA极轨气象卫星提供的OLR产品,分析门源MS6.9 地震震前红外辐射演化。NOAA卫星自1974年发射至今已累计存档近50年的长波辐射OLR观测资料,提供白天和夜间两组数据,时间分辨率高,适合开展逐日的变化趋势研究,在全球气象分析业务和研究领域取得了较好的应用(任静等,2015; 孟庆岩等,2014)。
采用基于数字图像增强技术的周均涡度距平法处理2021年12月22—28日覆盖甘肃青海地区的OLR产品(空间分辨率1°×1°,时间分辨率1天),提取的红外长波辐射显著增强区结果如图2所示,并依据此图,预测未来地震危险区域。由图2可见,青海门源地区出现明显的OLR辐射异常增强区,增强沿青海德令哈—西宁—甘肃武威一带分布,呈“哑铃”状近WE向水平展布,“哑铃”两端有2个辐射增强涡度中心。辐射增强条带西起青海德令哈,在德令哈周边凝聚一个涡度值中心,最大涡度值20W/m2,自西向东OLR涡度值逐渐减小至12W/m2,辐射增强条带东侧在甘肃武威周边凝聚了一个涡度值16W/m2的弱中心。门源MS6.9 地震震中位于德令哈涡度中心至武威涡度中心过渡的弱涡度增强区,辐射增强维持在16W/m2。
图2 2021年12月22—28日甘肃青海地区OLR辐射涡度距平空间分布
针对震前预测的异常区域,进一步开展OLR异常日变跟踪(2022年1月1—11日),如图3所示。由图可见在震前1月2日,震中东南部红外异常开始出现,1月4—6日震中东南部异常持续并进一步加剧,1月7日,震前1天异常现象消退并消失,至1月8日震中东南部异常再次出现后主震发生。震后1月9—10日,异常逐渐衰退,1月11日红外异常全部消失。OLR异常空间演化经历了起始增温—异常加强—高峰—衰减—平静的演化过程,与岩石在受力破裂过程中的空间演化(挤压—岩石微破裂—岩石破裂扩展—应力闭锁—破裂终止)具有基本一致的特征(吴立新等,2004),基本反映出岩石在构造应力加载—准静态成核—动态破裂—应力重分布—断层强度再恢复的时间演化过程。
图3 门源 MS6.9 地震OLR辐射涡度距平时空演化
此次门源MS6.9 地震OLR异常变化过程与2013年4月20日芦山MS7.0、2021年5月22日玛多MS7.4、2020年7月12日唐山MS5.1、2020年7月23日西藏尼玛MS6.6 地震过程中热辐射变化规律基本一致(马未宇等,2014、2020; 王淑艳等,2021; 姜祥华等,2021)。这5次地震震前一天均出现OLR异常衰退消失的现象,与吴立新等(2018)在应力加载红外观测实验中发现的岩石应力闭锁阶段相对应,可能预示着震源区岩石介质地应力闭锁的出现,为地震的最终到来提供指示。
以上分析了OLR异常的空间演化与岩石破裂过程演化的对应性。为研究震前预测区OLR增强是否与地震形变区域相关,本文进一步利用InSAR同震形变场分析门源MS6.9 地震前的预测区红外异常是否由地震构造应力变化引起。
选取欧洲航天局的Sentinel-1A升、降轨影像进行雷达差分干涉处理。Sentinel-1A数据为C波段,图像分辨率5m×20m,VV极化,幅宽250km。升轨干涉对时间基线12天,空间基线38.24m; 降轨干涉对时间基线12天,空间基线55.56m(表2)。较小的时间基线和空间基线保证了干涉相对较高的相干性。
表2 Sentinal-1A干涉对信息
基于二轨差分干涉法,采用瑞士GAMMA软件生成差分干涉图。对于Sentinel-1A SAR数据,利用欧洲航天局提供的POD精密定轨星历数据进行SLC数据轨道校正,去除轨道误差,利用30m分辨率的SRTM DEM数据去除地形相位并进行地理编码。为使Sentinel-1A干涉像对保持较高的空间分辨率并达到去除斑点噪声的目的,采用10︰2的距离向、方位向视数比。采用加权功率谱法对干涉图进行3次滤波,滤波窗口依次设置为128、64、32个像素,以提升干涉图的相干性。采用最小费用流法(Eineder et al,1998)进行相位解缠,相干性阈值设置为0.3,对于相干性阈值以下的区域不进行解缠计算。对于2个差分干涉图中的大气垂直分层延迟相位,基于各自的数字高程模型建立延迟模型,并从原始差分干涉相位中去除,获得去除各种相位误差、清晰的干涉条纹图(图4(a)、4(b))。最后,将解缠后的干涉相位转化为地表沿雷达视线方向形变量,获得门源MS6.9 地震升、降轨同震形变场(图4(c)、4(d))。
图4 门源 MS6.9 地震升、降轨同震形变分布
图4(c)、4(d)分别为Sentinel-1A升、降轨数据提取的门源MS6.9 地震雷达视线向同震形变,形变信息连续完整,升轨和降轨同震形变的空间分布基本一致。同震形变主要位于肃南—祁连断裂(俄堡段)、托莱山断裂和冷龙岭断裂的交汇区,形态呈“蝴蝶状”,长轴方向为NWW向,长度约40km,短轴方向为NNE向,长度约30km,即形变区的范围约40km×30km。形变场被高分7号影像解译出的NWW-SEE走向的断层地表破裂迹线从中间分割,地表破裂迹线与托莱山断裂—冷龙岭断裂基本重合,与多家机构测定的震源机制解走向参数(表1)基本一致。图4中形变正值表示靠近卫星飞行方向,即LOS向缩短,形变负值表示远离卫星飞行方向,即LOS向拉伸。T128升轨和T33降轨LOS向缩短的最大形变量分别为0.48m和0.62m;LOS向拉伸的最大形变量分别为0.61m和0.53m。升、降轨极震区的同震形变均存在较严重的失相干现象,主要是因为发震断层冷龙岭断裂破裂到达地表,导致地表发生了大梯度形变。此次地震造成的最大形变量约0.62m,位于冷龙岭断裂的南盘。升轨和降轨数据提取的同震形变在冷龙岭断裂上、下盘表现为运动方向相反,形变量相当,同时同一轨道的形变在上、下盘也表现为相反的运动特征,表明门源MS6.9 地震引起的地表形变主要以水平方向为主,断层的运动性质具有典型的走滑变形特征(李振洪等,2022;Zhang et al,2022)。
进一步分析震前预测的OLR辐射异常与InSAR同震形变场的关系可以发现,2021年12月22—28日震中周边出现明显的热辐射异常增强区,沿青海德令哈—西宁—甘肃武威一带呈“哑铃”状近WE向水平展布,OLR呈现顶峰值—弱增强—次峰值的空间分布,震中位于青海德令哈涡度中心至甘肃武威涡度中心过渡的弱涡度增强区。其可能原因是地震较易发生在与构造活动有关联的OLR变化高梯度带上,芦山MS7.0 和唐山MS5.1 地震的OLR分布与震中也有类似的位置关系(马未宇等,2014、2020)。对照重力学科,重力变化高梯度带也是重力异常判断的一项重要指标(祝意青等,2022)。2022年1月8日门源MS6.9 地震发震后,利用InSAR技术提取的极震区同震破裂位于肃南—祁连断裂(俄堡段)、托莱山断裂和冷龙岭断裂的交汇区,形变区范围40km×30km,同震形变场被托莱山断裂和冷龙岭断裂从中间分割,南、北盘呈现相反的运动特征。升、降轨同震形变集中分布于“哑铃”状红外辐射异常区的内部,具体位于德令哈涡度中心至武威涡度中心过渡的弱涡度增强区,东侧临近武威涡度中心(图5)。由于OLR长波辐射是通过卫星平台被动获取的下垫面射出长波辐射,辐射能在传输过程受到大气辐射增温进一步发散,且OLR图像空间分辨率(1°×1°)远远低于SAR图像空间分辨率(30m×30m),二者在一定程度上导致OLR异常在空间分布上比InSAR同震形变范围大,但震前预测的红外辐射增强区与InSAR同震破裂形变区的空间位置基本吻合、扩展形式基本相似(同震破裂形变区分布在红外辐射异常区内部),间接反映了门源地震前热异常与地震构造地应力释放具有显著关联性,震前的红外辐射预测异常是此次门源地震前构造应力状态变化的遥感影像反映。
图5 震前OLR显著增强区与同震形变场的对应关系
通过分析震前预测的青海门源MS6.9 地震的射出长波辐射短期异常时空演化特征,进而利用InSAR技术提取门源地震升、降轨同震形变场,探讨了震前红外预测异常与同震地表形变(构造应力)的关联性,得到主要结论如下:
(1)震前OLR的空间演化特征显示,在极震区周边出现了明显的热辐射增强异常,青海德令哈—西宁—甘肃武威一带出现了呈“哑铃”状近WE向水平展布的OLR热辐射增强区,空间演化从顶峰值—弱增强—次峰值,具有唯一性和高空间可辨识性的特征,证明了强震前震中附近存在热红外辐射增强现象。虽然暂无法对此次地震前的红外异常机理做出明确解释,但通过卫星红外辐射场的动态化数值分析研究强震前异常指征,对地震预报具有一定程度的指示意义。
(2)利用Sentinle-1A升、降轨数据提取的同震形变主要位于肃南—祁连断裂(俄堡段)、托莱山断裂和冷龙岭断裂的交汇区,形态呈“蝴蝶状”,发震断裂冷龙岭南盘最大形变量约0.62m,北盘最大形变量约0.5m。南、北盘运动方向相反,量级相当,揭示门源MS6.9 地震引起的地表形变主要以水平方向为主,断层运动具有典型的走滑变形特征,与震源机制解相符。
(3)门源MS6.9 地震OLR短期异常时空演化过程遵循了岩石应力加载破裂过程中的热异常规律,显示热异常变化与应力变化存在关联。多个震例研究显示的震前一天出现OLR异常衰退消失现象可能表征岩石介质应力闭锁的出现,可为地震的最终到来提供指示。震前热红外预测异常区域与InSAR技术提取的同震破裂形变区域高度对应,反映了门源震前热异常与地震构造地应力释放具有显著关联性,震前的红外辐射预测异常是门源MS6.9 地震前构造应力状态变化的遥感物理参量的反映。此外,本研究综合应用InSAR和热红外遥感两种空间对地观测技术,分析门源地震前热异常与地震构造应力变化之间的关联,为后续建立空间、陆地、海洋相结合的综合协同观测体系提供了探索性尝试。
需要指出的是,此次门源MS6.9 地震的红外异常与地震构造应力变化存在关联是通过InSAR同震形变场间接说明的,未来仍需要深入研究InSAR震前形变参量,解算出能够真实反映震前地表构造运动的形变信息,进而与震前红外参量进行同步立体式耦合,建立地球物理多参量异常监测模型,提高特定区域强震前的早期预测能力。
致谢:审稿专家对本文提出宝贵的修改意见,欧洲航天局提供免费的Sentinel-1A数据和精密轨道数据,NOAA提供了OLR数据,在此一并表示感谢。