利用InSAR技术研究新疆阿什库勒火山群现今活动性

季灵运 许建东 赵 波 万 园 申欢欢

1)中国地震局地质研究所，活动构造与火山重点实验室，北京 100029

2)中国地震局第二监测中心，西安 710054

3)中国地震局地质研究所，北京 100029

0 引言

阿什库勒火山群位于新疆于田县以南250km处的阿什库勒盆地内部(图1)。阿什库勒盆地位于康西瓦断裂和阿尔金断裂的交会处，构造运动活跃，中强地震频繁发生(图1)。阿什库勒火山群由10余座主火山和数十个子火山组成，区域海拔超过4 600m，是世界上海拔最高的火山区之一，熔岩面积达200km2(刘若新，2000)。新疆日报曾报道阿什库勒火山群于1951年5月27日爆发过(新疆日报，1951)，但存在争议(邓万明，1989;刘嘉麒等，1990;刘若新，2000)。2011年5月中国地震局地质研究所和新疆维吾尔自治区地震局对阿什库勒火山群开展了科学考察，在阿什火山北侧火口沿发现了疑似1951年喷发遗留的喷气孔(许建东等，2011)。如果1951年的喷发事件真的存在，那么阿什库勒火山群现今是否还在活动?是否具有潜在的喷发危险性?另外，2008年3月21日，阿什库勒火山群南部山区发生了MS7.3强烈地震(陈学忠等，2008;图1)，本文称之为2008年于田地震，此次强震事件对火山活动是否有影响?

图1 阿什库勒火山群区域构造与地震震中分布图Fig.1 Tectonic and epicentral distribution map around Ashikule volcanic group.

阿什库勒火山群地区自然条件极其恶劣，研究程度较低，目前的研究主要集中在火山地质与年代学方面(Liu et al.，1989;刘嘉麒等，1990;许建东等，2011)，迄今为止，还没有地表形变方面的报道。恶劣的气候环境使得人类难以接近，因而以地面控制点为基础的常规地表形变监测手段如GPS、水准测量很难实施。近年来，作为一种全新的空间对地观测手段，InSAR(Interferometric Synthetic Aperture Radar)技术因其具有不需要地面控制点，且能够以高精度、低成本、快速度的优势完成大范围、高空间分辨率的“面”的地表形变信息，而广泛应用于地学研究的各个领域(Massonnet et al.，1993，1995;单新建等，2006;孙建宝等，2008)。

本文基于Envisat ASAR和ALOS PALSAR卫星影像资料，利用干涉图堆叠技术(Stacking)和短基线集合成孔径雷达干涉测量技术(SBAS-InSAR，Small BAseline Subsets-Interferometric Synthetic Aperture Radar)提取阿什库勒火山群地区的地表形变，了解并分析其在2008年于田地震前后的活动特征，填补该火山群现今活动性未知的空白，进而为该火山群的喷发预测、灾害评估等工作提供基础资料。

1 SAR影像资料与InSAR数据处理

1.1 SAR影像简介

为了检测阿什库勒火山群地区在2008年3月于田地震前后可能存在的形变信息，我们购买和收集到了Envisat ASAR和ALOS PALSAR两种卫星影像资料。其中Envisat ASAR影像共41景(1个航迹，降轨)，其中于田地震之前20景，震后21景，时间跨越2003—2010年;ALOS PALSAR影像共8景(1个航迹，升轨)，其中震前3景，震后5景，时间跨越2007—2009年。

1.2 InSAR数据处理方法

首先，对于所有影像，不指定公共主影像，进行自由组合干涉。考虑到阿什库勒火山群地区地处高原无人区，地表无植被覆盖，长时间间隔仍能保持较高的相干性，但山区地形起伏较大，宜采用较短垂直基线的干涉图，所以给定干涉限制条件为:时间间隔短于5a，垂直基线＜300m(Envisat ASAR)与＜3 000m(ALOS PALSAR)。

对于每个符合条件的干涉对，进行常规两轨法D－InSAR处理，使用美国NASA公布的SRTM DEM(分辨率约为90m)去除地形相位部分。目标是提取火山区可能存在的大范围地表形变异常，所以干涉图的多视因子选定4×20(距离向×方位向)。则对于Envisat ASAR影像生成的干涉对，每个像元约为80×80m2，对于ALOS PALSAR，每个像元约为60×60m2。干涉图相位解缠采用最小费用流算法(Minimum Cost Flow)，仅解缠相干系数值＞0.5的像元。对含有基线误差的干涉对，基于高程信息和解缠相位，采用非线性最小二乘平差方法去除基线误差(Rosen et al.，1996;Lu，2007)。对明显含有与地形相关的大气相位延迟(Atmospheric Phase Screen，APS)的干涉对，根据无形变区域的解缠相位和高程信息，建立大气相位延迟模型，去除了与地形相关的大气相位延迟干扰。

削弱各种噪声，提取缓慢变形信号的一种常用的方法是干涉图堆叠技术(Zebker et al.，1994，1997;Peltzer et al.，2001)。干涉图堆叠技术通过对覆盖同一地区的多个干涉图进行加权平均(一般以每幅干涉图的时间间隔为权重)，实现削弱每个干涉对所含的大气延迟相位以及轨道残余误差，达到提高最终形变结果可靠性的目的。

为了揭示阿什库勒火山群地区在于田地震后地表形变随时间的演化特征，本文采用了短基线集算法求解地表形变时间序列。短基线集算法是由Berardino等(2002)提出来的，后来由Jung等(2008)进行了精化改进。短基线集算法以分布式散射体(DS，Distribute Scatterers)为研究对象，允许以相对短的空间基线和时间间隔对多景影像进行自由组合干涉，增加了时间采样率并提高了干涉图的相干性，削弱了空间失相干的影响，降低了外部DEM引入的误差，从而达到提高形变结果精度的目的。由于大气延迟相位在空间域上高度相关，而在时间域则是随机的(Ferretti et al.，2000;Berardino et al.，2002;Jung et al.，2008)，短基线集算法通过时间域的高通滤波和空间域的低通滤波来去除大气相位延迟。精化的短基线集算法在以下4个方面提高了时间序列形变结果的可靠性:1)选取高质量的干涉图(不含相位解缠误差，且相干性好)估计外部DEM误差和线性形变相位;2)对形变时间序列和误差进行迭代计算;3)利用有限差分平滑方法削弱时间域的噪声;4)对参考点本身存在的误差进行了校正(Jung et al.，2008)。

2 InSAR形变结果分析

2.1 2008年于田地震前阿什库勒火山群地区地表形变特征

通过逐一分析2008年于田地震前的所有干涉图，没有发现明显的形变信息，选取18个相干性较好的干涉图，采用干涉图堆叠技术，计算得到了阿什库勒火山群地区2003—2007年的平均形变速率图(图2)，其中黑色三角形表示各个小火山口的位置。可以看出，阿什库勒火山群地区相干性较好，南部山区由于终年积雪导致失相干现象较严重，但沟谷地带仍能保持较好的相干性。形变为沿雷达视线方向(LOS，line of sight)，可以看出在2008年3月于田地震发生之前，阿什库勒火山群地区没有显示出明显的地表形变，表明该区在2003—2007年期间相对比较稳定。

图2 阿什库勒火山群2003—2007年平均形变速率图(LOS)Fig.2 Average deformation rate of Ashikule volcanic group during 2003 －2007(Line of Sight).黑色三角形与图1相同;红色圆点表示参考点

2.2 2008年于田地震后阿什库勒火山群地区地表形变特征

通过逐一分析2008年于田地震后的所有干涉对，发现阿什库勒火山群和于田地震发震地区均存在不同程度的形变信息。通过进一步对比，发现形变与时间呈明显的线性关系。我们首先选取了Envisat 43个、ALOS PALSAR 6个相干性较好的干涉图，采用干涉图堆叠技术获得了震后的平均形变速率图(图3)。可以看出，阿什库勒火山群地区整体有沿雷达视线方向隆升的态势，但量级不大，＜1cm/a。南部山区没有明显的形变趋势，但是于田地震同震地表破裂带(徐锡伟等，2011)北端地区呈现了明显的沿雷达视线方向下沉的趋势，且量级较大，最大达到约4cm/a(雷达视线方向)。

考虑到大气相位延迟、外部DEM误差，选择短基线集算法求解形变时间序列，便于直观分析和更好地揭示研究区域震后形变随时间的演化特征。选取Envisat相干性较好的43幅干涉图，ALOS PALSAR的6幅干涉图，分别获取了研究区域基于2种卫星影像资料的形变时间序列特征。图4给出了所使用的干涉对的时间与空间垂直基线，实线代表Envisat ASAR影像干涉对，虚线代表ALOS PALSAR影像干涉对。可以看出，2种卫星影像资料各自生成了一个干涉图子集。而ALOS PALSAR的时间被包含在Envisat ASAR干涉对的时间内。图5列出了根据Envisat影像获取的阿什库勒火山群地区震后形变时间序列，即所有影像成像时间相对于参考影像(Envisat为20080420，ALOS PALSAR为20080711)的累积形变。从整体上看，阿什库勒火山群地区震后有沿雷达视线方向隆升的态势。值得注意的是，于田地震同震地表破裂带北端存在明显的震后形变，表现为明显的沿雷达视线方向下沉的形变趋势，形变与时间呈较强的线性关系。

图3 阿什库勒火山群于田地震震后平均形变速率图(LOS)Fig.3 Average deformation rate of Ashikule volcanic group after the Yutian earthquake(Line of Sight).

为了更加清晰直观地表现相应研究目标的形变时间序列演化，以阿什山火山以及于田地震同震地表破裂带北端为例，研究它们的形变演化特征。图6a，b分别列出了相应目标的形变时间序列，结合图3，5，可以得到如下结论:

(1)邻近参考点地区在时间上比较稳定，表明参考点选择合理。

(2)阿什山火山在Envisat(降轨，即卫星由N向S飞行)干涉结果上显示沿雷达视线方向隆升(即朝向卫星方向)，而在 ALOS PALSAR(升轨，即卫星由S向N飞行)干涉结果上则显示沿雷达视线方向下沉(即远离卫星方向)。然而，InSAR仅能观测到一维形变信息，即雷达视线方向形变是地表3方向(E、N、上)形变的共同贡献。根据雷达入射角以及卫星飞行方位角，可以计算得到雷达视线方向位移与地表3方向位移的几何关系，本文中，Envisat为［0.3719，－0.080，0.922］，ALOS PALSAR 为［－0.616，－0.110，0.780］。可以看出，2种卫星影像均对垂向位移最敏感，而对SN方向位移的敏感性最差。另外，注意到2种卫星影像对EW方向位移的表现是相反的，因而可以推断阿什山火山的形变以水平方向分量为主。通过对比图5，6可以得出，以阿什山火山为代表的阿什库勒火山群地区在震后有整体向E运动的趋势，量级较小，累积最大约1cm。这种形变分布不属于火山活动的形变特征，应为区域块体水平构造运动的反映。

图4 干涉对基线组合图Fig.4 Temporal and perpendicular baselines of InSAR pairs.

图5 阿什库勒火山群于田地震震后形变时间序列(LOS)Fig.5 Deformation time series in Ashikule volcanic group after Yutian earthquake.

(3)于田地震同震地表破裂带北端附近在2种卫星影像资料上均表现出了较强的下沉形变态势(沿雷达视线方向)，且与时间呈较强的相关性。考虑到图5中于田地震发震断层两盘在震后的位移，沿垂直于断层的走向做了一条剖面(位置见图3)，以详细分析其形变的空间特征。由升降轨道成像几何关系可知，SN方向形变在雷达视线方向上最不敏感，所以假定SN方向位移为0，联合2种卫星的成像几何关系(Wright et al.，2004)，求解得到于田地震同震地表破裂带北端的E向、垂向的形变速率分量，如图7所示。

图6 几个局部区域的形变时间序列(位置见图3a)Fig.6 Deformation time series at certain areas(see Fig.3a for locations).a Envisat;b PALSAR

由图7可以看出，于田地震地表破裂带北端的震后形变特征为水平方向整体向E运动，最大约5cm/a;垂直方向下沉，且距离断层越近下沉速率越大，最大约6cm/a。同时，形变速率分量与高程也无明显的相关关系。我们认为，此局部地表下沉不属于火山活动的影响范围:1)该区域距离火山群较远;2)如果确属岩浆运移引起的地表下沉，那么在未喷发的情况下，岩浆到达新的存储地点后，会引起局部地表隆升。因而，阿什库勒火山群在2008年于田地震后没有表现出明显的火山活动引起的地表形变异常。同震地表破裂带北端的形变异常应该是震后块体应力平衡调整的结果。

我们收集了阿什库勒火山群地区1970—2011年的地震目录(图1)，虽然区域内地震活动比较频繁，但没有明显的空间分布规律。火山群内部地震分布较少，与各个子火山也没有空间的相关关系。这表明，研究区域的地震活动可能与阿什库勒火山群无直接的联系，应该是区域内部发育的断层活动引发的构造地震。因此，在过去的40余年间，阿什库勒火山群地区没有发生与火山活动有关的地震事件。

图7 于田地震地表破裂带北端剖面(图3)的E向与垂向速度关系曲线Fig.7 Cross section of average deformation rates along the rupture caused by Yutian earthquake.

3 结论

(1)2008年于田地震之前，阿什库勒火山群地区总体形变不明显，表明阿什库勒火山群在2003—2007年期间相对较稳定。

(2)2008年于田地震之后，阿什库勒火山群地区整体有E向运动的趋势，到2010年7月，累积最大位移约1cm，为震后块体应力调整的结果，非火山活动引起。地表同震破裂带北端在震后显示出了较强的E向水平运动和垂向下沉运动的形变态势，累积最大约6cm，与火山活动没有必然的联系。

(3)1970—2011年的地震目录表明，阿什库勒火山群在过去的40多年以来没有发生火山地震活动事件。

陈学忠，蒋长胜，李燕娥.2008.2008年3月21日新疆于田7.3级地震［J］.国际地震动态，325(4):18—28.

CHEN Xue-zhong，JIANG Chang-sheng，LI Yan-e.2008.The March 21，2008 Yutian，Xinjiang MS7.3 earthquake［J］.Recent Development in World Seismology，325(4):18—28(in Chinese).

邓万明.1989.新疆乌鲁克库勒火山群地质考察:“1951年1号火山喷发”质疑［J］.自然资源学报，4(4):349—354.

DENG Wan-ming.1989.Geological investigation on the volcanic group in the Wulukule(lake)region，Xinjiang:Query the validity of“the eruption of No.1 volcano in 1951”［J］.Acta Natural Resources，4(4):349—354(in Chinese).

刘嘉麒，买买提依明.1990.西昆仑山第四纪火山的分布与K－Ar年龄［J］.中国科学(D辑)，(2):180—187.

LIU Jia-qi，YIMING Maimaiti.1990.The distribution of Quaternary volcanoes and its K － Ar isotopic age in western Kunlun Mts［J］.Science in China(Ser D)，(2):180—187(in Chinese).

刘若新.2000.中国的活火山［M］.北京:地震出版社.

LIU Ruo-xin.2000.Active Volcanoes in China［M］.Seismological Press，Beijing(in Chinese).

单新建，张国宏.2006.孕震区震前D－InSAR干涉形变场动态演化图像分析［J］.地震地质，28(3):441—446.

SHAN Xin-jian，ZHANG Guo-hong.2006.An analysis of dynamic evolution of preseismic interferometric deformation fields in seismic area［J］.Seismology and Geology，28(3):441—446(in Chinese).

孙建宝，梁芳，沈正康，等.2008.汶川MS8.0地震InSAR形变观测及初步分析［J］.地震地质，30(3):789—795.

SUN Jian-bao，LIANG Fang，SHEN Zheng-kang，et al.2011.InSAR deformation observation and preliminary analysis of the MS8 Wenchuan earthquake［J］.Seismology and Geology，30(3):789—795(in Chinese).

新疆日报.1951.新疆昆仑大阪西沟发现火山喷发［Z］.7月5日，第二版.

Xinjiang Daily.1951.Volcano Eruption Occurred in the Western Valley of Kunlun Mountains，Southern Xinjiang［Z］.May 7，page 2(in Chinese).

许建东，赵波，张柳毅，等.2011.新疆阿什库勒火山群野外地质科学考察［J］.地震地质，33(3):499—505.doi:10.3969/j.issn.0253 －4967.2011.03.001.

XU Jian-dong，ZHAO Bo，ZHANG Liu-yi，et al.2011.Field geological exploration of Ashikule volcano group in western Kunlun Mountains［J］.Seismology and Geology，33(3):499—505(in Chinese).

徐锡伟，谭锡斌，吴国栋，等.2011.2008年于田MS7.3地震地表破裂带特征及其构造属性讨论［J］.地震地质，33(2):462—471.doi:10.3969/j.issn.0253 －4967.2011.02.019.

XU Xi-wei，TAN Xi-bin，WU Guo-dong，et al.2011.Surface rupture features of the 2008 Yutian MS7.3 earthquake and its tectonic nature［J］.Seismology and Geology，33(2):462—471(in Chinese).

Berardino P，Fornaro G，Lanari R，et al.2002.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms［J］.Institute of Electrical and Electronics Engineers Transactions on Geosciences and Remote Sensing，40(11):2375—2383.

Ferretti A，Prati C，Rocca F.2000.Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry［J］.IEEE Transactions on Geoscience and Remote Sensing，38(5):2202—2212.doi:10.1109/36.868878.

Jung H S，Lee C W，Park J W，et al.2008.Improvement of small baseline subset(SBAS)algorithm for measuring time-series surface deformations from differential SAR interferograms［J］.Korean Journal of Remote Sensing，24:165—177.

Liu J Q，Maimaiti Y M.1989.Distribution and ages of Ashikule volcanoes on the West Kunlun Mountains，west China［J］.Bulletin of Glacier Research，7:187—190.

Lu Z.2007.InSAR imaging of volcanic deformation over cloud-prone areas:Aleutian Islands［J］.Photogrammetric Engineering ＆ Remote Sensing，73(3):245—257.

Massonnet D，Rossi M，Carmona C，et al.1993.The displacement field of the Landers earthquake mapped by radar interferometry［J］.Nature，364:138—142.

Massonnet D，Briole P，Arnaud A.1995.Deflation of Mount Etna monitored by spaceborne radar interferometry［J］.Nature，375:567—570.

Peltzer G，CrampéF，Hensley S，et al.2001.Transient strain accumulation and fault interaction in the eastern California shear zone［J］.Geology，29:975—978.

Rosen P，Hensley S，Zebker H，et al.1996.Surface deformation and coherence measurements of Kilauea Volcano，Hawaii，from SIR － C radar interferometry［J］.Journal of Geophysical Rearch，101:23109—23125.

Wright T J，Parsons B E，Lu Z.2004.Toward mapping surface deformation in three dimensions using InSAR［J］.Geophysical Research Letter，31(1):L01607.

Zebker H A，Rosen P A，Goldstein R M，et al.1994.On the derivation of coseismic displacement fields using differential radar interferometry—the Landers earthquake［J］.Journal of Geophysical Research，99(B10):19617—19634.

Zebker H A，Rosen P A，Hensley S.1997.Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps［J］.Journal of Geophysical Research，102:7547—7563.