2. 华东交通大学 土木建筑学院,江西 南昌 330013
2. School of Civil Engineering and Architecture,East China Jiaotong University,Nanchang 330013,China
1 引 言
汶川地震在龙门山地区引发的地质灾害在短期内难以消除，尤其是震后的5～10年，滑坡和泥石流将成为主要的高风险性灾害源[1, 2]。及时准确调查主断裂带地区的滑坡分布特征，可为震后滑坡灾害监测、斜坡稳定性评估以及灾后规划重建提供重要的基础观测数据。雷达干涉测量(interfer-ometric synthetic aperture radar,InSAR)在滑坡探测及地表形变监测中已初步展现出极大的应用潜力[3, 4]。然而，由于InSAR受雷达回波信号失相关和大气相位延迟等负面干扰，多种改进的干涉测量方法，如永久散射体、短基线集、角反射器等InSAR技术应运而生[5, 6, 7, 8, 9, 10, 11]，为区域地质灾害的时空演变分析提供丰富的形变观测资料。
在地质构造特征复杂的龙门山地区开展震后滑坡探测，除了受到SAR回波信号失相关影响外，还受到龙门山特殊地形地貌引起的干涉相位负面影响。一方面，由于龙门山具有剧烈起伏的地形特征，突变的高程变化容易导致干涉相位信号的整周跳变，使得从干涉相位中提取有用的形变信息更加困难；另一方面，过大的地形高差使得山区近地表的大气湿度存在显著差异，不同成像时间的大气变化导致SAR信号差异性的相位延迟[12, 13, 14, 15, 16, 17]，山区大气在时空尺度上20%的相对湿度变化，可导致10～14cm的形变监测误差。
针对InSAR在龙门山高山险峻地区监测滑坡灾害所面临的技术难题，本文从以下3方面开展了InSAR监测滑坡形变灾害的理论与方法探索。首先，为克服龙门山剧烈地形起伏对干涉相位信号数值的周跳影响，选取具有超短基线特征的两幅雷达影像构建干涉像对；其次，为克服山区近地表大气湿度时空分布不均匀引起的SAR相位信延迟问题，提出采用区域三维空间因子的大气相位建模方法，建立区域自适应大气修正模型以分离长波量级的大气相位干扰；第三，采用雷达视线方向形变阈值扣除大气短波扰动和噪声相位影响，联合DEM地形坡度因子和形变阈值解译滑坡灾害的空间分布特征。
2 短基线InSAR探测滑坡的计算模型与方法 2.1 干涉相位信号分量解析对重复轨道星载SAR影像进行干涉测量，各像元的干涉相位由以下5个相位分量组成[6, 7, 8, 9]
式中，Φflat为参考椭球面相位；Φtop表示地形起伏相位；Φdef为地表形变相位；Φatm表示大气延迟相位；Φnoi代表随机噪声相位。
借助卫星轨道姿态数据和干涉几何关系，可从干涉相位中计算并去除参考相位分量Φflat，利用已有的DEM可以部分去除地形相位Φtop分量(由于所用的DEM可能存在高程误差)，采用滤波方法减弱随机噪声Φnoi，则各像素单元的差分干涉相位Φdint可表示为如下3个主要贡献分量[9, 10, 11]
式中，δφres_top表示由于DEM高程误差δh引入的地形残留相位；λ为雷达波波长；B⊥为垂直空间基线长度；R为雷达天线至地面目标的斜距；θ为雷达波入射角；Δr为雷达视线方向的位移量。
地形残留相位δφres_top与干涉对的有效空间基线长度B⊥(即垂直基线分量)成正比关系[9]。在龙门山震后滑坡形变监测中，为尽量减弱剧烈地形起伏和DEM高程误差可能引入的残留相位分量，选取具有超短基线特征的干涉像对开展滑坡形变监测，例如65m基线长度的干涉像对中20m高程误差引入的相位误差约为0.13rad，相对于滑坡形变相位数值来说，残留地形相位贡献微小，因此超短基线InSAR有利于对形变信号的准确提取。
2.2 基于三维空间因子的山区大气相位建模SAR电磁波信号在大气层传播过程中，主要受到对流层湿分量延迟影响[12, 13, 14]。对于不同成像时刻t1和t2，设δqpt1和δqpt2为地面上不同高度p、q两点之间电磁波传播路径延迟，则对应的大气相位延迟差异为
式中，λ为雷达波长；θ为雷达波入射角；Φp、Φq分别表示p、q两点的大气延迟相位。
对流层大气折射指标N可以表示为三维空间因子(x,y,z)的函数，是水平方向大气紊流δN(x,y,z)以及温度、湿度和气压等因素垂直分层N(z)两种物理过程的共同作用[12, 15, 16, 17]。将对流层垂直分成若干薄层，则雷达信号传播路径延迟可由不同时序大气折射指标Nti计算
代入式(3)，得到干涉相位中的大气延迟分量
龙门山地区缺乏与SAR成像时间同步的高分辨率地面气象数据，对上式折射指标的量化估算相当困难。值得注意的是，在地形起伏较大的复杂山区，大气垂直分层效应尤其明显，大气附加相位与地形起伏具有较强的相关性[12, 15]。为削弱大气相位延迟误差，国内外学者建立了以地形高程因子h为自变量的大气相位修正模型，其线性模型和非线性模型分别为
式中，Φatm为大气相位值；a1、a2和b1、b2、…、b5分别为线性和非线性模型系数。
顾及山区不同平面区域大气的不均匀特性，本文在上述地形高程因子h的基础上，提出增加顾及平面位置因子(x,y)的大气建模方法，建立以三维空间因子(x,y,h)为自变量的大气相位函数模型
式中，(x,y)为像素点的雷达方位向和距离向坐标，表示三维空间因子中的平面位置；h为像素点的地形高程；c1、c2、…、c7为模型参数。将解缠干涉图中以大气延迟为主分量的相位值视为观测量，在解缠干涉图中提取均匀分布的相位值，联合像元坐标因子(x,y,h)代入式(8)，采用最小二乘方法解算出模型系数ci，获得大气相位修正模型。 2.3 联合InSAR形变和DEM坡度因子探测滑坡
通过上述的三维地理空间因子大气建模，可以计算出雷达干涉图中的大气相位，从短基线解缠干涉图中减去大气相位贡献分量，获得地表沿雷达视线方向的位移量r。为提高滑坡探测结果的可靠性，从干涉图中提取高信噪比像元的地表位移量r进行统计分析，计算形变标准偏差mr。滑坡形变显著大于该地区平均地表形变r，以形变偏差大于3倍标准差的单元为滑坡探测位置，以此剔除大气局部短波相位扰动。设形变阈值Td=r+3mr，提取InSAR形变观测量大于该阈值Td的区域为滑坡初选位置。
已有实测滑坡资料表明[2]，龙门山地区滑坡多发生于斜坡超过10°的区域，为消除上述单一形变阈值可能引入的滑坡位置误判，进一步引入地形坡度因子进行滑坡区域的综合识别[18, 19]。地面坡度因子S是地形曲面z=f(x,y)在东西向和南北向高程变化率的函数，坡度因子可根据DEM数据计算得到，即
式中，fx是南北方向高程变化率；fy是东西方向高程变化率。坡度因子可采用局部范围数值微分方法进行计算[18, 19]，如图 1所示的3×3移动窗口，本文采用如下3阶不带权差分计算模型 式中，zi(i=1,2,…,8)为中心点周围格网点的高程；d为DEM格网分辨率。
|
| 图 1 DEM3×3移动计算窗口 Fig. 1 Gradient calculation window |
设定某一坡度阈值Ts，对InSAR初选滑坡区域进一步识别，在初选滑坡中判定坡度因子大于给定阈值Ts的地面目标为滑坡位置，以此消除在平坦地区可能引入的滑坡误判信息。
3 研究区域与数据分析根据上述InSAR探测滑坡的计算模型与方法，开展汶川震后龙门山主断裂两侧降雨期滑坡分布的调查研究。选取垂直基线长度为65m的ALOS/PALSAR影像建立超短基线干涉对，两幅影像成像时间分别为2010-07-27和2010-09-11，时间间隔为46d，地震高烈度地区的映秀、绵竹和北川于2010年7月至9月遭受了多次强降雨过程，期间的8月12日至14日为强降雨期，诱发了严重的滑坡和泥石流地质灾害。
PALSAR影像的斜距向和方位向分辨率分别为4.68m和3.15m，两幅影像重叠区域的中心经纬度为(104.04°E，31.44°N)，影像的地面覆盖范围约为70km×70km(见图 2)。影像覆盖了龙门山西北方向的汶川和茂县部分地区，以及东南方向的德阳和什邡等区域，该区域呈现西北高而东南低的地形特征。
|
| F1:汶川-茂县断裂;F2:映秀-北川断裂;F3:灌县-江油断裂;★:主震震中;●:余震序列 图 2 ALOS/PALSAR影像覆盖区域 Fig. 2 Ground coverage of ALOS/PALSAR images |
对上述两幅PALSAR影像进行干涉处理，图 3为去除参考相位后的干涉相位图，主要为地形起伏、地表形变和大气效应的综合体现。在65m超短基线长度条件下，干涉图条纹的模糊度高为980m，条纹对地形高程变化的敏感度非常低，有利于大气信息和形变信息的准确提取。采用美国NASA Terra卫星ASTER传感器获得的DEM去除地形相位[20]，图 4为研究区域的ASTER DEM。该区域西北地区为龙门山推覆构造带，群山起伏，海拔高程为1000～5000m，东南地区为成都平原，海拔高程为480～1000m。图 5为解缠后的差分干涉相位图，主要为地表形变和大气相位的综合体现[21, 22]。
|
| 图 3 缠绕的干涉相位图 Fig. 3 Wrapped phase |
|
| 图 4 影像覆盖区的ASTER DEM Fig. 4 ASTER DEM over the study area |
|
| 图 5 解缠相位图 Fig. 5 Unwrapped phase |
从图 4和图 5可以看出，解缠干涉相位与DEM呈现出较强的空间相关性，表明解缠相位中的主分量为受地形高程和平面位置影响显著的大气延迟相位。进一步选取图 5中的AB剖面作为定量分析单元，提取高信噪比像元(相干系数大于0.3)的解缠相位和对应像素的DEM高程进行分析。如图 6所示，左侧坐标纵轴为解缠相位数据，右侧坐标纵轴为DEM高程信息，两条曲线变化规律呈现出较高的一致性，表明解缠相位与地形高程表现出较强的空间相关性。
|
| 图 6 解缠干涉相位和DEM高程数值序列 Fig. 6 Relationship between unwrapped phase and topographic elevation |
对上述获得的解缠干涉相位数值进行区域大气相位建模。以单一地形高程因子h为自变量的线性和非线性大气建模结果，分别如图 7中的(a)和(b)所示，其中红色离散点为InSAR解缠相位，图 7(a)蓝色直线为线性建模拟合结果，建模精度(即建模后的相位均方根误差RMSE)为0.67rad；图 7(b)为非线性建模结果(蓝色曲线)，建模精度为0.66rad，线性与非线性建模结果无显著差异。由于没有考虑山区大气水平状态分布不均匀的特性，这两种建模结果与InSAR实际观测相位(解缠相位)差异性较大。
|
| 图 7 大气相位空间建模 Fig. 7 Modeling of atmospheric phase effect |
根据前述提出的顾及大气相位在水平方向的扰动影响，在顾及高程因子的基础上，增加平面区域二维位置因子(x,y)对该地区大气相位进行建模，计算结果如图 7(c)所示，图中蓝色离散点为建模计算相位。与图 7(a)和图 7(b)进行比较，图 7(c)中的建模相位为高程因子的多值函数，具有相同高程但位于不同平面区域的大气可能具有不同的相位延迟。从图 7(c)中可以看出，基于三维空间因子的建模相位与InSAR实测大气相位吻合性较好。统计表明，顾及三维空间因子大气相位建模的精度(RMSE)达到0.45rad，相比于单一高程因子的建模方法，相位建模精度提高了约32%。图 7(d)为基于三维空间因子(x,y,h)模型所得到的区域大气相位图。
3.3 滑坡空间分布特征从差分干涉图中分离和去除大气相位数值后，进一步提取地表形变，设相关系数阈值为0.3，提取SAR影像中相干系数高于该阈值的像元进行统计计算，提取形变量超过形变阈值的像元区域为滑坡初选区域，再依据各像素的ASTER DEM坡度因子，进一步判别坡度因子大于阈值Ts=10°的地表区域为滑坡位置，由此提取出该区域滑坡群的空间展布特征，如图 8所示。
|
| 图 8 探测出的滑坡区域及其空间分布特征 Fig. 8 Spatial distribution pattern of the detected landslides |
在图 8的中间图幅中，背景图为PALSAR的振幅影像，图中的黄色离散区域为探测得到的滑坡群。从图中滑坡的空间分布来看，这些滑坡群密集分布于主断裂带即映秀-北川断裂带(黑色虚线)的两侧，沿NE-SW方向呈条带状展布，距离发震断层约3～15km，其中85%的滑坡位于主断裂带两侧10km范围内。图中所示滑坡的数量和面积在断裂带两侧呈现出明显的不对称性，其中70%的滑坡位于逆冲断层的上盘，而下盘滑坡仅占约30%，主要原因可能在于汶川地震的发震断层破裂方式为单侧破裂[1, 6, 25]，上盘逆冲的近场同震地表位移达到约9.0m，导致震后滑坡受发震断层上盘逆冲特性的滞后效应较为显著。
图 8中间图幅内的红色三角形标识为面积较大的滑坡群，编号分别为(A-L)，图 8中的四周图幅则分别表示这些滑坡群的详细位置及其沿雷达LOS方向的形变量，其背景图为DEM坡度因子，其中的多边形区域为滑坡灾害位置，这些大型滑坡中约有2/3的滑坡分布于断层上盘，仍然呈现出显著的上盘逆冲后效应现象。
龙门山主断裂带两侧产生密集的滑坡分布，一方面受汶川地震导致的龙门山部分斜坡岩体破碎，多数山体产生了大量松散固体物质，山体整体性遭到反复破坏，大量斜坡体处于高位、高危和高发的不稳定状态；另一方面，由于震后地质环境更加脆弱，降雨诱发滑坡的敏感性极高。在2010年8月12日至14日，主断裂带地区遭遇强降雨天气[23]，灾区3天总降水量达229.6mm，最大小时雨强为32.2mm，为震后两年来的极强降雨，由此诱发大量的山体滑坡。
龙门山震后降雨期滑坡灾害的空间分布特征呈现出上述显著的上盘效应、主断裂带距离效应，此外，进一步的统计分析还表明，震后滑坡还与局部地形高程、坡度因子等因素密切相关。对上述活动滑坡的地形高程和坡度因子进行数据统计，约有90%的滑坡群分布于1000～3000m的中低海拔山区，坡度因子则多集中于15°～35°的中小坡度范围，而高海拔地区和陡峭区域的滑坡分布相对较少，表明龙门山震后滑坡与区域地形的高程和坡度因子也存在较强的空间联系。
4 结 论为探测龙门山地区汶川震后强降雨诱发大面积滑坡的空间分布特征，本文引入超短基线InSAR技术开展主断裂带两侧的震后滑坡探测，提出基于三维空间因子的InSAR大气相位修正模型去除大气长波量级的信号干扰，采用地表形变阈值扣除大气短波相位扰动，结合地形坡度因子，综合提取震后强降雨期滑坡分布的空间特征。
超短基线InSAR计算结果表明，龙门山主断裂带两侧的震后滑坡灾害基本沿发震断层北川-映秀断裂带呈条带状密集分布，距离发震断层约为3～15km，滑坡密度与数量呈现出显著的上盘效应，其中70%滑坡分布于逆冲断层的上盘区，而下盘滑坡仅占约30%，此外90%滑坡分布于高程1000～3000m的中低海拔山区和坡度因子为15°～35°的中小坡度范围。上述特征表明，震后滑坡一方面受局部地形高程、坡度因子和强降雨因素所激发；另一方面，震后滑坡还与汶川地震发震断层的逆冲特性相关，使得主断裂带两侧的震后滑坡受上盘逆冲运动的滞后效应较为显著。
| [1] | ZHANG Peizhen, XU Xiwei, WEN Xueze, et al. Slip Rates and Recurrence Intervals of Longmen Shan Active Fault Zone, and Tectonic Implications for the Mechanism of the May 12 Wenchuan Earthquake, 2008, Sichuan, China [J]. Chinese Journal of Geophysics, 2008, 51(4): 1066-1073. (张培震, 徐锡伟, 闻学泽, 等. 2008年汶川8.0级地震发震断裂的滑动速率、复发周期和构造成因[J]. 地球物理学报, 2008, 51(4): 1066-1073.) |
| [2] | HUANG Runqiu. After Effect of Geohazards Induced by the Wenchuan Earthquake [J]. Journal of Engineering Geology, 2011, 19(2): 145-150. (黄润秋. 汶川地震地质灾害后效应分析[J]. 工程地质学报, 2011, 19(2): 145-150.) |
| [3] | XU Caujun, JIANG Guoyan, WANG Jianjun, et al. A Seismic Hazard Evaluating System for Active Faults Based on GNSS/InSAR/GIS [J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 661-669. (许才军, 江国焰, 汪建军, 等. 基于GNSS/InSAR/GIS的活动断层地震危险性评估系统[J]. 测绘学报, 2012, 41(5): 661-669.) |
| [4] | LONG Sichun, ZHANG Shiyu, FENG Tao, et al. A New Approach of Weighted Stacking Based on Common Master Image and Its Application in Ground Subsidence Monitoring [J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(6): 844-850. (龙四春, 张诗玉, 冯涛, 等. 公用主影像干涉图加权叠加方法及其在地面沉降监测中的应用[J]. 测绘学报, 2012, 41(6): 844-850.) |
| [5] | WANG Yan, LIAO Mingsheng, LI Deren, et al. Subsidence Velocity Retrieval from Long-term Coherent Targets in Radar Interferometric Stacks[J]. Chinese Journal of Geophysics, 2007, 50(2): 598-604. (王艳, 廖明生, 李德仁, 等. 利用长时间序列相干目标获取地面沉降场[J]. 地球物理学报, 2007, 50(2): 598-604.) |
| [6] | CHEN Qiang, LIU Guoxiang, LI Yongshu, et al. Automated Detection of Permanent Scatter in Radar Interferometry: Algorithm and Testing Results[J]. Acta Geodeatica et Cartographica Sinica, 2006, 35(2): 112-117. (陈强, 刘国祥, 李永树, 等. 干涉雷达永久散射体自动探测:算法与实验结果[J]. 测绘学报, 2006, 35(2): 112-117.) |
| [7] | ZHANG Qin, ZHAO Chaoying, DING Xiaoli, et al. Research on Recent Characteristics of Spatio-temporal Evolution and Mechanism of Xi’an Land Subsidence and Ground Fissure by Using GPS and InSAR Techniques [J]. Chinese Journal of Geophysics, 2009, 52(5): 1214-1222. (张勤, 赵超英, 丁晓利, 等. 利用GPS与InSAR研究西安现今地面沉降与地裂缝时空演化特征[J]. 地球物理学报, 2009, 52 (5): 1214-1222.) |
| [8] | CHEN Qiang, DING Xiaoli, LIU Guoxiang, et al. Baseline Recognition and Parameter Estimation of Persistent-scatterer Network in Radar Interferometry[J]. Chinese Journal of Geophysics, 2009, 52(9): 2229-2236. (陈强, 丁晓利, 刘国祥, 等. 雷达干涉PS网络的基线识别与解算方法[J]. 地球物理学报, 2009, 52(9): 2229-2236.) |
| [9] | TANG Yixian, ZHANG Hong, WANG Chao. Subsidence Retrieval in Suzhou City Using Permanent Scatterers(PS) InSAR [J]. Progress in Natural Science, 2006, 16(8): 97-102. (汤益先, 张红, 王超. 基于永久散射体雷达干涉测量的苏州地区沉降研究[J]. 自然科学进展, 2006, 16(8): 97-102.) |
| [10] | BERARDINO P, FORNARO G, LANARI R, et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40: 2375-2383. |
| [11] | FERRETTI A, PRATI C, ROCCA F. Permanent Scatterers in SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39: 8-20. |
| [12] | LI Z W, DING X L, HUANG C, et al. Modeling of Atmospheric Effects on InSAR Measurements by Incorporating Terrain Elevation Information[J]. Journal of Atmospheric Solar-Terrestrial Physics, 2006, 68(11): 1189-1194. |
| [13] | LI Z W, DING X L, LIU G X. Modeling Atmospheric Effects on InSAR with Meteorological and Continuous GPS Observations: Algorithm and Some Test Results[J]. Journal of Atmospheric and Solar-Terrestrial Physics, 2004, 66(11): 907-917. |
| [14] | ONN F, ZEBKER H A. Correction for Interferometric Synthetic Aperture Radar Atmospheric Phase Artifacts Using Time Series of Zenith Wet Delay Observations from a GPS Network [J]. Journal of Geophysical Research, 2006, 111( B9): 102-110. |
| [15] | LI Z H, FIELDING E J, CROSS P, et al. Interferometric Synthetic Aperture Radar Atmospheric Correction: GPS Topography-dependent Turbulence Model[J]. Journal of Geophysical Research, 2006, 111(B2): 224-237. |
| [16] | DING X L, LI Z W, ZHU J J, et al. Atmospheric Effects on InSAR Measurements and Their Mitigation [J]. Sensors, 2008(8): 5426-5448. |
| [17] | REMY D, BONVALOT S, BRIOLE P, et al. Accurate Measurements of Tropospheric Effects in Volcanic Areas from SAR Interferometry Data: Application to Sakurajima Volcano (Japan) [J]. Earth and Planetary Science Letters, 2003, 213(3): 299-310. |
| [18] | ZHAO C Y, LU Z, ZHANG Q, et al. Large-area Landslide Detection and Monitoring with ALOS/PALSAR Imagery Data over Northen California and Southen Oregon, USA [J]. Remote Sensing of Environment, 2012 (124): 348-359. |
| [19] | LIU Xuejun, GONG Jianya, ZHOU Qiming, et al. A Study of Accuracy and Algorithms for Calculating Slope and Aspect Based on Grid Digital Elevation Model (DEM) [J]. Acta Geodaetica et Cartographica Sinica, 2004, 33(3): 258-263. (刘学军, 龚健雅, 周启鸣, 等. 基于DEM坡度坡向算法精度的分析研究[J]. 测绘学报, 2004, 33(3): 258-263.) |
| [20] | NICKI F S, HAROLD G, PETER J M. NASA EOS Terra ASTER: Volcanic Topographic Mapping and Capability[J]. Remote Sensing of Environment, 2004, 90(3): 405-414. |
| [21] | CHEN C W, ZEBKER H A. Phase Unwrapping for Large SAR Interferograms: Statistical Segmentation and Generalized Network Models[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1709-1719. |
| [22] | CHEN Qiang, YANG Yinghui, LIU Guoxiang, et al. InSAR Phase Unwrapping Using Least Squares Method with Integer Ambiguity Resolution and Edge Detection[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 441-448. (陈强, 杨莹辉, 刘国祥, 等. 基于边界探测的InSAR最小二乘整周相位解缠方法[J]. 测绘学报, 2012, 41(3): 441-448.) |
| [23] | TANG Chuan, LI Weile, DING Jun, et al. Field Investigation and Research on Giant Debris Flow on August 14, 2010 in Yingxiu Town, Epicenter of Wenchuan Earthquake[J]. Earth Science: Journal of China University of Geosciences, 2011, 36(1): 172-180. (唐川, 李为乐, 丁军, 等. 汶川震区映秀镇“8.14”特大泥石流灾害调查[J]. 地球科学: 中国地质大学学报, 2011, 36(1):172-180.) |